Line | Count | Source |
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 | 672k | ((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 | 192k | { |
43 | 192k | unsigned nbits_in_r; |
44 | 192k | mp_limb_t r; |
45 | 192k | mp_size_t i; |
46 | | |
47 | 192k | if (bi <= nbits) |
48 | 318 | { |
49 | 318 | return p[0] & (((mp_limb_t) 1 << bi) - 1); |
50 | 318 | } |
51 | 192k | else |
52 | 192k | { |
53 | 192k | bi -= nbits; /* bit index of low bit to extract */ |
54 | 192k | i = bi / GMP_NUMB_BITS; /* word index of low bit to extract */ |
55 | 192k | bi %= GMP_NUMB_BITS; /* bit index in low word */ |
56 | 192k | r = p[i] >> bi; /* extract (low) bits */ |
57 | 192k | nbits_in_r = GMP_NUMB_BITS - bi; /* number of bits now in r */ |
58 | 192k | if (nbits_in_r < nbits) /* did we get enough bits? */ |
59 | 17.7k | r += p[i + 1] << nbits_in_r; /* prepend bits from higher word */ |
60 | 192k | return r & (((mp_limb_t) 1 << nbits) - 1); |
61 | 192k | } |
62 | 192k | } |
63 | | |
64 | | static inline unsigned |
65 | | win_size (mp_bitcnt_t eb) |
66 | 622 | { |
67 | 622 | unsigned k; |
68 | 622 | static mp_bitcnt_t x[] = {7,25,81,241,673,1793,4609,11521,28161,~(mp_bitcnt_t)0}; |
69 | 622 | ASSERT (eb > 1); |
70 | 3.79k | for (k = 0; eb > x[k++];) |
71 | 3.17k | ; |
72 | 622 | return k; |
73 | 622 | } |
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 | 622 | { |
85 | 622 | unsigned cnt; |
86 | 622 | mp_bitcnt_t ebi; |
87 | 622 | unsigned windowsize, this_windowsize; |
88 | 622 | mp_limb_t expbits; |
89 | 622 | mp_limb_t *pp; |
90 | 622 | long i; |
91 | 622 | int flipflop; |
92 | 622 | TMP_DECL; |
93 | | |
94 | 622 | ASSERT (en > 1 || (en == 1 && ep[0] > 1)); |
95 | | |
96 | 622 | TMP_MARK; |
97 | | |
98 | 622 | MPN_SIZEINBASE_2EXP(ebi, ep, en, 1); |
99 | | |
100 | 622 | windowsize = win_size (ebi); |
101 | 622 | if (windowsize > 1) |
102 | 607 | { |
103 | 607 | mp_limb_t *this_pp, *last_pp; |
104 | 607 | ASSERT (windowsize < ebi); |
105 | | |
106 | 607 | pp = TMP_ALLOC_LIMBS ((n << (windowsize - 1))); |
107 | | |
108 | 607 | this_pp = pp; |
109 | | |
110 | 607 | MPN_COPY (this_pp, bp, n); |
111 | | |
112 | | /* Store b^2 in tp. */ |
113 | 607 | mpn_sqrlo (tp, bp, n); |
114 | | |
115 | | /* Precompute odd powers of b and put them in the temporary area at pp. */ |
116 | 607 | i = (1 << (windowsize - 1)) - 1; |
117 | 607 | do |
118 | 29.6k | { |
119 | 29.6k | last_pp = this_pp; |
120 | 29.6k | this_pp += n; |
121 | 29.6k | mpn_mullo_n (this_pp, last_pp, tp, n); |
122 | 29.6k | } while (--i != 0); |
123 | | |
124 | 607 | expbits = getbits (ep, ebi, windowsize); |
125 | 607 | ebi -= windowsize; |
126 | | |
127 | | /* THINK: Should we initialise the case expbits % 4 == 0 with a mullo? */ |
128 | 607 | count_trailing_zeros (cnt, expbits); |
129 | 607 | ebi += cnt; |
130 | 607 | expbits >>= cnt; |
131 | | |
132 | 607 | MPN_COPY (rp, pp + n * (expbits >> 1), n); |
133 | 607 | } |
134 | 15 | else |
135 | 15 | { |
136 | 15 | pp = tp + n; |
137 | 15 | MPN_COPY (pp, bp, n); |
138 | 15 | MPN_COPY (rp, bp, n); |
139 | 15 | --ebi; |
140 | 15 | } |
141 | | |
142 | 622 | flipflop = 0; |
143 | | |
144 | 622 | do |
145 | 192k | { |
146 | 672k | while (getbit (ep, ebi) == 0) |
147 | 479k | { |
148 | 479k | mpn_sqrlo (tp, rp, n); |
149 | 479k | MP_PTR_SWAP (rp, tp); |
150 | 479k | flipflop = ! flipflop; |
151 | 479k | if (--ebi == 0) |
152 | 373 | goto done; |
153 | 479k | } |
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 | 192k | expbits = getbits (ep, ebi, windowsize); |
159 | 192k | this_windowsize = MIN (windowsize, ebi); |
160 | | |
161 | 192k | count_trailing_zeros (cnt, expbits); |
162 | 192k | this_windowsize -= cnt; |
163 | 192k | ebi -= this_windowsize; |
164 | 192k | expbits >>= cnt; |
165 | | |
166 | 699k | while (this_windowsize > 1) |
167 | 507k | { |
168 | 507k | mpn_sqrlo (tp, rp, n); |
169 | 507k | mpn_sqrlo (rp, tp, n); |
170 | 507k | this_windowsize -= 2; |
171 | 507k | } |
172 | | |
173 | 192k | if (this_windowsize != 0) |
174 | 126k | mpn_sqrlo (tp, rp, n); |
175 | 65.7k | else |
176 | 65.7k | { |
177 | 65.7k | MP_PTR_SWAP (rp, tp); |
178 | 65.7k | flipflop = ! flipflop; |
179 | 65.7k | } |
180 | | |
181 | 192k | mpn_mullo_n (rp, tp, pp + n * (expbits >> 1), n); |
182 | 192k | } while (ebi != 0); |
183 | | |
184 | 622 | done: |
185 | 622 | if (flipflop) |
186 | 309 | MPN_COPY (tp, rp, n); |
187 | 622 | TMP_FREE; |
188 | 622 | } |