[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