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