--- /dev/null
+/* mpz_powm(res,base,exp,mod) -- Set RES to (base**exp) mod MOD.
+
+Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 Free Software Foundation, Inc.
+Contributed by Paul Zimmermann.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+#ifdef BERKELEY_MP
+#include "mp.h"
+#endif
+
+
+/* set c <- (a*b)/R^n mod m c has to have at least (2n) allocated limbs */
+static void
+#if __STDC__
+mpz_redc (mpz_ptr c, mpz_srcptr a, mpz_srcptr b, mpz_srcptr m, mp_limb_t Nprim)
+#else
+mpz_redc (c, a, b, m, Nprim)
+ mpz_ptr c;
+ mpz_srcptr a;
+ mpz_srcptr b;
+ mpz_srcptr m;
+ mp_limb_t Nprim;
+#endif
+{
+ mp_ptr cp, mp = PTR (m);
+ mp_limb_t cy, cout = 0;
+ mp_limb_t q;
+ size_t j, n = ABSIZ (m);
+
+ ASSERT (ALLOC (c) >= 2 * n);
+
+ mpz_mul (c, a, b);
+ cp = PTR (c);
+ j = ABSIZ (c);
+ MPN_ZERO (cp + j, 2 * n - j);
+ for (j = 0; j < n; j++)
+ {
+ q = cp[0] * Nprim;
+ cy = mpn_addmul_1 (cp, mp, n, q);
+ cout += mpn_add_1 (cp + n, cp + n, n - j, cy);
+ cp++;
+ }
+ cp -= n;
+ if (cout)
+ {
+ cy = cout - mpn_sub_n (cp, cp + n, mp, n);
+ while (cy)
+ cy -= mpn_sub_n (cp, cp, mp, n);
+ }
+ else
+ MPN_COPY (cp, cp + n, n);
+ MPN_NORMALIZE (cp, n);
+ SIZ (c) = SIZ (c) < 0 ? -n : n;
+}
+
+/* average number of calls to redc for an exponent of n bits
+ with the sliding window algorithm of base 2^k: the optimal is
+ obtained for the value of k which minimizes 2^(k-1)+n/(k+1):
+
+ n\k 4 5 6 7 8
+ 128 156* 159 171 200 261
+ 256 309 307* 316 343 403
+ 512 617 607* 610 632 688
+ 1024 1231 1204 1195* 1207 1256
+ 2048 2461 2399 2366 2360* 2396
+ 4096 4918 4787 4707 4665* 4670
+*/
+\f
+#ifndef BERKELEY_MP
+void
+#if __STDC__
+mpz_powm (mpz_ptr res, mpz_srcptr base, mpz_srcptr e, mpz_srcptr mod)
+#else
+mpz_powm (res, base, e, mod)
+ mpz_ptr res;
+ mpz_srcptr base;
+ mpz_srcptr e;
+ mpz_srcptr mod;
+#endif
+#else /* BERKELEY_MP */
+void
+#if __STDC__
+pow (mpz_srcptr base, mpz_srcptr e, mpz_srcptr mod, mpz_ptr res)
+#else
+pow (base, e, mod, res)
+ mpz_srcptr base;
+ mpz_srcptr e;
+ mpz_srcptr mod;
+ mpz_ptr res;
+#endif
+#endif /* BERKELEY_MP */
+{
+ mp_limb_t invm, *ep, c, mask;
+ mpz_t xx, *g;
+ mp_size_t n, i, K, j, l, k;
+ int sh;
+ int use_redc;
+
+#ifdef POWM_DEBUG
+ mpz_t exp;
+ mpz_init (exp);
+#endif
+
+ n = ABSIZ (mod);
+
+ if (n == 0)
+ DIVIDE_BY_ZERO;
+
+ if (SIZ (e) == 0)
+ {
+ /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0
+ depending on if MOD equals 1. */
+ SIZ(res) = (ABSIZ (mod) == 1 && (PTR(mod))[0] == 1) ? 0 : 1;
+ PTR(res)[0] = 1;
+ return;
+ }
+
+ /* Use REDC instead of usual reduction for sizes < POWM_THRESHOLD.
+ In REDC each modular multiplication costs about 2*n^2 limbs operations,
+ whereas using usual reduction it costs 3*K(n), where K(n) is the cost of a
+ multiplication using Karatsuba, and a division is assumed to cost 2*K(n),
+ for example using Burnikel-Ziegler's algorithm. This gives a theoretical
+ threshold of a*KARATSUBA_SQR_THRESHOLD, with a=(3/2)^(1/(2-ln(3)/ln(2))) ~
+ 2.66. */
+ /* For now, also disable REDC when MOD is even, as the inverse can't
+ handle that. */
+
+#ifndef POWM_THRESHOLD
+#define POWM_THRESHOLD ((8 * KARATSUBA_SQR_THRESHOLD) / 3)
+#endif
+
+ use_redc = (n < POWM_THRESHOLD && PTR(mod)[0] % 2 != 0);
+ if (use_redc)
+ {
+ /* invm = -1/m mod 2^BITS_PER_MP_LIMB, must have m odd */
+ modlimb_invert (invm, PTR(mod)[0]);
+ invm = -invm;
+ }
+
+ /* determines optimal value of k */
+ l = ABSIZ (e) * BITS_PER_MP_LIMB; /* number of bits of exponent */
+ k = 1;
+ K = 2;
+ while (2 * l > K * (2 + k * (3 + k)))
+ {
+ k++;
+ K *= 2;
+ }
+
+ g = (mpz_t *) (*_mp_allocate_func) (K / 2 * sizeof (mpz_t));
+ /* compute x*R^n where R=2^BITS_PER_MP_LIMB */
+ mpz_init (g[0]);
+ if (use_redc)
+ {
+ mpz_mul_2exp (g[0], base, n * BITS_PER_MP_LIMB);
+ mpz_mod (g[0], g[0], mod);
+ }
+ else
+ mpz_mod (g[0], base, mod);
+
+ /* compute xx^g for odd g < 2^k */
+ mpz_init (xx);
+ if (use_redc)
+ {
+ _mpz_realloc (xx, 2 * n);
+ mpz_redc (xx, g[0], g[0], mod, invm); /* xx = x^2*R^n */
+ }
+ else
+ {
+ mpz_mul (xx, g[0], g[0]);
+ mpz_mod (xx, xx, mod);
+ }
+ for (i = 1; i < K / 2; i++)
+ {
+ mpz_init (g[i]);
+ if (use_redc)
+ {
+ _mpz_realloc (g[i], 2 * n);
+ mpz_redc (g[i], g[i - 1], xx, mod, invm); /* g[i] = x^(2i+1)*R^n */
+ }
+ else
+ {
+ mpz_mul (g[i], g[i - 1], xx);
+ mpz_mod (g[i], g[i], mod);
+ }
+ }
+
+ /* now starts the real stuff */
+ mask = (mp_limb_t) ((1<<k) - 1);
+ ep = PTR (e);
+ i = ABSIZ (e) - 1; /* current index */
+ c = ep[i]; /* current limb */
+ count_leading_zeros (sh, c);
+ sh = BITS_PER_MP_LIMB - sh; /* significant bits in ep[i] */
+ sh -= k; /* index of lower bit of ep[i] to take into account */
+ if (sh < 0)
+ { /* k-sh extra bits are needed */
+ if (i > 0)
+ {
+ i--;
+ c = (c << (-sh)) | (ep[i] >> (BITS_PER_MP_LIMB + sh));
+ sh += BITS_PER_MP_LIMB;
+ }
+ }
+ else
+ c = c >> sh;
+#ifdef POWM_DEBUG
+ printf ("-1/m mod 2^%u = %lu\n", BITS_PER_MP_LIMB, invm);
+ mpz_set_ui (exp, c);
+#endif
+ j=0;
+ while (c % 2 == 0)
+ {
+ j++;
+ c = (c >> 1);
+ }
+ mpz_set (xx, g[c >> 1]);
+ while (j--)
+ {
+ if (use_redc)
+ mpz_redc (xx, xx, xx, mod, invm);
+ else
+ {
+ mpz_mul (xx, xx, xx);
+ mpz_mod (xx, xx, mod);
+ }
+ }
+
+#ifdef POWM_DEBUG
+ printf ("x^"); mpz_out_str (0, 10, exp);
+ printf ("*2^%u mod m = ", n * BITS_PER_MP_LIMB); mpz_out_str (0, 10, xx);
+ putchar ('\n');
+#endif
+
+ while (i > 0 || sh > 0)
+ {
+ c = ep[i];
+ sh -= k;
+ l = k; /* number of bits treated */
+ if (sh < 0)
+ {
+ if (i > 0)
+ {
+ i--;
+ c = (c << (-sh)) | (ep[i] >> (BITS_PER_MP_LIMB + sh));
+ sh += BITS_PER_MP_LIMB;
+ }
+ else
+ {
+ l += sh; /* may be less bits than k here */
+ c = c & ((1<<l) - 1);
+ }
+ }
+ else
+ c = c >> sh;
+ c = c & mask;
+
+ /* this while loop implements the sliding window improvement */
+ while ((c & (1 << (k - 1))) == 0 && (i > 0 || sh > 0))
+ {
+ if (use_redc) mpz_redc (xx, xx, xx, mod, invm);
+ else
+ {
+ mpz_mul (xx, xx, xx);
+ mpz_mod (xx, xx, mod);
+ }
+ if (sh)
+ {
+ sh--;
+ c = (c<<1) + ((ep[i]>>sh) & 1);
+ }
+ else
+ {
+ i--;
+ sh = BITS_PER_MP_LIMB - 1;
+ c = (c<<1) + (ep[i]>>sh);
+ }
+ }
+
+#ifdef POWM_DEBUG
+ printf ("l=%u c=%lu\n", l, c);
+ mpz_mul_2exp (exp, exp, k);
+ mpz_add_ui (exp, exp, c);
+#endif
+
+ /* now replace xx by xx^(2^k)*x^c */
+ if (c != 0)
+ {
+ j = 0;
+ while (c % 2 == 0)
+ {
+ j++;
+ c = c >> 1;
+ }
+ /* c0 = c * 2^j, i.e. xx^(2^k)*x^c = (A^(2^(k - j))*c)^(2^j) */
+ l -= j;
+ while (l--)
+ if (use_redc) mpz_redc (xx, xx, xx, mod, invm);
+ else
+ {
+ mpz_mul (xx, xx, xx);
+ mpz_mod (xx, xx, mod);
+ }
+ if (use_redc)
+ mpz_redc (xx, xx, g[c >> 1], mod, invm);
+ else
+ {
+ mpz_mul (xx, xx, g[c >> 1]);
+ mpz_mod (xx, xx, mod);
+ }
+ }
+ else
+ j = l; /* case c=0 */
+ while (j--)
+ {
+ if (use_redc)
+ mpz_redc (xx, xx, xx, mod, invm);
+ else
+ {
+ mpz_mul (xx, xx, xx);
+ mpz_mod (xx, xx, mod);
+ }
+ }
+#ifdef POWM_DEBUG
+ printf ("x^"); mpz_out_str (0, 10, exp);
+ printf ("*2^%u mod m = ", n * BITS_PER_MP_LIMB); mpz_out_str (0, 10, xx);
+ putchar ('\n');
+#endif
+ }
+
+ /* now convert back xx to xx/R^n */
+ if (use_redc)
+ {
+ mpz_set_ui (g[0], 1);
+ mpz_redc (xx, xx, g[0], mod, invm);
+ if (mpz_cmp (xx, mod) >= 0)
+ mpz_sub (xx, xx, mod);
+ }
+ mpz_set (res, xx);
+
+ mpz_clear (xx);
+ for (i = 0; i < K / 2; i++)
+ mpz_clear (g[i]);
+ (*_mp_free_func) (g, K / 2 * sizeof (mpz_t));
+}
projects / ghc-hetmet.git / blobdiff