source/libs/gmp/gmp-src/mpn/generic/powlo.c

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