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

mercurial at gmplib.org mercurial at gmplib.org
Wed Mar 8 18:57:49 UTC 2017


details:   /var/hg/gmp/rev/6e859993b5b4
changeset: 17324:6e859993b5b4
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Wed Mar 08 19:13:47 2017 +0100
description:
tests/devel/sqrtrem_1_2.c (main): Interface improvement.

details:   /var/hg/gmp/rev/bba7e03de01a
changeset: 17325:bba7e03de01a
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Wed Mar 08 19:14:48 2017 +0100
description:
mpn/generic/sqrtrem.c: Direct use of sqrtrem2 when n==2.

details:   /var/hg/gmp/rev/a5fb41d654ff
changeset: 17326:a5fb41d654ff
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Wed Mar 08 19:30:42 2017 +0100
description:
mpn/generic/sqrtrem.c: Indent and comment.

details:   /var/hg/gmp/rev/6393a907d1b5
changeset: 17327:6393a907d1b5
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Wed Mar 08 19:57:44 2017 +0100
description:
mpn/generic/sqrtrem.c: Support SQRTREM2_INPLACE.

diffstat:

 mpn/generic/sqrtrem.c     |  105 +++++++++++++++++++++++++++++++--------------
 tests/devel/sqrtrem_1_2.c |   25 +++++++++-
 2 files changed, 92 insertions(+), 38 deletions(-)

diffs (210 lines):

diff -r f3cce34d1228 -r 6393a907d1b5 mpn/generic/sqrtrem.c
--- a/mpn/generic/sqrtrem.c	Tue Mar 07 16:55:42 2017 +0100
+++ b/mpn/generic/sqrtrem.c	Wed Mar 08 19:57:44 2017 +0100
@@ -166,12 +166,24 @@
 
 
 #define Prec (GMP_NUMB_BITS >> 1)
+#if ! defined(SQRTREM2_INPLACE)
+#define SQRTREM2_INPLACE 0
+#endif
 
 /* same as mpn_sqrtrem, but for size=2 and {np, 2} normalized
    return cc such that {np, 2} = sp[0]^2 + cc*2^GMP_NUMB_BITS + rp[0] */
+#if SQRTREM2_INPLACE
+#define CALL_SQRTREM2_INPLACE(sp,rp) mpn_sqrtrem2 (sp, rp)
+static mp_limb_t
+mpn_sqrtrem2 (mp_ptr sp, mp_ptr rp)
+{
+  mp_srcptr np = rp;
+#else
+#define CALL_SQRTREM2_INPLACE(sp,rp) mpn_sqrtrem2 (sp, rp, rp)
 static mp_limb_t
 mpn_sqrtrem2 (mp_ptr sp, mp_ptr rp, mp_srcptr np)
 {
+#endif
   mp_limb_t q, u, np0, sp0, rp0, q2;
   int cc;
 
@@ -222,44 +234,43 @@
   int c, b;			/* carry out of remainder */
   mp_size_t l, h;
 
+  ASSERT (n > 1);
   ASSERT (np[2 * n - 1] >= GMP_NUMB_HIGHBIT / 2);
 
-  if (n == 1)
-    c = mpn_sqrtrem2 (sp, np, np);
+  l = n / 2;
+  h = n - l;
+  if (h == 1)
+    q = CALL_SQRTREM2_INPLACE (sp + l, np + 2 * l);
   else
+    q = mpn_dc_sqrtrem (sp + l, np + 2 * l, h, 0, scratch);
+  if (q != 0)
+    ASSERT_CARRY (mpn_sub_n (np + 2 * l, np + 2 * l, sp + l, h));
+  TRACE(printf("tdiv_qr(,,,,%u,,%u) -> %u\n", (unsigned) n, (unsigned) h, (unsigned) (n - h + 1)));
+  mpn_tdiv_qr (scratch, np + l, 0, np + l, n, sp + l, h);
+  q += scratch[l];
+  c = scratch[0] & 1;
+  mpn_rshift (sp, scratch, l, 1);
+  sp[l - 1] |= (q << (GMP_NUMB_BITS - 1)) & GMP_NUMB_MASK;
+  if (UNLIKELY ((sp[0] & approx) != 0)) /* (sp[0] & mask) > 1 */
+    return 1; /* Remainder is non-zero */
+  q >>= 1;
+  if (c != 0)
+    c = mpn_add_n (np + l, np + l, sp + l, h);
+  TRACE(printf("sqr(,,%u)\n", (unsigned) l));
+  mpn_sqr (np + n, sp, l);
+  b = q + mpn_sub_n (np, np, np + n, 2 * l);
+  c -= (l == h) ? b : mpn_sub_1 (np + 2 * l, np + 2 * l, 1, (mp_limb_t) b);
+
+  if (c < 0)
     {
-      l = n / 2;
-      h = n - l;
-      q = mpn_dc_sqrtrem (sp + l, np + 2 * l, h, 0, scratch);
-      if (q != 0)
-	ASSERT_CARRY (mpn_sub_n (np + 2 * l, np + 2 * l, sp + l, h));
-      TRACE(printf("tdiv_qr(,,,,%u,,%u) -> %u\n", (unsigned) n, (unsigned) h, (unsigned) (n - h + 1)));
-      mpn_tdiv_qr (scratch, np + l, 0, np + l, n, sp + l, h);
-      q += scratch[l];
-      c = scratch[0] & 1;
-      mpn_rshift (sp, scratch, l, 1);
-      sp[l - 1] |= (q << (GMP_NUMB_BITS - 1)) & GMP_NUMB_MASK;
-      if (UNLIKELY ((sp[0] & approx) != 0)) /* (sp[0] & mask) > 1 */
-	return 1; /* Remainder is non-zero */
-      q >>= 1;
-      if (c != 0)
-	c = mpn_add_n (np + l, np + l, sp + l, h);
-      TRACE(printf("sqr(,,%u)\n", (unsigned) l));
-      mpn_sqr (np + n, sp, l);
-      b = q + mpn_sub_n (np, np, np + n, 2 * l);
-      c -= (l == h) ? b : mpn_sub_1 (np + 2 * l, np + 2 * l, 1, (mp_limb_t) b);
-
-      if (c < 0)
-	{
-	  q = mpn_add_1 (sp + l, sp + l, h, q);
+      q = mpn_add_1 (sp + l, sp + l, h, q);
 #if HAVE_NATIVE_mpn_addlsh1_n_ip1 || HAVE_NATIVE_mpn_addlsh1_n
-	  c += mpn_addlsh1_n_ip1 (np, sp, n) + 2 * q;
+      c += mpn_addlsh1_n_ip1 (np, sp, n) + 2 * q;
 #else
-	  c += mpn_addmul_1 (np, sp, n, CNST_LIMB(2)) + 2 * q;
+      c += mpn_addmul_1 (np, sp, n, CNST_LIMB(2)) + 2 * q;
 #endif
-	  c -= mpn_sub_1 (np, np, n, CNST_LIMB(1));
-	  q -= mpn_sub_1 (sp, sp, n, CNST_LIMB(1));
-	}
+      c -= mpn_sub_1 (np, np, n, CNST_LIMB(1));
+      q -= mpn_sub_1 (sp, sp, n, CNST_LIMB(1));
     }
 
   return c;
@@ -419,7 +430,7 @@
 mp_size_t
 mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr np, mp_size_t nn)
 {
-  mp_limb_t *tp, s0[1], cc, high, rl;
+  mp_limb_t cc, high, rl;
   int c;
   mp_size_t rn, tn;
   TMP_DECL;
@@ -458,6 +469,32 @@
       }
     return rl != 0;
   }
+  if (nn == 2) {
+    mp_limb_t tp [2];
+    if (rp == NULL) rp = tp;
+    if (c == 0)
+      {
+#if SQRTREM2_INPLACE
+	rp[1] = high;
+	rp[0] = np[0];
+	cc = CALL_SQRTREM2_INPLACE (sp, rp);
+#else
+	cc = mpn_sqrtrem2 (sp, rp, np);
+#endif
+	rp[1] = cc;
+	return ((rp[0] | cc) != 0) + cc;
+      }
+    else
+      {
+	rl = np[0];
+	rp[1] = (high << (2*c)) | (rl >> (GMP_NUMB_BITS - 2*c));
+	rp[0] = rl << (2*c);
+	CALL_SQRTREM2_INPLACE (sp, rp);
+	cc = sp[0] >>= c;	/* c != 0, the higest bit of the root cc is 0. */
+	rp[0] = rl -= cc*cc;	/* Computed modulo 2^GMP_LIMB_BITS, because it's smaller. */
+	return rl != 0;
+      }
+  }
   tn = (nn + 1) / 2; /* 2*tn is the smallest even integer >= nn */
 
   if ((rp == NULL) && (nn > 8))
@@ -465,8 +502,8 @@
   TMP_MARK;
   if (((nn & 1) | c) != 0)
     {
-      mp_limb_t mask;
-      mp_ptr scratch;
+      mp_limb_t s0[1], mask;
+      mp_ptr tp, scratch;
       TMP_ALLOC_LIMBS_2 (tp, 2 * tn, scratch, tn / 2 + 1);
       tp[0] = 0;	     /* needed only when 2*tn > nn, but saves a test */
       if (c != 0)
diff -r f3cce34d1228 -r 6393a907d1b5 tests/devel/sqrtrem_1_2.c
--- a/tests/devel/sqrtrem_1_2.c	Tue Mar 07 16:55:42 2017 +0100
+++ b/tests/devel/sqrtrem_1_2.c	Wed Mar 08 19:57:44 2017 +0100
@@ -18,13 +18,13 @@
 
 /* Usage:
 
-   ./sqrtrem_1_2
+   ./sqrtrem_1_2 x
 
      Checks mpn_sqrtrem() exhaustively, starting from 0, incrementing
      the operand by a single unit, until all values handled by
      sqrtrem_{1,2} are tested. SLOW.
 
-   ./sqrtrem_1_2 any_parameter
+   ./sqrtrem_1_2 c
 
      Checks all corner cased for mpn_sqrtrem(). I.e. values of the form
      i*i and (i+1)*(i+1)-1, for each value of i, until all such values,
@@ -52,7 +52,7 @@
 
 int something_wrong (mp_limb_t er, mp_limb_t ec, mp_limb_t es)
 {
-  printf ("root = %lu , rem = {%lu , %lu}\n", (long unsigned) es,(long unsigned) ec,(long unsigned) er);
+  fprintf (stderr, "root = %lu , rem = {%lu , %lu}\n", (long unsigned) es,(long unsigned) ec,(long unsigned) er);
   return -1;
 }
 
@@ -199,7 +199,24 @@
 int
 main (int argc, char **argv)
 {
-  if (argc < 2)
+  int mode = 0;
+
+  if (argc > 1)
+    {
+      if (*argv[1] == 'x')
+	mode = 0;
+      else if (*argv[1] == 'c')
+	mode = 1;
+      else
+	mode = -1;
+      if (argc > 2 || mode == -1)
+	{
+	  fprintf (stderr, "usage: sqrtrem_1_2 [x|c]\n");
+	  exit (1);
+	}
+    }
+
+  if (mode == 0)
     return check_all_values ();
   else
     return check_corner_cases ();


More information about the gmp-commit mailing list