Moved the Karatsuba macros into a separate file for better sharing.
[u/mdw/catacomb] / mpmont-mexp.c
CommitLineData
4cc02d06 1/* -*-c-*-
2 *
ef5f4810 3 * $Id: mpmont-mexp.c,v 1.3 1999/12/10 23:18:39 mdw Exp $
4cc02d06 4 *
5 * Multiplle simultaneous exponentiations
6 *
7 * (c) 1999 Straylight/Edgeware
8 */
9
10/*----- Licensing notice --------------------------------------------------*
11 *
12 * This file is part of Catacomb.
13 *
14 * Catacomb is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU Library General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
18 *
19 * Catacomb is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU Library General Public License for more details.
23 *
24 * You should have received a copy of the GNU Library General Public
25 * License along with Catacomb; if not, write to the Free
26 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
27 * MA 02111-1307, USA.
28 */
29
30/*----- Revision history --------------------------------------------------*
31 *
32 * $Log: mpmont-mexp.c,v $
ef5f4810 33 * Revision 1.3 1999/12/10 23:18:39 mdw
34 * Change interface for suggested destinations.
35 *
79a34029 36 * Revision 1.2 1999/11/21 11:35:10 mdw
37 * Performance improvement: use @mp_sqr@ and @mpmont_reduce@ instead of
38 * @mpmont_mul@ for squaring in exponentiation.
39 *
4cc02d06 40 * Revision 1.1 1999/11/19 13:19:29 mdw
41 * Simultaneous exponentiation support.
42 *
43 */
44
45/*----- Header files ------------------------------------------------------*/
46
47#include "mp.h"
48#include "mpmont.h"
49
50/*----- Main code ---------------------------------------------------------*/
51
52/* --- @mpmont_mexpr@ --- *
53 *
54 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
ef5f4810 55 * @mp *d@ = fake destination
4cc02d06 56 * @mpmont_factor *f@ = pointer to array of factors
57 * @size_t n@ = number of factors supplied
58 *
59 * Returns: If the bases are %$g_0, g_1, \ldots, g_{n-1}$% and the
60 * exponents are %$e_0, e_1, \ldots, e_{n-1}$% then the result
61 * is:
62 *
63 * %$g_0^{e_0} g_1^{e_1} \ldots g_{n-1}^{e_{n-1}} R \bmod m$%
64 */
65
66typedef struct scan {
67 size_t len;
68 mpw w;
69} scan;
70
ef5f4810 71mp *mpmont_mexpr(mpmont *mm, mp *d, mpmont_factor *f, size_t n)
4cc02d06 72{
73 size_t vn = 1 << n;
74 mp **v = xmalloc(vn * sizeof(mp *));
75 scan *s;
76 size_t o;
77 unsigned b;
78 mp *a = MP_COPY(mm->r);
79 mp *spare = MP_NEW;
80
81 /* --- Perform the precomputation --- */
82
83 {
84 size_t i, j;
85 size_t mask;
86
87 /* --- Fill in the rest of the array --- *
88 *
89 * Zero never gets used.
90 */
91
92 j = 0;
93 mask = 0;
94 for (i = 1; i < vn; i++) {
95
96 /* --- Check for a new bit entering --- *
97 *
98 * If a bit gets set that wasn't set before, then all the lower bits
99 * are zeroes and I've got to introduce a new base into the array.
100 */
101
102 if ((i & mask) == 0) {
103 v[i] = mpmont_mul(mm, MP_NEW, f[j++].base, mm->r2);
104 mask = i;
105 }
106
107 /* --- Otherwise I can get away with a single multiplication --- *
108 *
109 * In particular, if %$i$% has more than one bit set, then I only need
110 * to calculate %$v_i = v_{\mathit{mask}} v_{i - \mathit{mask}}$%.
111 * Since both are less than %$i$%, they must have already been
112 * computed.
113 */
114
115 else
116 v[i] = mpmont_mul(mm, MP_NEW, v[mask], v[i & ~mask]);
117 }
118 }
119
120 /* --- Set up the bitscanners --- *
121 *
122 * I must scan the exponents from left to right, which is a shame. It
123 * means that I can't use the standard @mpscan@ stuff, in particular.
124 */
125
126 {
127 size_t i;
128
129 s = xmalloc(n * sizeof(scan));
130 o = 0;
131 for (i = 0; i < n; i++) {
132 s[i].len = MP_LEN(f[i].exp);
133 if (s[i].len > o)
134 o = s[i].len;
135 }
136 b = 0;
137 }
138
139 /* --- Now do the actual calculation --- */
140
141 b = 0;
142 for (;;) {
143 size_t i;
144 size_t j;
145 mp *dd;
146
147 /* --- If no more bits, get some more --- */
148
149 if (!b) {
150 if (!o)
151 break;
152 o--;
153 b = MPW_BITS;
154 }
155
156 /* --- Work out the next index --- */
157
158 j = 0;
159 b--;
160 for (i = 0; i < n; i++) {
161 if (o < s[i].len)
162 j |= (((f[i].exp->v[o] >> b) & 1) << i);
163 }
164
165 /* --- Accumulate the result --- */
166
167 if (spare) {
79a34029 168 dd = mp_sqr(spare, a);
169 dd = mpmont_reduce(mm, dd, dd);
4cc02d06 170 spare = a;
171 a = dd;
172 }
173
174 if (j) {
175 dd = mpmont_mul(mm, spare, a, v[j]);
176 spare = a;
177 a = dd;
178 }
179 }
180
181 /* --- Tidy up afterwards --- */
182
183 {
184 size_t i;
185 for (i = 1; i < vn; i++)
186 MP_DROP(v[i]);
187 if (spare)
188 MP_DROP(spare);
189 free(v);
190 free(s);
191 }
192
ef5f4810 193 if (d != MP_NEW)
194 MP_DROP(d);
195
4cc02d06 196 return (a);
197}
198
199/* --- @mpmont_mexp@ --- *
200 *
201 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
ef5f4810 202 * @mp *d@ = fake destination
4cc02d06 203 * @mpmont_factor *f@ = pointer to array of factors
204 * @size_t n@ = number of factors supplied
205 *
206 * Returns: Product of bases raised to exponents, all mod @m@.
207 *
208 * Use: Convenient interface over @mpmont_mexpr@.
209 */
210
ef5f4810 211mp *mpmont_mexp(mpmont *mm, mp *d, mpmont_factor *f, size_t n)
4cc02d06 212{
ef5f4810 213 d = mpmont_mexpr(mm, d, f, n);
4cc02d06 214 d = mpmont_reduce(mm, d, d);
215 return (d);
216}
217
218/*----- Test rig ----------------------------------------------------------*/
219
220#ifdef TEST_RIG
221
222#include <mLib/testrig.h>
223
224static int verify(size_t n, dstr *v)
225{
226 mp *m = *(mp **)v[0].buf;
227 mpmont_factor *f = xmalloc(n * sizeof(*f));
228 mp *r, *rr;
229 size_t i, j;
230 mpmont mm;
231 int ok = 1;
232
233 j = 1;
234 for (i = 0; i < n; i++) {
235 f[i].base = *(mp **)v[j++].buf;
236 f[i].exp = *(mp **)v[j++].buf;
237 }
238
239 rr = *(mp **)v[j].buf;
240 mpmont_create(&mm, m);
ef5f4810 241 r = mpmont_mexp(&mm, MP_NEW, f, n);
4cc02d06 242 if (MP_CMP(r, !=, rr)) {
243 fputs("\n*** mexp failed\n", stderr);
244 fputs("m = ", stderr); mp_writefile(m, stderr, 10);
245 for (i = 0; i < n; i++) {
246 fprintf(stderr, "\ng_%u = ", i);
247 mp_writefile(f[i].base, stderr, 10);
248 fprintf(stderr, "\ne_%u = ", i);
249 mp_writefile(f[i].exp, stderr, 10);
250 }
251 fputs("\nr = ", stderr); mp_writefile(r, stderr, 10);
252 fputs("\nR = ", stderr); mp_writefile(rr, stderr, 10);
253 fputc('\n', stderr);
254 ok = 0;
255 }
256
257 for (i = 0; i < n; i++) {
258 MP_DROP(f[i].base);
259 MP_DROP(f[i].exp);
260 }
261 MP_DROP(m);
262 MP_DROP(r);
263 MP_DROP(rr);
264 mpmont_destroy(&mm);
ef5f4810 265 assert(mparena_count(MPARENA_GLOBAL) == 0);
4cc02d06 266 return (ok);
267}
268
269static int t1(dstr *v) { return verify(1, v); }
270static int t2(dstr *v) { return verify(2, v); }
271static int t3(dstr *v) { return verify(3, v); }
272static int t4(dstr *v) { return verify(4, v); }
273static int t5(dstr *v) { return verify(5, v); }
274
275static test_chunk tests[] = {
276 { "mexp-1", t1, { &type_mp,
277 &type_mp, &type_mp,
278 &type_mp, 0 } },
279 { "mexp-2", t2, { &type_mp,
280 &type_mp, &type_mp,
281 &type_mp, &type_mp,
282 &type_mp, 0 } },
283 { "mexp-3", t3, { &type_mp,
284 &type_mp, &type_mp,
285 &type_mp, &type_mp,
286 &type_mp, &type_mp,
287 &type_mp, 0 } },
288 { "mexp-4", t4, { &type_mp,
289 &type_mp, &type_mp,
290 &type_mp, &type_mp,
291 &type_mp, &type_mp,
292 &type_mp, &type_mp,
293 &type_mp, 0 } },
294 { "mexp-5", t5, { &type_mp,
295 &type_mp, &type_mp,
296 &type_mp, &type_mp,
297 &type_mp, &type_mp,
298 &type_mp, &type_mp,
299 &type_mp, &type_mp,
300 &type_mp, 0 } },
301 { 0, 0, { 0 } }
302};
303
304int main(int argc, char *argv[])
305{
306 sub_init();
307 test_run(argc, argv, tests, SRCDIR "/tests/mpmont");
308 return (0);
309}
310
311#endif
312
313/*----- That's all, folks -------------------------------------------------*/