1 /* mpn_powlo -- Compute R = U^E mod B^n, where B is the limb base.
2
3 Copyright 2007-2009, 2012 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
getbits(const mp_limb_t * p,mp_bitcnt_t bi,int nbits)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
win_size(mp_bitcnt_t eb)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 for (k = 0; eb > x[k]; k++)
70 ;
71 return k;
72 }
73
74 /* rp[n-1..0] = bp[n-1..0] ^ ep[en-1..0] mod B^n, B is the limb base.
75 Requires that ep[en-1] is non-zero.
76 Uses scratch space tp[3n-1..0], i.e., 3n words. */
77 void
mpn_powlo(mp_ptr rp,mp_srcptr bp,mp_srcptr ep,mp_size_t en,mp_size_t n,mp_ptr tp)78 mpn_powlo (mp_ptr rp, mp_srcptr bp,
79 mp_srcptr ep, mp_size_t en,
80 mp_size_t n, mp_ptr tp)
81 {
82 int cnt;
83 mp_bitcnt_t ebi;
84 int windowsize, this_windowsize;
85 mp_limb_t expbits;
86 mp_limb_t *pp, *this_pp, *last_pp;
87 mp_limb_t *b2p;
88 long i;
89 TMP_DECL;
90
91 ASSERT (en > 1 || (en == 1 && ep[0] > 1));
92
93 TMP_MARK;
94
95 MPN_SIZEINBASE_2EXP(ebi, ep, en, 1);
96
97 windowsize = win_size (ebi);
98
99 pp = TMP_ALLOC_LIMBS ((n << (windowsize - 1)) + n); /* + n is for mullo ign part */
100
101 this_pp = pp;
102
103 MPN_COPY (this_pp, bp, n);
104
105 b2p = tp + 2*n;
106
107 /* Store b^2 in b2. */
108 mpn_sqr (tp, bp, n); /* FIXME: Use "mpn_sqrlo" */
109 MPN_COPY (b2p, tp, 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, b2p, n);
117 }
118
119 expbits = getbits (ep, ebi, windowsize);
120 if (ebi < windowsize)
121 ebi = 0;
122 else
123 ebi -= windowsize;
124
125 count_trailing_zeros (cnt, expbits);
126 ebi += cnt;
127 expbits >>= cnt;
128
129 MPN_COPY (rp, pp + n * (expbits >> 1), n);
130
131 while (ebi != 0)
132 {
133 while (getbit (ep, ebi) == 0)
134 {
135 mpn_sqr (tp, rp, n); /* FIXME: Use "mpn_sqrlo" */
136 MPN_COPY (rp, tp, n);
137 ebi--;
138 if (ebi == 0)
139 goto done;
140 }
141
142 /* The next bit of the exponent is 1. Now extract the largest block of
143 bits <= windowsize, and such that the least significant bit is 1. */
144
145 expbits = getbits (ep, ebi, windowsize);
146 this_windowsize = windowsize;
147 if (ebi < windowsize)
148 {
149 this_windowsize -= windowsize - ebi;
150 ebi = 0;
151 }
152 else
153 ebi -= windowsize;
154
155 count_trailing_zeros (cnt, expbits);
156 this_windowsize -= cnt;
157 ebi += cnt;
158 expbits >>= cnt;
159
160 do
161 {
162 mpn_sqr (tp, rp, n);
163 MPN_COPY (rp, tp, n);
164 this_windowsize--;
165 }
166 while (this_windowsize != 0);
167
168 mpn_mullo_n (tp, rp, pp + n * (expbits >> 1), n);
169 MPN_COPY (rp, tp, n);
170 }
171
172 done:
173 TMP_FREE;
174 }
175