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

数学计算

开发平台:

Unix_Linux

  1. /* mpz_powm(res,base,exp,mod) -- Set R to (U^E) mod M.
  2.    Contributed to the GNU project by Torbjorn Granlund.
  3. Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2005, 2008, 2009
  4. Free Software Foundation, Inc.
  5. This file is part of the GNU MP Library.
  6. The GNU MP Library is free software; you can redistribute it and/or modify
  7. it under the terms of the GNU Lesser General Public License as published by
  8. the Free Software Foundation; either version 3 of the License, or (at your
  9. option) any later version.
  10. The GNU MP Library is distributed in the hope that it will be useful, but
  11. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  13. License for more details.
  14. You should have received a copy of the GNU Lesser General Public License
  15. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  16. #include "gmp.h"
  17. #include "gmp-impl.h"
  18. #include "longlong.h"
  19. #ifdef BERKELEY_MP
  20. #include "mp.h"
  21. #endif
  22. /* TODO
  23.  * Improve handling of buffers.  It is pretty ugly now.
  24.  * For even moduli, we compute a binvert of its odd part both here and in
  25.    mpn_powm.  How can we avoid this recomputation?
  26. */
  27. /*
  28.   b ^ e mod m   res
  29.   0   0     0    ?
  30.   0   e     0    ?
  31.   0   0     m    ?
  32.   0   e     m    0
  33.   b   0     0    ?
  34.   b   e     0    ?
  35.   b   0     m    1 mod m
  36.   b   e     m    b^e mod m
  37. */
  38. #define HANDLE_NEGATIVE_EXPONENT 1
  39. void
  40. #ifndef BERKELEY_MP
  41. mpz_powm (mpz_ptr r, mpz_srcptr b, mpz_srcptr e, mpz_srcptr m)
  42. #else /* BERKELEY_MP */
  43. pow (mpz_srcptr b, mpz_srcptr e, mpz_srcptr m, mpz_ptr r)
  44. #endif /* BERKELEY_MP */
  45. {
  46.   mp_size_t n, nodd, ncnt;
  47.   int cnt;
  48.   mp_ptr rp, tp;
  49.   mp_srcptr bp, ep, mp;
  50.   mp_size_t rn, bn, es, en, itch;
  51.   TMP_DECL;
  52.   n = ABSIZ(m);
  53.   if (n == 0)
  54.     DIVIDE_BY_ZERO;
  55.   mp = PTR(m);
  56.   TMP_MARK;
  57.   es = SIZ(e);
  58.   if (UNLIKELY (es <= 0))
  59.     {
  60.       mpz_t new_b;
  61.       if (es == 0)
  62. {
  63.   /* b^0 mod m,  b is anything and m is non-zero.
  64.      Result is 1 mod m, i.e., 1 or 0 depending on if m = 1.  */
  65.   SIZ(r) = n != 1 || mp[0] != 1;
  66.   PTR(r)[0] = 1;
  67.   TMP_FREE; /* we haven't really allocated anything here */
  68.   return;
  69. }
  70. #if HANDLE_NEGATIVE_EXPONENT
  71.       MPZ_TMP_INIT (new_b, n + 1);
  72.       if (! mpz_invert (new_b, b, m))
  73. DIVIDE_BY_ZERO;
  74.       b = new_b;
  75.       es = -es;
  76. #else
  77.       DIVIDE_BY_ZERO;
  78. #endif
  79.     }
  80.   en = es;
  81.   bn = ABSIZ(b);
  82.   if (UNLIKELY (bn == 0))
  83.     {
  84.       SIZ(r) = 0;
  85.       TMP_FREE;
  86.       return;
  87.     }
  88.   ep = PTR(e);
  89.   /* Handle (b^1 mod m) early, since mpn_pow* do not handle that case.  */
  90.   if (UNLIKELY (en == 1 && ep[0] == 1))
  91.     {
  92.       rp = TMP_ALLOC_LIMBS (n);
  93.       bp = PTR(b);
  94.       if (bn >= n)
  95. {
  96.   mp_ptr qp = TMP_ALLOC_LIMBS (bn - n + 1);
  97.   mpn_tdiv_qr (qp, rp, 0L, bp, bn, mp, n);
  98.   rn = n;
  99.   MPN_NORMALIZE (rp, rn);
  100.   if (SIZ(b) < 0 && rn != 0)
  101.     {
  102.       mpn_sub (rp, mp, n, rp, rn);
  103.       rn = n;
  104.       MPN_NORMALIZE (rp, rn);
  105.     }
  106. }
  107.       else
  108. {
  109.   if (SIZ(b) < 0)
  110.     {
  111.       mpn_sub (rp, mp, n, bp, bn);
  112.       rn = n;
  113.       rn -= (rp[rn - 1] == 0);
  114.     }
  115.   else
  116.     {
  117.       MPN_COPY (rp, bp, bn);
  118.       rn = bn;
  119.     }
  120. }
  121.       goto ret;
  122.     }
  123.   /* Remove low zero limbs from M.  This loop will terminate for correctly
  124.      represented mpz numbers.  */
  125.   ncnt = 0;
  126.   while (UNLIKELY (mp[0] == 0))
  127.     {
  128.       mp++;
  129.       ncnt++;
  130.     }
  131.   nodd = n - ncnt;
  132.   cnt = 0;
  133.   if (mp[0] % 2 == 0)
  134.     {
  135.       mp_ptr new = TMP_ALLOC_LIMBS (nodd);
  136.       count_trailing_zeros (cnt, mp[0]);
  137.       mpn_rshift (new, mp, nodd, cnt);
  138.       nodd -= new[nodd - 1] == 0;
  139.       mp = new;
  140.       ncnt++;
  141.     }
  142.   if (ncnt != 0)
  143.     {
  144.       /* We will call both mpn_powm and mpn_powlo.  */
  145.       /* rp needs n, mpn_powlo needs 4n, the 2 mpn_binvert might need more */
  146.       mp_size_t n_largest_binvert = MAX (ncnt, nodd);
  147.       mp_size_t itch_binvert = mpn_binvert_itch (n_largest_binvert);
  148.       itch = 3 * n + MAX (itch_binvert, 2 * n);
  149.     }
  150.   else
  151.     {
  152.       /* We will call just mpn_powm.  */
  153.       mp_size_t itch_binvert = mpn_binvert_itch (nodd);
  154.       itch = n + MAX (itch_binvert, 2 * n);
  155.     }
  156.   tp = TMP_ALLOC_LIMBS (itch);
  157.   rp = tp;  tp += n;
  158.   bp = PTR(b);
  159.   mpn_powm (rp, bp, bn, ep, en, mp, nodd, tp);
  160.   rn = n;
  161.   if (ncnt != 0)
  162.     {
  163.       mp_ptr r2, xp, yp, odd_inv_2exp;
  164.       unsigned long t;
  165.       int bcnt;
  166.       if (bn < ncnt)
  167. {
  168.   mp_ptr new = TMP_ALLOC_LIMBS (ncnt);
  169.   MPN_COPY (new, bp, bn);
  170.   MPN_ZERO (new + bn, ncnt - bn);
  171.   bp = new;
  172. }
  173.       r2 = tp;
  174.       if (bp[0] % 2 == 0)
  175. {
  176.   if (en > 1)
  177.     {
  178.       MPN_ZERO (r2, ncnt);
  179.       goto zero;
  180.     }
  181.   ASSERT (en == 1);
  182.   t = (ncnt - (cnt != 0)) * GMP_NUMB_BITS + cnt;
  183.   /* Count number of low zero bits in B, up to 3.  */
  184.   bcnt = (0x1213 >> ((bp[0] & 7) << 1)) & 0x3;
  185.   /* Note that ep[0] * bcnt might overflow, but that just results
  186.      in a missed optimization.  */
  187.   if (ep[0] * bcnt >= t)
  188.     {
  189.       MPN_ZERO (r2, ncnt);
  190.       goto zero;
  191.     }
  192. }
  193.       mpn_powlo (r2, bp, ep, en, ncnt, tp + ncnt);
  194.     zero:
  195.       if (nodd < ncnt)
  196. {
  197.   mp_ptr new = TMP_ALLOC_LIMBS (ncnt);
  198.   MPN_COPY (new, mp, nodd);
  199.   MPN_ZERO (new + nodd, ncnt - nodd);
  200.   mp = new;
  201. }
  202.       odd_inv_2exp = tp + n;
  203.       mpn_binvert (odd_inv_2exp, mp, ncnt, tp + 2 * n);
  204.       mpn_sub (r2, r2, ncnt, rp, nodd > ncnt ? ncnt : nodd);
  205.       xp = tp + 2 * n;
  206.       mpn_mullo_n (xp, odd_inv_2exp, r2, ncnt);
  207.       if (cnt != 0)
  208. xp[ncnt - 1] &= (CNST_LIMB(1) << cnt) - 1;
  209.       yp = tp;
  210.       if (ncnt > nodd)
  211. mpn_mul (yp, xp, ncnt, mp, nodd);
  212.       else
  213. mpn_mul (yp, mp, nodd, xp, ncnt);
  214.       mpn_add (rp, yp, n, rp, nodd);
  215.       ASSERT (nodd + ncnt >= n);
  216.       ASSERT (nodd + ncnt <= n + 1);
  217.     }
  218.   MPN_NORMALIZE (rp, rn);
  219.   if ((ep[0] & 1) && SIZ(b) < 0 && rn != 0)
  220.     {
  221.       mpn_sub (rp, PTR(m), n, rp, rn);
  222.       rn = n;
  223.       MPN_NORMALIZE (rp, rn);
  224.     }
  225.  ret:
  226.   MPZ_REALLOC (r, rn);
  227.   SIZ(r) = rn;
  228.   MPN_COPY (PTR(r), rp, rn);
  229.   TMP_FREE;
  230. }