289 lines
7.0 KiB
C
289 lines
7.0 KiB
C
/* mpz_powm_ui(res,base,exp,mod) -- Set R to (B^E) mod M.
|
|
|
|
Contributed to the GNU project by Torbjörn Granlund.
|
|
|
|
Copyright 1991, 1993, 1994, 1996, 1997, 2000-2002, 2005, 2008, 2009, 2011-2013
|
|
Free Software Foundation, Inc.
|
|
|
|
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 either:
|
|
|
|
* the GNU Lesser General Public License as published by the Free
|
|
Software Foundation; either version 3 of the License, or (at your
|
|
option) any later version.
|
|
|
|
or
|
|
|
|
* the GNU General Public License as published by the Free Software
|
|
Foundation; either version 2 of the License, or (at your option) any
|
|
later version.
|
|
|
|
or both in parallel, as here.
|
|
|
|
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 General Public License
|
|
for more details.
|
|
|
|
You should have received copies of the GNU General Public License and the
|
|
GNU Lesser General Public License along with the GNU MP Library. If not,
|
|
see https://www.gnu.org/licenses/. */
|
|
|
|
|
|
#include "mpir.h"
|
|
#include "gmp-impl.h"
|
|
#include "longlong.h"
|
|
|
|
|
|
/* This code is very old, and should be rewritten to current GMP standard. It
|
|
is slower than mpz_powm for large exponents, but also for small exponents
|
|
when the mod argument is small.
|
|
|
|
As an intermediate solution, we now deflect to mpz_powm for exponents >= 20.
|
|
*/
|
|
|
|
/*
|
|
b ^ e mod m res
|
|
0 0 0 ?
|
|
0 e 0 ?
|
|
0 0 m ?
|
|
0 e m 0
|
|
b 0 0 ?
|
|
b e 0 ?
|
|
b 0 m 1 mod m
|
|
b e m b^e mod m
|
|
*/
|
|
|
|
static void
|
|
mod (mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, mp_limb_t dinv,
|
|
mp_ptr tp)
|
|
{
|
|
mp_ptr qp;
|
|
TMP_DECL;
|
|
TMP_MARK;
|
|
|
|
qp = tp;
|
|
|
|
if (dn == 1)
|
|
{
|
|
np[0] = mpn_divrem_1 (qp, (mp_size_t) 0, np, nn, dp[0]); /* JPF */
|
|
}
|
|
else if (dn == 2)
|
|
{
|
|
mpn_divrem_2 (qp, 0L, np, nn, dp);
|
|
}
|
|
else if (BELOW_THRESHOLD (dn, DC_DIV_QR_THRESHOLD) ||
|
|
BELOW_THRESHOLD (nn - dn, DC_DIV_QR_THRESHOLD))
|
|
{
|
|
mpn_sb_div_qr (qp, np, nn, dp, dn, dinv); /* JPF: no gmp_pi1_t, two limbs pi */
|
|
}
|
|
/* Different conditions */
|
|
else if (BELOW_THRESHOLD (dn, INV_DIV_QR_THRESHOLD) || /* fast condition */
|
|
BELOW_THRESHOLD (nn, 2 * INV_DIV_QR_THRESHOLD)) /* fast condition */
|
|
{
|
|
mpn_dc_div_qr (qp, np, nn, dp, dn, dinv); /* JPF: no gmp_pi1_t */
|
|
}
|
|
else
|
|
{
|
|
mp_ptr dinv2 = TMP_ALLOC_LIMBS(dn); /* JPF: ... */
|
|
mpn_invert(dinv2, dp, dn); /* JPF: ... */
|
|
mpn_inv_div_qr (qp, np, nn, dp, dn, dinv2); /* JPF: need nn+1 for new np? I don't think so as everyting is already normalized */
|
|
}
|
|
|
|
TMP_FREE;
|
|
}
|
|
|
|
/* Compute t = a mod m, a is defined by (ap,an), m is defined by (mp,mn), and
|
|
t is defined by (tp,mn). */
|
|
static void
|
|
reduce (mp_ptr tp, mp_srcptr ap, mp_size_t an, mp_srcptr mp, mp_size_t mn,
|
|
mp_limb_t dinv)
|
|
{
|
|
mp_ptr rp, scratch;
|
|
TMP_DECL;
|
|
TMP_MARK;
|
|
|
|
rp = TMP_ALLOC_LIMBS (an);
|
|
scratch = TMP_ALLOC_LIMBS (an - mn + 1);
|
|
MPN_COPY (rp, ap, an);
|
|
mod (rp, an, mp, mn, dinv, scratch);
|
|
MPN_COPY (tp, rp, mn);
|
|
|
|
TMP_FREE;
|
|
}
|
|
|
|
void
|
|
mpz_powm_ui (mpz_ptr r, mpz_srcptr b, mpir_ui el, mpz_srcptr m)
|
|
{
|
|
if (el < 20) /* JPF */
|
|
{
|
|
mp_ptr xp, tp, mp, bp, scratch;
|
|
mp_size_t xn, tn, mn, bn;
|
|
int m_zero_cnt;
|
|
int c;
|
|
mp_limb_t e, m2;
|
|
mp_limb_t dinv;
|
|
TMP_DECL;
|
|
|
|
mp = PTR(m);
|
|
mn = ABSIZ(m);
|
|
if (UNLIKELY (mn == 0))
|
|
DIVIDE_BY_ZERO;
|
|
|
|
if (el == 0)
|
|
{
|
|
/* Exponent is zero, result is 1 mod M, i.e., 1 or 0 depending on if
|
|
M equals 1. */
|
|
SIZ(r) = (mn == 1 && mp[0] == 1) ? 0 : 1;
|
|
PTR(r)[0] = 1;
|
|
return;
|
|
}
|
|
|
|
TMP_MARK;
|
|
|
|
/* Normalize m (i.e. make its most significant bit set) as required by
|
|
division functions below. */
|
|
count_leading_zeros (m_zero_cnt, mp[mn - 1]);
|
|
m_zero_cnt -= GMP_NAIL_BITS;
|
|
if (m_zero_cnt != 0)
|
|
{
|
|
mp_ptr new_mp = TMP_ALLOC_LIMBS (mn);
|
|
mpn_lshift (new_mp, mp, mn, m_zero_cnt);
|
|
mp = new_mp;
|
|
}
|
|
|
|
m2 = mn == 1 ? 0 : mp[mn - 2];
|
|
mpir_invert_pi1 (dinv, mp[mn - 1], m2); /* JPF: don't use gmp_pi1_t */
|
|
|
|
bn = ABSIZ(b);
|
|
bp = PTR(b);
|
|
if (bn > mn)
|
|
{
|
|
/* Reduce possibly huge base. Use a function call to reduce, since we
|
|
don't want the quotient allocation to live until function return. */
|
|
mp_ptr new_bp = TMP_ALLOC_LIMBS (mn);
|
|
reduce (new_bp, bp, bn, mp, mn, dinv); /* JPF */
|
|
bp = new_bp;
|
|
bn = mn;
|
|
/* Canonicalize the base, since we are potentially going to multiply with
|
|
it quite a few times. */
|
|
MPN_NORMALIZE (bp, bn);
|
|
}
|
|
|
|
if (bn == 0)
|
|
{
|
|
SIZ(r) = 0;
|
|
TMP_FREE;
|
|
return;
|
|
}
|
|
|
|
tp = TMP_ALLOC_LIMBS (2 * mn + 1);
|
|
xp = TMP_ALLOC_LIMBS (mn);
|
|
scratch = TMP_ALLOC_LIMBS (mn + 1);
|
|
|
|
MPN_COPY (xp, bp, bn);
|
|
xn = bn;
|
|
|
|
e = el;
|
|
count_leading_zeros (c, e);
|
|
e = (e << c) << 1; /* shift the exp bits to the left, lose msb */
|
|
c = GMP_LIMB_BITS - 1 - c;
|
|
|
|
if (c == 0)
|
|
{
|
|
/* If m is already normalized (high bit of high limb set), and b is
|
|
the same size, but a bigger value, and e==1, then there's no
|
|
modular reductions done and we can end up with a result out of
|
|
range at the end. */
|
|
if (xn == mn && mpn_cmp (xp, mp, mn) >= 0)
|
|
mpn_sub_n (xp, xp, mp, mn);
|
|
}
|
|
else
|
|
{
|
|
/* Main loop. */
|
|
do
|
|
{
|
|
mpn_sqr (tp, xp, xn);
|
|
tn = 2 * xn; tn -= tp[tn - 1] == 0;
|
|
if (tn < mn)
|
|
{
|
|
MPN_COPY (xp, tp, tn);
|
|
xn = tn;
|
|
}
|
|
else
|
|
{
|
|
mod (tp, tn, mp, mn, dinv, scratch); /* JPF */
|
|
MPN_COPY (xp, tp, mn);
|
|
xn = mn;
|
|
}
|
|
|
|
if ((mp_limb_signed_t) e < 0)
|
|
{
|
|
mpn_mul (tp, xp, xn, bp, bn);
|
|
tn = xn + bn; tn -= tp[tn - 1] == 0;
|
|
if (tn < mn)
|
|
{
|
|
MPN_COPY (xp, tp, tn);
|
|
xn = tn;
|
|
}
|
|
else
|
|
{
|
|
mod (tp, tn, mp, mn, dinv, scratch); /* JPF */
|
|
MPN_COPY (xp, tp, mn);
|
|
xn = mn;
|
|
}
|
|
}
|
|
e <<= 1;
|
|
c--;
|
|
}
|
|
while (c != 0);
|
|
}
|
|
|
|
/* We shifted m left m_zero_cnt steps. Adjust the result by reducing it
|
|
with the original M. */
|
|
if (m_zero_cnt != 0)
|
|
{
|
|
mp_limb_t cy;
|
|
cy = mpn_lshift (tp, xp, xn, m_zero_cnt);
|
|
tp[xn] = cy; xn += cy != 0;
|
|
|
|
if (xn < mn)
|
|
{
|
|
MPN_COPY (xp, tp, xn);
|
|
}
|
|
else
|
|
{
|
|
mod (tp, xn, mp, mn, dinv, scratch); /* JPF */
|
|
MPN_COPY (xp, tp, mn);
|
|
xn = mn;
|
|
}
|
|
mpn_rshift (xp, xp, xn, m_zero_cnt);
|
|
}
|
|
MPN_NORMALIZE (xp, xn);
|
|
|
|
if ((el & 1) != 0 && SIZ(b) < 0 && xn != 0)
|
|
{
|
|
mp = PTR(m); /* want original, unnormalized m */
|
|
mpn_sub (xp, mp, mn, xp, xn);
|
|
xn = mn;
|
|
MPN_NORMALIZE (xp, xn);
|
|
}
|
|
MPZ_REALLOC (r, xn);
|
|
SIZ (r) = xn;
|
|
MPN_COPY (PTR(r), xp, xn);
|
|
|
|
TMP_FREE;
|
|
}
|
|
else /* e >= 20 */
|
|
{
|
|
/* For large exponents, fake a mpz_t exponent and deflect to the more
|
|
sophisticated mpz_powm. */
|
|
mpz_t e;
|
|
mp_limb_t ep[LIMBS_PER_UI]; /* JPF: no ulong in MPIR */
|
|
MPZ_FAKE_UI (e, ep, el);
|
|
mpz_powm (r, b, e, m);
|
|
}
|
|
}
|