1 /* mpn_powlo -- Compute R = U^E mod B^n, where B is the limb base.
2
3 Copyright 2007-2009, 2012, 2015, 2016, 2018, 2020 Free Software
4 Foundation, Inc.
5
6 This file is part of the GNU MP Library.
7
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of either:
10
11 * the GNU Lesser General Public License as published by the Free
12 Software Foundation; either version 3 of the License, or (at your
13 option) any later version.
14
15 or
16
17 * the GNU General Public License as published by the Free Software
18 Foundation; either version 2 of the License, or (at your option) any
19 later version.
20
21 or both in parallel, as here.
22
23 The GNU MP Library is distributed in the hope that it will be useful, but
24 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
25 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
26 for more details.
27
28 You should have received copies of the GNU General Public License and the
29 GNU Lesser General Public License along with the GNU MP Library. If not,
30 see https://www.gnu.org/licenses/. */
31
32
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, unsigned nbits)
42 {
43 unsigned 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 unsigned
65 win_size (mp_bitcnt_t eb)
66 {
67 unsigned k;
68 static mp_bitcnt_t x[] = {7,25,81,241,673,1793,4609,11521,28161,~(mp_bitcnt_t)0};
69 ASSERT (eb > 1);
70 for (k = 0; eb > x[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 unsigned cnt;
86 mp_bitcnt_t ebi;
87 unsigned windowsize, this_windowsize;
88 mp_limb_t expbits;
89 mp_limb_t *pp;
90 long i;
91 int flipflop;
92 TMP_DECL;
93
94 ASSERT (en > 1 || (en == 1 && ep[0] > 1));
95
96 TMP_MARK;
97
98 MPN_SIZEINBASE_2EXP(ebi, ep, en, 1);
99
100 windowsize = win_size (ebi);
101 if (windowsize > 1)
102 {
103 mp_limb_t *this_pp, *last_pp;
104 ASSERT (windowsize < ebi);
105
106 pp = TMP_ALLOC_LIMBS ((n << (windowsize - 1)));
107
108 this_pp = pp;
109
110 MPN_COPY (this_pp, bp, n);
111
112 /* Store b^2 in tp. */
113 mpn_sqrlo (tp, bp, n);
114
115 /* Precompute odd powers of b and put them in the temporary area at pp. */
116 i = (1 << (windowsize - 1)) - 1;
117 do
118 {
119 last_pp = this_pp;
120 this_pp += n;
121 mpn_mullo_n (this_pp, last_pp, tp, n);
122 } while (--i != 0);
123
124 expbits = getbits (ep, ebi, windowsize);
125 ebi -= windowsize;
126
127 /* THINK: Should we initialise the case expbits % 4 == 0 with a mullo? */
128 count_trailing_zeros (cnt, expbits);
129 ebi += cnt;
130 expbits >>= cnt;
131
132 MPN_COPY (rp, pp + n * (expbits >> 1), n);
133 }
134 else
135 {
136 pp = tp + n;
137 MPN_COPY (pp, bp, n);
138 MPN_COPY (rp, bp, n);
139 --ebi;
140 }
141
142 flipflop = 0;
143
144 do
145 {
146 while (getbit (ep, ebi) == 0)
147 {
148 mpn_sqrlo (tp, rp, n);
149 MP_PTR_SWAP (rp, tp);
150 flipflop = ! flipflop;
151 if (--ebi == 0)
152 goto done;
153 }
154
155 /* The next bit of the exponent is 1. Now extract the largest block of
156 bits <= windowsize, and such that the least significant bit is 1. */
157
158 expbits = getbits (ep, ebi, windowsize);
159 this_windowsize = MIN (windowsize, ebi);
160
161 count_trailing_zeros (cnt, expbits);
162 this_windowsize -= cnt;
163 ebi -= this_windowsize;
164 expbits >>= cnt;
165
166 while (this_windowsize > 1)
167 {
168 mpn_sqrlo (tp, rp, n);
169 mpn_sqrlo (rp, tp, n);
170 this_windowsize -= 2;
171 }
172
173 if (this_windowsize != 0)
174 mpn_sqrlo (tp, rp, n);
175 else
176 {
177 MP_PTR_SWAP (rp, tp);
178 flipflop = ! flipflop;
179 }
180
181 mpn_mullo_n (rp, tp, pp + n * (expbits >> 1), n);
182 } while (ebi != 0);
183
184 done:
185 if (flipflop)
186 MPN_COPY (tp, rp, n);
187 TMP_FREE;
188 }