[Gmp-commit] /var/hg/gmp: 3 new changesets

[mercurial at gmplib.org]

details:   /var/hg/gmp/rev/cc2a98d41eab
changeset: 16309:cc2a98d41eab
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Fri Feb 21 19:28:25 2014 +0100
description:
mini-gmp: New function mpn_sqrtrem.

details:   /var/hg/gmp/rev/5cbbb43c9217
changeset: 16310:5cbbb43c9217
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Fri Feb 21 19:30:14 2014 +0100
description:
mini-gmp/tests/t-sqrt.c: Test mpn_sqrtrem.

details:   /var/hg/gmp/rev/45bd52435e58
changeset: 16311:45bd52435e58
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Fri Feb 21 19:32:38 2014 +0100
description:
ChangeLog

diffstat:

 ChangeLog               |   6 ++++
 mini-gmp/mini-gmp.c     |  21 ++++++++++++++
 mini-gmp/mini-gmp.h     |   1 +
 mini-gmp/tests/t-sqrt.c |  70 ++++++++++++++++++++++++++++++++++++++++++++----
 4 files changed, 92 insertions(+), 6 deletions(-)

diffs (170 lines):

diff -r 4560d88e082a -r 45bd52435e58 ChangeLog
--- a/ChangeLog	Mon Feb 17 22:13:07 2014 +0100
+++ b/ChangeLog	Fri Feb 21 19:32:38 2014 +0100
@@ -1,3 +1,9 @@
+2014-02-21 Marco Bodrato <bodrato at mail.dm.unipi.it>
+
+	* mini-gmp/mini-gmp.c (mpn_sqrtrem): New function.
+	* mini-gmp/mini-gmp.h: Declare it.
+	* mini-gmp/tests/t-sqrt.c: Test it.
+
 2014-02-17  Niels Möller  <nisse at lysator.liu.se>
 
 	* mpn/generic/div_qr_1.c (mpn_div_qr_1): Revert yesterday's fix.
diff -r 4560d88e082a -r 45bd52435e58 mini-gmp/mini-gmp.c
--- a/mini-gmp/mini-gmp.c	Mon Feb 17 22:13:07 2014 +0100
+++ b/mini-gmp/mini-gmp.c	Fri Feb 21 19:32:38 2014 +0100
@@ -3285,6 +3285,27 @@
   assert (p [n-1] != 0);
   return mpz_root (NULL, mpz_roinit_n (t, p, n), 2);
 }
+
+mp_size_t
+mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n)
+{
+  mpz_t s, r, u;
+  mp_size_t res;
+
+  assert (n > 0);
+  assert (p [n-1] != 0);
+
+  mpz_init (r);
+  mpz_init (s);
+  mpz_rootrem (s, r, mpz_roinit_n (u, p, n), 2);
+  mpn_copyd (sp, s->_mp_d, s->_mp_size);
+  mpz_clear (s);
+  res = r->_mp_size;
+  if (rp)
+    mpn_copyd (rp, r->_mp_d, res);
+  mpz_clear (r);
+  return res;
+}
 

 /* Combinatorics */
 
diff -r 4560d88e082a -r 45bd52435e58 mini-gmp/mini-gmp.h
--- a/mini-gmp/mini-gmp.h	Mon Feb 17 22:13:07 2014 +0100
+++ b/mini-gmp/mini-gmp.h	Fri Feb 21 19:32:38 2014 +0100
@@ -99,6 +99,7 @@
 void mpn_mul_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
 void mpn_sqr (mp_ptr, mp_srcptr, mp_size_t);
 int mpn_perfect_square_p (mp_srcptr, mp_size_t);
+mp_size_t mpn_sqrtrem (mp_ptr, mp_ptr, mp_srcptr, mp_size_t);
 
 mp_limb_t mpn_lshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int);
 mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int);
diff -r 4560d88e082a -r 45bd52435e58 mini-gmp/tests/t-sqrt.c
--- a/mini-gmp/tests/t-sqrt.c	Mon Feb 17 22:13:07 2014 +0100
+++ b/mini-gmp/tests/t-sqrt.c	Fri Feb 21 19:32:38 2014 +0100
@@ -24,7 +24,7 @@
 #include "testutils.h"
 
 #define MAXBITS 400
-#define COUNT 10000
+#define COUNT 9000
 
 static void
 dump (const char *label, const mpz_t x)
@@ -61,24 +61,71 @@
 }
 
 void
+mpz_mpn_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
+{
+  mp_limb_t *sp, *rp;
+  mp_size_t un, sn, ret;
+
+  un = mpz_size (u);
+
+  mpz_xor (s, s, u);
+  sn = (un + 1) / 2;
+  sp = mpz_limbs_write (s, sn + 1);
+  sp [sn] = 11;
+
+  if (un & 1)
+    rp = NULL; /* Exploits the fact that r already is correct. */
+  else {
+    mpz_add (r, u, s);
+    rp = mpz_limbs_write (r, un + 1);
+    rp [un] = 19;
+  }
+
+  ret = mpn_sqrtrem (sp, rp, mpz_limbs_read (u), un);
+
+  if (sp [sn] != 11)
+    {
+      fprintf (stderr, "mpn_sqrtrem buffer overrun on sp.\n");
+      abort ();
+    }
+  if (un & 1) {
+    if ((ret != 0) != (mpz_size (r) != 0)) {
+      fprintf (stderr, "mpn_sqrtrem wrong return value with NULL.\n");
+      abort ();
+    }
+  } else {
+    mpz_limbs_finish (r, ret);
+    if (ret != mpz_size (r)) {
+      fprintf (stderr, "mpn_sqrtrem wrong return value.\n");
+      abort ();
+    }
+    if (rp [un] != 19)
+      {
+	fprintf (stderr, "mpn_sqrtrem buffer overrun on rp.\n");
+	abort ();
+      }
+  } 
+  
+  mpz_limbs_finish (s, (un + 1) / 2);
+}
+
+void
 testmain (int argc, char **argv)
 {
   unsigned i;
   mpz_t u, s, r;
 
-  mpz_init (u);
   mpz_init (s);
   mpz_init (r);
 
-  mpz_set_si (u, -1);
+  mpz_init_set_si (u, -1);
   if (mpz_perfect_square_p (u))
     {
       fprintf (stderr, "mpz_perfect_square_p failed on -1.\n");
       abort ();
     }
 
-  mpz_set_ui (u, 0);
-  if (!mpz_perfect_square_p (u))
+  if (!mpz_perfect_square_p (s))
     {
       fprintf (stderr, "mpz_perfect_square_p failed on 0.\n");
       abort ();
@@ -86,7 +133,7 @@
 
   for (i = 0; i < COUNT; i++)
     {
-      mini_rrandomb (u, MAXBITS);
+      mini_rrandomb (u, MAXBITS - (i & 0xFF));
       mpz_sqrtrem (s, r, u);
 
       if (!sqrtrem_valid_p (u, s, r))
@@ -98,6 +145,17 @@
 	  abort ();
 	}
 
+      mpz_mpn_sqrtrem (s, r, u);
+
+      if (!sqrtrem_valid_p (u, s, r))
+	{
+	  fprintf (stderr, "mpn_sqrtrem failed:\n");
+	  dump ("u", u);
+	  dump ("sqrt", s);
+	  dump ("rem", r);
+	  abort ();
+	}
+
       if (mpz_sgn (r) == 0) {
 	mpz_neg (u, u);
 	mpz_sub_ui (u, u, 1);