Use sliding-window exponentiation.
[u/mdw/catacomb] / mpcrt.c
CommitLineData
5ee4c893 1/* -*-c-*-
2 *
22bab86c 3 * $Id: mpcrt.c,v 1.3 2000/10/08 12:11:22 mdw Exp $
5ee4c893 4 *
5 * Chinese Remainder Theorem computations (Gauss's algorithm)
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: mpcrt.c,v $
22bab86c 33 * Revision 1.3 2000/10/08 12:11:22 mdw
34 * Use @MP_EQ@ instead of @MP_CMP@.
35 *
5bda60bd 36 * Revision 1.2 1999/12/10 23:22:32 mdw
37 * Interface changes for suggested destinations. Use Barrett reduction.
38 *
5ee4c893 39 * Revision 1.1 1999/11/22 20:50:57 mdw
40 * Add support for solving Chinese Remainder Theorem problems.
41 *
42 */
43
44/*----- Header files ------------------------------------------------------*/
45
46#include "mp.h"
47#include "mpcrt.h"
5bda60bd 48#include "mpbarrett.h"
5ee4c893 49
50/*----- Main code ---------------------------------------------------------*/
51
52/* --- @mpcrt_create@ --- *
53 *
54 * Arguments: @mpcrt *c@ = pointer to CRT context
55 * @mpcrt_mod *v@ = pointer to vector of moduli
56 * @size_t k@ = number of moduli
57 * @mp *n@ = product of all moduli (@MP_NEW@ if unknown)
58 *
59 * Returns: ---
60 *
61 * Use: Initializes a context for solving Chinese Remainder Theorem
62 * problems. The vector of moduli can be incomplete. Omitted
63 * items must be left as null pointers. Not all combinations of
64 * missing things can be coped with, even if there is
65 * technically enough information to cope. For example, if @n@
66 * is unspecified, all the @m@ values must be present, even if
67 * there is one modulus with both @m@ and @n@ (from which the
68 * product of all moduli could clearly be calculated).
69 */
70
71void mpcrt_create(mpcrt *c, mpcrt_mod *v, size_t k, mp *n)
72{
5ee4c893 73 size_t i;
74
75 /* --- Simple initialization things --- */
76
77 c->k = k;
78 c->v = v;
79
80 /* --- Work out @n@ if I don't have it already --- */
81
5bda60bd 82 if (n != MP_NEW)
83 n = MP_COPY(n);
84 else {
5ee4c893 85 n = MP_COPY(v[0].m);
5bda60bd 86 for (i = 1; i < k; i++)
87 n = mp_mul(n, n, v[i].m);
88 }
89
90 /* --- A quick hack if %$k = 2$% --- */
91
92 if (k == 2) {
93
94 /* --- The %$n / n_i$% values are trivial in this case --- */
95
96 if (!v[0].n)
97 v[0].n = MP_COPY(v[1].m);
98 if (!v[1].n)
99 v[1].n = MP_COPY(v[0].m);
100
101 /* --- Now sort out the inverses --- *
102 *
103 * @mp_gcd@ will ensure that the first argument is negative.
104 */
105
106 if (!v[0].ni && !v[1].ni) {
107 mp_gcd(0, &v[0].ni, &v[1].ni, v[0].n, v[1].n);
108 v[0].ni = mp_add(v[0].ni, v[0].ni, v[1].n);
109 } else {
110 int i, j;
111 mp *x;
112
113 if (!v[0].ni)
114 i = 0, j = 1;
115 else
116 i = 1, j = 0;
117
118 x = mp_mul(MP_NEW, v[j].n, v[j].ni);
119 x = mp_sub(x, x, MP_ONE);
120 mp_div(&x, 0, x, v[i].n);
121 v[i].ni = x;
5ee4c893 122 }
123 }
124
5bda60bd 125 /* --- Set up the Barrett context --- */
5ee4c893 126
5bda60bd 127 mpbarrett_create(&c->mb, n);
5ee4c893 128
129 /* --- Walk through filling in @n@, @ni@ and @nnir@ --- */
130
131 for (i = 0; i < k; i++) {
132 if (!v[i].n)
133 mp_div(&v[i].n, 0, n, v[i].m);
134 if (!v[i].ni)
135 mp_gcd(0, &v[i].ni, 0, v[i].n, v[i].m);
5bda60bd 136 if (!v[i].nni)
137 v[i].nni = mp_mul(MP_NEW, v[i].n, v[i].ni);
5ee4c893 138 }
139
140 /* --- Done --- */
141
5bda60bd 142 mp_drop(n);
5ee4c893 143}
144
145/* --- @mpcrt_destroy@ --- *
146 *
147 * Arguments: @mpcrt *c@ - pointer to CRT context
148 *
149 * Returns: ---
150 *
151 * Use: Destroys a CRT context, releasing all the resources it holds.
152 */
153
154void mpcrt_destroy(mpcrt *c)
155{
156 size_t i;
157
158 for (i = 0; i < c->k; i++) {
159 if (c->v[i].m) mp_drop(c->v[i].m);
160 if (c->v[i].n) mp_drop(c->v[i].n);
161 if (c->v[i].ni) mp_drop(c->v[i].ni);
5bda60bd 162 if (c->v[i].nni) mp_drop(c->v[i].nni);
5ee4c893 163 }
5bda60bd 164 mpbarrett_destroy(&c->mb);
5ee4c893 165}
166
167/* --- @mpcrt_solve@ --- *
168 *
169 * Arguments: @mpcrt *c@ = pointer to CRT context
5bda60bd 170 * @mp *d@ = fake destination
5ee4c893 171 * @mp **v@ = array of residues
172 *
173 * Returns: The unique solution modulo the product of the individual
174 * moduli, which leaves the given residues.
175 *
176 * Use: Constructs a result given its residue modulo an array of
177 * coprime integers. This can be used to improve performance of
178 * RSA encryption or Blum-Blum-Shub generation if the factors
179 * of the modulus are known, since results can be computed mod
180 * each of the individual factors and then combined at the end.
181 * This is rather faster than doing the full-scale modular
182 * exponentiation.
183 */
184
5bda60bd 185mp *mpcrt_solve(mpcrt *c, mp *d, mp **v)
5ee4c893 186{
187 mp *a = MP_ZERO;
188 mp *x = MP_NEW;
189 size_t i;
190
191 for (i = 0; i < c->k; i++) {
5bda60bd 192 x = mp_mul(x, c->v[i].nni, v[i]);
193 x = mpbarrett_reduce(&c->mb, x, x);
5ee4c893 194 a = mp_add(a, a, x);
195 }
196 if (x)
5bda60bd 197 MP_DROP(x);
198 a = mpbarrett_reduce(&c->mb, a, a);
199 if (d != MP_NEW)
200 MP_DROP(d);
5ee4c893 201 return (a);
202}
203
204/*----- Test rig ----------------------------------------------------------*/
205
206#ifdef TEST_RIG
207
208static int verify(size_t n, dstr *v)
209{
210 mpcrt_mod *m = xmalloc(n * sizeof(mpcrt_mod));
211 mp **r = xmalloc(n * sizeof(mp *));
212 mpcrt c;
213 mp *a, *b;
214 size_t i;
215 int ok = 1;
216
217 for (i = 0; i < n; i++) {
218 r[i] = *(mp **)v[2 * i].buf;
219 m[i].m = *(mp **)v[2 * i + 1].buf;
220 m[i].n = 0;
221 m[i].ni = 0;
5bda60bd 222 m[i].nni = 0;
5ee4c893 223 }
224 a = *(mp **)v[2 * n].buf;
225
226 mpcrt_create(&c, m, n, 0);
5bda60bd 227 b = mpcrt_solve(&c, MP_NEW, r);
5ee4c893 228
22bab86c 229 if (!MP_EQ(a, b)) {
5ee4c893 230 fputs("\n*** failed\n", stderr);
231 fputs("n = ", stderr);
5bda60bd 232 mp_writefile(c.mb.m, stderr, 10);
5ee4c893 233 for (i = 0; i < n; i++) {
234 fprintf(stderr, "\nr[%u] = ", i);
235 mp_writefile(r[i], stderr, 10);
236 fprintf(stderr, "\nm[%u] = ", i);
237 mp_writefile(m[i].m, stderr, 10);
238 fprintf(stderr, "\nN[%u] = ", i);
239 mp_writefile(m[i].n, stderr, 10);
240 fprintf(stderr, "\nM[%u] = ", i);
241 mp_writefile(m[i].ni, stderr, 10);
242 }
243 fputs("\nresult = ", stderr);
244 mp_writefile(b, stderr, 10);
245 fputs("\nexpect = ", stderr);
246 mp_writefile(a, stderr, 10);
247 fputc('\n', stderr);
248 ok = 0;
249 }
250
5bda60bd 251 for (i = 0; i < n; i++)
252 mp_drop(r[i]);
5ee4c893 253 mp_drop(a);
254 mp_drop(b);
255 mpcrt_destroy(&c);
256 free(m);
257 free(r);
5bda60bd 258 assert(mparena_count(MPARENA_GLOBAL) == 0);
5ee4c893 259 return (ok);
260}
261
262static int crt1(dstr *v) { return verify(1, v); }
263static int crt2(dstr *v) { return verify(2, v); }
264static int crt3(dstr *v) { return verify(3, v); }
265static int crt4(dstr *v) { return verify(4, v); }
266static int crt5(dstr *v) { return verify(5, v); }
267
268static test_chunk tests[] = {
269 { "crt-1", crt1, { &type_mp, &type_mp,
270 &type_mp, 0 } },
271 { "crt-2", crt2, { &type_mp, &type_mp,
272 &type_mp, &type_mp,
273 &type_mp, 0 } },
274 { "crt-3", crt3, { &type_mp, &type_mp,
275 &type_mp, &type_mp,
276 &type_mp, &type_mp,
277 &type_mp, 0 } },
278 { "crt-4", crt4, { &type_mp, &type_mp,
279 &type_mp, &type_mp,
280 &type_mp, &type_mp,
281 &type_mp, &type_mp,
282 &type_mp, 0 } },
283 { "crt-5", crt5, { &type_mp, &type_mp,
284 &type_mp, &type_mp,
285 &type_mp, &type_mp,
286 &type_mp, &type_mp,
287 &type_mp, &type_mp,
288 &type_mp, 0 } },
289 { 0, 0, { 0 } }
290};
291
292int main(int argc, char *argv[])
293{
294 sub_init();
295 test_run(argc, argv, tests, SRCDIR "/tests/mpcrt");
296 return (0);
297}
298
299#endif
300
301/*----- That's all, folks -------------------------------------------------*/