[Gmp-commit] /var/hg/gmp: 2 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Fri Jan 31 16:32:12 UTC 2020
details: /var/hg/gmp/rev/3faec61228bf
changeset: 18025:3faec61228bf
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Fri Jan 31 17:20:44 2020 +0100
description:
Set version so that anyone knows it is experimental.
details: /var/hg/gmp/rev/63e53ddfd210
changeset: 18026:63e53ddfd210
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Fri Jan 31 17:31:42 2020 +0100
description:
mpn/generic/powm.c (MPN_REDC_0): Subtractive redc
(mpn_2powm): New, static, function
diffstat:
gmp-h.in | 2 +-
mpn/generic/powm.c | 403 +++++++++++++++++++++++++++++++++++++++++++++++++++-
2 files changed, 391 insertions(+), 14 deletions(-)
diffs (truncated from 461 to 300 lines):
diff -r 72a89e24a793 -r 63e53ddfd210 gmp-h.in
--- a/gmp-h.in Tue Jan 28 04:56:54 2020 +0100
+++ b/gmp-h.in Fri Jan 31 17:31:42 2020 +0100
@@ -2329,7 +2329,7 @@
/* Major version number is the value of __GNU_MP__ too, above. */
#define __GNU_MP_VERSION 6
#define __GNU_MP_VERSION_MINOR 2
-#define __GNU_MP_VERSION_PATCHLEVEL 0
+#define __GNU_MP_VERSION_PATCHLEVEL 99
#define __GNU_MP_RELEASE (__GNU_MP_VERSION * 10000 + __GNU_MP_VERSION_MINOR * 100 + __GNU_MP_VERSION_PATCHLEVEL)
#define __GMP_H__
diff -r 72a89e24a793 -r 63e53ddfd210 mpn/generic/powm.c
--- a/mpn/generic/powm.c Tue Jan 28 04:56:54 2020 +0100
+++ b/mpn/generic/powm.c Fri Jan 31 17:31:42 2020 +0100
@@ -85,17 +85,17 @@
#include "longlong.h"
#undef MPN_REDC_0
-#define MPN_REDC_0(rp, up, mp, invm) \
+#define MPN_REDC_0(r0, u1, u0, m0, invm) \
do { \
- mp_limb_t p1, r0, u0, _dummy; \
- u0 = *(up); \
- umul_ppmm (p1, _dummy, *(mp), (u0 * (invm)) & GMP_NUMB_MASK); \
- ASSERT (((u0 + _dummy) & GMP_NUMB_MASK) == 0); \
- p1 += (u0 != 0); \
- r0 = (up)[1] + p1; \
- if (p1 > r0) \
- r0 -= *(mp); \
- *(rp) = r0; \
+ mp_limb_t _p1, _u1, _u0, _m0, _r0, _dummy; \
+ _u0 = (u0); \
+ _m0 = (m0); \
+ umul_ppmm (_p1, _dummy, _m0, (_u0 * (invm)) & GMP_NUMB_MASK); \
+ ASSERT (((_u0 - _dummy) & GMP_NUMB_MASK) == 0); \
+ _u1 = (u1); \
+ _r0 = _u1 - _p1; \
+ _r0 = _u1 < _p1 ? _r0 + _m0 : _r0; /* _u1 < _r0 */ \
+ (r0) = _r0 & GMP_NUMB_MASK; \
} while (0)
#undef MPN_REDC_1
@@ -185,6 +185,377 @@
TMP_FREE;
}
+#if ! HAVE_NATIVE_mpn_rsblsh1_n_ip2
+#undef mpn_rsblsh1_n_ip2
+#if HAVE_NATIVE_mpn_rsblsh1_n
+#define mpn_rsblsh1_n_ip2(a,b,n) mpn_rsblsh1_n(a,b,a,n)
+#else
+#define mpn_rsblsh1_n_ip2(a,b,n) \
+ do \
+ { \
+ mpn_lshift (a, a, n, 1); \
+ mpn_sub_n (a, a, b, n); \
+ } while (0)
+#endif
+#endif
+
+#define INNERLOOP2 \
+ do \
+ { \
+ MPN_SQR (tp, rp, n); \
+ MPN_REDUCE (rp, tp, mp, n, mip); \
+ if (mpn_cmp (rp, mp, n) >= 0) \
+ ASSERT_NOCARRY (mpn_sub_n (rp, rp, mp, n)); \
+ if (getbit (ep, ebi) != 0) \
+ { \
+ if (rp[n - 1] >> (mbi - 1) % GMP_LIMB_BITS == 0) \
+ ASSERT_NOCARRY (mpn_lshift (rp, rp, n, 1)); \
+ else \
+ mpn_rsblsh1_n_ip2 (rp, mp, n); \
+ } \
+ } while (--ebi != 0)
+
+/* rp[n-1..0] = 2 ^ ep[en-1..0] mod mp[n-1..0]
+ Requires that mp[n-1..0] is odd and > 1.
+ Requires that ep[en-1..0] is > 1.
+ Uses scratch space at tp of MAX(mpn_binvert_itch(n),2n) limbs. */
+static void
+mpn_2powm (mp_ptr rp, mp_srcptr ep, mp_size_t en,
+ mp_srcptr mp, mp_size_t n, mp_ptr tp)
+{
+ mp_limb_t ip[2], *mip;
+ mp_bitcnt_t ebi, mbi, tbi;
+ mp_size_t tn;
+ int count;
+ TMP_DECL;
+
+ ASSERT (en > 1 || (en == 1 && ep[0] > 1));
+ ASSERT (n > 0 && (mp[0] & 1) != 0);
+
+ MPN_SIZEINBASE_2EXP(ebi, ep, en, 1);
+ MPN_SIZEINBASE_2EXP(mbi, mp, n, 1);
+
+ if (LIKELY (mbi <= GMP_NUMB_MAX))
+ {
+ count_leading_zeros(count, (mp_limb_t) mbi);
+ count = GMP_NUMB_BITS - (count - GMP_NAIL_BITS);
+ }
+ else
+ {
+ mp_bitcnt_t tc = mbi;
+
+ count = 0;
+ do { ++count; } while ((tc >>= 1) != 0);
+ }
+
+ tbi = getbits (ep, ebi, count);
+ if (tbi >= mbi)
+ {
+ --count;
+ ASSERT ((tbi >> count) == 1);
+ tbi >>= 1;
+ ASSERT (tbi < mbi);
+ ASSERT (ebi > count);
+ }
+ else if (ebi <= count)
+ {
+ MPN_FILL (rp, n, 0);
+ rp[tbi / GMP_LIMB_BITS] = CNST_LIMB (1) << (tbi % GMP_LIMB_BITS);
+ return;
+ }
+ ebi -= count;
+
+ if (n == 1)
+ {
+ mp_limb_t r0, m0, invm;
+ m0 = *mp;
+
+ /* redcify (rp, tp, tn + 1, mp, n); */
+ /* TODO: test direct use of udiv_qrnnd */
+ ASSERT (tbi < GMP_LIMB_BITS);
+ tp[1] = CNST_LIMB (1) << tbi;
+ tp[0] = CNST_LIMB (0);
+ r0 = mpn_mod_1 (tp, 2, m0);
+
+ binvert_limb (invm, m0);
+ do
+ {
+ mp_limb_t t0, t1, t2;
+ /* MPN_SQR (tp, rp, n); */
+ umul_ppmm (t1, t0, r0, r0);
+ /* MPN_REDUCE (rp, tp, mp, n, mip); */
+ MPN_REDC_0(r0, t1, t0, m0, invm);
+
+ t2 = r0 << 1;
+ t2 = r0 > (m0 >> 1) ? t2 - m0 : t2;
+ r0 = getbit (ep, ebi) != 0 ? t2 : r0;
+ } while (--ebi != 0);
+
+ /* tp[1] = 0; tp[0] = r0; */
+ /* MPN_REDUCE (rp, tp, mp, n, mip); */
+ MPN_REDC_0(*rp, 0, r0, m0, invm);
+
+ return;
+ }
+
+ TMP_MARK;
+
+#if WANT_REDC_2
+ if (BELOW_THRESHOLD (n, REDC_1_TO_REDC_2_THRESHOLD))
+ {
+ mip = ip;
+ binvert_limb (ip[0], mp[0]);
+ ip[0] = -ip[0];
+ }
+ else if (BELOW_THRESHOLD (n, REDC_2_TO_REDC_N_THRESHOLD))
+ {
+ mip = ip;
+ mpn_binvert (ip, mp, 2, tp);
+ ip[0] = -ip[0]; ip[1] = ~ip[1];
+ }
+#else
+ if (BELOW_THRESHOLD (n, REDC_1_TO_REDC_N_THRESHOLD))
+ {
+ mip = ip;
+ binvert_limb (ip[0], mp[0]);
+ ip[0] = -ip[0];
+ }
+#endif
+ else
+ {
+ mip = TMP_ALLOC_LIMBS (n);
+ mpn_binvert (mip, mp, n, tp);
+ }
+
+ tn = tbi / GMP_LIMB_BITS;
+ MPN_ZERO (tp, tn);
+ tp[tn] = CNST_LIMB (1) << (tbi % GMP_LIMB_BITS);
+
+ redcify (rp, tp, tn + 1, mp, n);
+
+#if WANT_REDC_2
+ if (REDC_1_TO_REDC_2_THRESHOLD < MUL_TOOM22_THRESHOLD)
+ {
+ if (BELOW_THRESHOLD (n, REDC_1_TO_REDC_2_THRESHOLD))
+ {
+ if (REDC_1_TO_REDC_2_THRESHOLD < SQR_BASECASE_THRESHOLD
+ || BELOW_THRESHOLD (n, SQR_BASECASE_THRESHOLD))
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_mul_basecase (r,a,n,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) MPN_REDC_1 (rp, tp, mp, n, mip[0])
+ INNERLOOP2;
+ }
+ else
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_sqr_basecase (r,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) MPN_REDC_1 (rp, tp, mp, n, mip[0])
+ INNERLOOP2;
+ }
+ }
+ else if (BELOW_THRESHOLD (n, MUL_TOOM22_THRESHOLD))
+ {
+ if (MUL_TOOM22_THRESHOLD < SQR_BASECASE_THRESHOLD
+ || BELOW_THRESHOLD (n, SQR_BASECASE_THRESHOLD))
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_mul_basecase (r,a,n,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) MPN_REDC_2 (rp, tp, mp, n, mip)
+ INNERLOOP2;
+ }
+ else
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_sqr_basecase (r,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) MPN_REDC_2 (rp, tp, mp, n, mip)
+ INNERLOOP2;
+ }
+ }
+ else if (BELOW_THRESHOLD (n, REDC_2_TO_REDC_N_THRESHOLD))
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_sqr (r,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) MPN_REDC_2 (rp, tp, mp, n, mip)
+ INNERLOOP2;
+ }
+ else
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_sqr (r,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) mpn_redc_n (rp, tp, mp, n, mip)
+ INNERLOOP2;
+ }
+ }
+ else
+ {
+ if (BELOW_THRESHOLD (n, MUL_TOOM22_THRESHOLD))
+ {
+ if (MUL_TOOM22_THRESHOLD < SQR_BASECASE_THRESHOLD
+ || BELOW_THRESHOLD (n, SQR_BASECASE_THRESHOLD))
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_mul_basecase (r,a,n,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) MPN_REDC_1 (rp, tp, mp, n, mip[0])
+ INNERLOOP2;
+ }
+ else
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_sqr_basecase (r,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) MPN_REDC_1 (rp, tp, mp, n, mip[0])
+ INNERLOOP2;
+ }
+ }
+ else if (BELOW_THRESHOLD (n, REDC_1_TO_REDC_2_THRESHOLD))
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_sqr (r,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) MPN_REDC_1 (rp, tp, mp, n, mip[0])
+ INNERLOOP2;
+ }
+ else if (BELOW_THRESHOLD (n, REDC_2_TO_REDC_N_THRESHOLD))
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_sqr (r,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) MPN_REDC_2 (rp, tp, mp, n, mip)
+ INNERLOOP2;
+ }
+ else
+ {
+#undef MPN_SQR
+#undef MPN_REDUCE
+#define MPN_SQR(r,a,n) mpn_sqr (r,a,n)
+#define MPN_REDUCE(rp,tp,mp,n,mip) mpn_redc_n (rp, tp, mp, n, mip)
+ INNERLOOP2;
More information about the gmp-commit
mailing list