powm_ui.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:4k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* mpz_powm_ui(res,base,exp,mod) -- Set RES to (base**exp) mod MOD.
  2. Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2005 Free Software
  3. Foundation, Inc.
  4. This file is part of the GNU MP Library.
  5. The GNU MP Library is free software; you can redistribute it and/or modify
  6. it under the terms of the GNU Lesser General Public License as published by
  7. the Free Software Foundation; either version 3 of the License, or (at your
  8. option) any later version.
  9. The GNU MP Library is distributed in the hope that it will be useful, but
  10. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  11. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  12. License for more details.
  13. You should have received a copy of the GNU Lesser General Public License
  14. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  15. #include "gmp.h"
  16. #include "gmp-impl.h"
  17. #include "longlong.h"
  18. /* Compute t = a mod m, a is defined by (ap,an), m is defined by (mp,mn), and
  19.    t is defined by (tp,mn).  */
  20. static void
  21. reduce (mp_ptr tp, mp_srcptr ap, mp_size_t an, mp_srcptr mp, mp_size_t mn)
  22. {
  23.   mp_ptr qp;
  24.   TMP_DECL;
  25.   TMP_MARK;
  26.   qp = TMP_ALLOC_LIMBS (an - mn + 1);
  27.   mpn_tdiv_qr (qp, tp, 0L, ap, an, mp, mn);
  28.   TMP_FREE;
  29. }
  30. void
  31. mpz_powm_ui (mpz_ptr r, mpz_srcptr b, unsigned long int el, mpz_srcptr m)
  32. {
  33.   mp_ptr xp, tp, qp, mp, bp;
  34.   mp_size_t xn, tn, mn, bn;
  35.   int m_zero_cnt;
  36.   int c;
  37.   mp_limb_t e;
  38.   TMP_DECL;
  39.   mp = PTR(m);
  40.   mn = ABSIZ(m);
  41.   if (mn == 0)
  42.     DIVIDE_BY_ZERO;
  43.   if (el == 0)
  44.     {
  45.       /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0
  46.  depending on if MOD equals 1.  */
  47.       SIZ(r) = (mn == 1 && mp[0] == 1) ? 0 : 1;
  48.       PTR(r)[0] = 1;
  49.       return;
  50.     }
  51.   TMP_MARK;
  52.   /* Normalize m (i.e. make its most significant bit set) as required by
  53.      division functions below.  */
  54.   count_leading_zeros (m_zero_cnt, mp[mn - 1]);
  55.   m_zero_cnt -= GMP_NAIL_BITS;
  56.   if (m_zero_cnt != 0)
  57.     {
  58.       mp_ptr new_mp = TMP_ALLOC_LIMBS (mn);
  59.       mpn_lshift (new_mp, mp, mn, m_zero_cnt);
  60.       mp = new_mp;
  61.     }
  62.   bn = ABSIZ(b);
  63.   bp = PTR(b);
  64.   if (bn > mn)
  65.     {
  66.       /* Reduce possibly huge base.  Use a function call to reduce, since we
  67.  don't want the quotient allocation to live until function return.  */
  68.       mp_ptr new_bp = TMP_ALLOC_LIMBS (mn);
  69.       reduce (new_bp, bp, bn, mp, mn);
  70.       bp = new_bp;
  71.       bn = mn;
  72.       /* Canonicalize the base, since we are potentially going to multiply with
  73.  it quite a few times.  */
  74.       MPN_NORMALIZE (bp, bn);
  75.     }
  76.   if (bn == 0)
  77.     {
  78.       SIZ(r) = 0;
  79.       TMP_FREE;
  80.       return;
  81.     }
  82.   tp = TMP_ALLOC_LIMBS (2 * mn + 1);
  83.   xp = TMP_ALLOC_LIMBS (mn);
  84.   qp = TMP_ALLOC_LIMBS (mn + 1);
  85.   MPN_COPY (xp, bp, bn);
  86.   xn = bn;
  87.   e = el;
  88.   count_leading_zeros (c, e);
  89.   e = (e << c) << 1; /* shift the exp bits to the left, lose msb */
  90.   c = GMP_LIMB_BITS - 1 - c;
  91.   /* Main loop. */
  92.   /* If m is already normalized (high bit of high limb set), and b is the
  93.      same size, but a bigger value, and e==1, then there's no modular
  94.      reductions done and we can end up with a result out of range at the
  95.      end. */
  96.   if (c == 0)
  97.     {
  98.       if (xn == mn && mpn_cmp (xp, mp, mn) >= 0)
  99.         mpn_sub_n (xp, xp, mp, mn);
  100.       goto finishup;
  101.     }
  102.   while (c != 0)
  103.     {
  104.       mpn_sqr (tp, xp, xn);
  105.       tn = 2 * xn; tn -= tp[tn - 1] == 0;
  106.       if (tn < mn)
  107. {
  108.   MPN_COPY (xp, tp, tn);
  109.   xn = tn;
  110. }
  111.       else
  112. {
  113.   mpn_tdiv_qr (qp, xp, 0L, tp, tn, mp, mn);
  114.   xn = mn;
  115. }
  116.       if ((mp_limb_signed_t) e < 0)
  117. {
  118.   mpn_mul (tp, xp, xn, bp, bn);
  119.   tn = xn + bn; tn -= tp[tn - 1] == 0;
  120.   if (tn < mn)
  121.     {
  122.       MPN_COPY (xp, tp, tn);
  123.       xn = tn;
  124.     }
  125.   else
  126.     {
  127.       mpn_tdiv_qr (qp, xp, 0L, tp, tn, mp, mn);
  128.       xn = mn;
  129.     }
  130. }
  131.       e <<= 1;
  132.       c--;
  133.     }
  134.  finishup:
  135.   /* We shifted m left m_zero_cnt steps.  Adjust the result by reducing
  136.      it with the original MOD.  */
  137.   if (m_zero_cnt != 0)
  138.     {
  139.       mp_limb_t cy;
  140.       cy = mpn_lshift (tp, xp, xn, m_zero_cnt);
  141.       tp[xn] = cy; xn += cy != 0;
  142.       if (xn < mn)
  143. {
  144.   MPN_COPY (xp, tp, xn);
  145. }
  146.       else
  147. {
  148.   mpn_tdiv_qr (qp, xp, 0L, tp, xn, mp, mn);
  149.   xn = mn;
  150. }
  151.       mpn_rshift (xp, xp, xn, m_zero_cnt);
  152.     }
  153.   MPN_NORMALIZE (xp, xn);
  154.   if ((el & 1) != 0 && SIZ(b) < 0 && xn != 0)
  155.     {
  156.       mp = PTR(m); /* want original, unnormalized m */
  157.       mpn_sub (xp, mp, mn, xp, xn);
  158.       xn = mn;
  159.       MPN_NORMALIZE (xp, xn);
  160.     }
  161.   MPZ_REALLOC (r, xn);
  162.   SIZ (r) = xn;
  163.   MPN_COPY (PTR(r), xp, xn);
  164.   TMP_FREE;
  165. }