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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_powlo -- Compute R = U^E mod B^n, where B is the limb base.
  2. Copyright 2007, 2008, 2009 Free Software Foundation, Inc.
  3. This file is part of the GNU MP Library.
  4. The GNU MP Library is free software; you can redistribute it and/or modify
  5. it under the terms of the GNU Lesser General Public License as published by
  6. the Free Software Foundation; either version 3 of the License, or (at your
  7. option) any later version.
  8. The GNU MP Library is distributed in the hope that it will be useful, but
  9. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  10. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  11. License for more details.
  12. You should have received a copy of the GNU Lesser General Public License
  13. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  14. #include "gmp.h"
  15. #include "gmp-impl.h"
  16. #include "longlong.h"
  17. #define getbit(p,bi) 
  18.   ((p[(bi - 1) / GMP_LIMB_BITS] >> (bi - 1) % GMP_LIMB_BITS) & 1)
  19. static inline mp_limb_t
  20. getbits (const mp_limb_t *p, mp_bitcnt_t bi, int nbits)
  21. {
  22.   int nbits_in_r;
  23.   mp_limb_t r;
  24.   mp_size_t i;
  25.   if (bi < nbits)
  26.     {
  27.       return p[0] & (((mp_limb_t) 1 << bi) - 1);
  28.     }
  29.   else
  30.     {
  31.       bi -= nbits; /* bit index of low bit to extract */
  32.       i = bi / GMP_NUMB_BITS; /* word index of low bit to extract */
  33.       bi %= GMP_NUMB_BITS; /* bit index in low word */
  34.       r = p[i] >> bi; /* extract (low) bits */
  35.       nbits_in_r = GMP_NUMB_BITS - bi; /* number of bits now in r */
  36.       if (nbits_in_r < nbits) /* did we get enough bits? */
  37. r += p[i + 1] << nbits_in_r; /* prepend bits from higher word */
  38.       return r & (((mp_limb_t ) 1 << nbits) - 1);
  39.     }
  40. }
  41. static inline int
  42. win_size (mp_bitcnt_t eb)
  43. {
  44.   int k;
  45.   static mp_bitcnt_t x[] = {1,7,25,81,241,673,1793,4609,11521,28161,~(mp_bitcnt_t)0};
  46.   for (k = 0; eb > x[k]; k++)
  47.     ;
  48.   return k;
  49. }
  50. /* rp[n-1..0] = bp[n-1..0] ^ ep[en-1..0] mod B^n, B is the limb base.
  51.    Requires that ep[en-1] is non-zero.
  52.    Uses scratch space tp[3n-1..0], i.e., 3n words.  */
  53. void
  54. mpn_powlo (mp_ptr rp, mp_srcptr bp,
  55.    mp_srcptr ep, mp_size_t en,
  56.    mp_size_t n, mp_ptr tp)
  57. {
  58.   int cnt;
  59.   mp_bitcnt_t ebi;
  60.   int windowsize, this_windowsize;
  61.   mp_limb_t expbits;
  62.   mp_limb_t *pp, *this_pp, *last_pp;
  63.   mp_limb_t *b2p;
  64.   long i;
  65.   TMP_DECL;
  66.   ASSERT (en > 1 || (en == 1 && ep[0] > 1));
  67.   TMP_MARK;
  68.   count_leading_zeros (cnt, ep[en - 1]);
  69.   ebi = (mp_bitcnt_t) en * GMP_LIMB_BITS - cnt;
  70.   windowsize = win_size (ebi);
  71.   pp = TMP_ALLOC_LIMBS ((n << (windowsize - 1)) + n); /* + n is for mullo ign part */
  72.   this_pp = pp;
  73.   MPN_COPY (this_pp, bp, n);
  74.   b2p = tp + 2*n;
  75.   /* Store b^2 in b2.  */
  76.   mpn_sqr (tp, bp, n); /* FIXME: Use "mpn_sqrlo" */
  77.   MPN_COPY (b2p, tp, n);
  78.   /* Precompute odd powers of b and put them in the temporary area at pp.  */
  79.   for (i = (1 << (windowsize - 1)) - 1; i > 0; i--)
  80.     {
  81.       last_pp = this_pp;
  82.       this_pp += n;
  83.       mpn_mullo_n (this_pp, last_pp, b2p, n);
  84.     }
  85.   expbits = getbits (ep, ebi, windowsize);
  86.   if (ebi < windowsize)
  87.     ebi = 0;
  88.   else
  89.     ebi -= windowsize;
  90.   count_trailing_zeros (cnt, expbits);
  91.   ebi += cnt;
  92.   expbits >>= cnt;
  93.   MPN_COPY (rp, pp + n * (expbits >> 1), n);
  94.   while (ebi != 0)
  95.     {
  96.       while (getbit (ep, ebi) == 0)
  97. {
  98.   mpn_sqr (tp, rp, n); /* FIXME: Use "mpn_sqrlo" */
  99.   MPN_COPY (rp, tp, n);
  100.   ebi--;
  101.   if (ebi == 0)
  102.     goto done;
  103. }
  104.       /* The next bit of the exponent is 1.  Now extract the largest block of
  105.  bits <= windowsize, and such that the least significant bit is 1.  */
  106.       expbits = getbits (ep, ebi, windowsize);
  107.       this_windowsize = windowsize;
  108.       if (ebi < windowsize)
  109. {
  110.   this_windowsize -= windowsize - ebi;
  111.   ebi = 0;
  112. }
  113.       else
  114. ebi -= windowsize;
  115.       count_trailing_zeros (cnt, expbits);
  116.       this_windowsize -= cnt;
  117.       ebi += cnt;
  118.       expbits >>= cnt;
  119.       do
  120. {
  121.   mpn_sqr (tp, rp, n);
  122.   MPN_COPY (rp, tp, n);
  123.   this_windowsize--;
  124. }
  125.       while (this_windowsize != 0);
  126.       mpn_mullo_n (tp, rp, pp + n * (expbits >> 1), n);
  127.       MPN_COPY (rp, tp, n);
  128.     }
  129.  done:
  130.   TMP_FREE;
  131. }