chiark / gitweb /
math/mpx-mul4-test.c: Set `dstr' length correctly in conversion function.
[catacomb] / math / mp-gcd.c
CommitLineData
d3409d5e 1/* -*-c-*-
d3409d5e 2 *
3 * Extended GCD calculation
4 *
5 * (c) 1999 Straylight/Edgeware
6 */
7
45c0fd36 8/*----- Licensing notice --------------------------------------------------*
d3409d5e 9 *
10 * This file is part of Catacomb.
11 *
12 * Catacomb is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU Library General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
45c0fd36 16 *
d3409d5e 17 * Catacomb is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Library General Public License for more details.
45c0fd36 21 *
d3409d5e 22 * You should have received a copy of the GNU Library General Public
23 * License along with Catacomb; if not, write to the Free
24 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
25 * MA 02111-1307, USA.
26 */
27
d3409d5e 28/*----- Header files ------------------------------------------------------*/
29
30#include "mp.h"
31
32/*----- Main code ---------------------------------------------------------*/
33
34/* --- @mp_gcd@ --- *
35 *
36 * Arguments: @mp **gcd, **xx, **yy@ = where to write the results
37 * @mp *a, *b@ = sources (must be nonzero)
38 *
39 * Returns: ---
40 *
41 * Use: Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
42 * @ax + by = gcd(a, b)@. This is useful for computing modular
ef5f4810 43 * inverses.
d3409d5e 44 */
45
46void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
47{
b2253c70 48 mp *x = MP_ONE, *X = MP_ZERO;
49 mp *y = MP_ZERO, *Y = MP_ONE;
d3409d5e 50 mp *u, *v;
3d6115d1 51 mp *q = MP_NEW, *t, *spare = MP_NEW;
ef5f4810 52 unsigned f = 0;
53
ceb3f0c0 54#define f_swap 1u
55#define f_aneg 2u
56#define f_bneg 4u
57#define f_ext 8u
ef5f4810 58
59 /* --- Sort out some initial flags --- */
d3409d5e 60
a7b6288f 61 if (xx || yy) f |= f_ext;
d3409d5e 62
a7b6288f
MW
63 if (MP_NEGP(a)) f |= f_aneg;
64 if (MP_NEGP(b)) f |= f_bneg;
ef5f4810 65
66 /* --- Ensure that @a@ is larger than @b@ --- *
67 *
68 * Use absolute values here!
69 */
70
71 if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
3d6115d1 72 t = a; a = b; b = t;
ef5f4810 73 f |= f_swap;
74 }
75
76 /* --- Check for zeroness --- */
77
3d6115d1 78 if (MP_ZEROP(b)) {
ef5f4810 79
80 /* --- Store %$|a|$% as the GCD --- */
81
82 if (gcd) {
83 if (*gcd) MP_DROP(*gcd);
84 a = MP_COPY(a);
a69a3efd 85 if (MP_NEGP(a)) {
ef5f4810 86 MP_SPLIT(a);
87 a->f &= ~MP_NEG;
88 f |= f_aneg;
89 }
90 *gcd = a;
91 }
92
93 /* --- Store %$1$% and %$0$% in the appropriate bins --- */
94
95 if (f & f_ext) {
96 if (f & f_swap) {
3d6115d1 97 mp **tt = xx; xx = yy; yy = tt;
ef5f4810 98 }
99 if (xx) {
100 if (*xx) MP_DROP(*xx);
a7b6288f
MW
101 if (MP_EQ(a, MP_ZERO)) *xx = MP_ZERO;
102 else if (f & f_aneg) *xx = MP_MONE;
103 else *xx = MP_ONE;
ef5f4810 104 }
105 if (yy) {
106 if (*yy) MP_DROP(*yy);
107 *yy = MP_ZERO;
108 }
109 }
110 return;
d3409d5e 111 }
112
cf76bcbb 113 /* --- Force the signs on the arguments and take copies --- */
d3409d5e 114
115 a = MP_COPY(a);
116 b = MP_COPY(b);
117
ef5f4810 118 MP_SPLIT(a); a->f &= ~MP_NEG;
119 MP_SPLIT(b); b->f &= ~MP_NEG;
120
d3409d5e 121 u = MP_COPY(a);
122 v = MP_COPY(b);
123
cf76bcbb
MW
124 /* --- Main extended Euclidean algorithm --- */
125
a69a3efd 126 while (!MP_ZEROP(v)) {
b2253c70 127 mp_div(&q, &u, u, v);
128 if (f & f_ext) {
3d6115d1 129 t = mp_mul(spare, X, q);
b2253c70 130 t = mp_sub(t, x, t);
3d6115d1
MW
131 spare = x; x = X; X = t;
132 t = mp_mul(spare, Y, q);
b2253c70 133 t = mp_sub(t, y, t);
3d6115d1 134 spare = y; y = Y; Y = t;
d3409d5e 135 }
b2253c70 136 t = u; u = v; v = t;
d3409d5e 137 }
138
3d6115d1 139 MP_DROP(q); if (spare) MP_DROP(spare);
ef5f4810 140 if (!gcd)
b2253c70 141 MP_DROP(u);
ef5f4810 142 else {
143 if (*gcd) MP_DROP(*gcd);
b2253c70 144 u->f &= ~MP_NEG;
145 *gcd = u;
ef5f4810 146 }
d3409d5e 147
148 /* --- Perform a little normalization --- *
149 *
150 * Ensure that the coefficient returned is positive, if there is only one.
ef5f4810 151 * If there are two, favour @y@. Of course, if the original arguments were
152 * negative then I'll need to twiddle their signs as well.
d3409d5e 153 */
154
ef5f4810 155 if (f & f_ext) {
156
157 /* --- If @a@ and @b@ got swapped, swap the coefficients back --- */
158
159 if (f & f_swap) {
3d6115d1 160 t = x; x = y; y = t;
da83a561 161 t = a; a = b; b = t;
d3409d5e 162 }
ef5f4810 163
164 /* --- Sort out the signs --- *
165 *
166 * Note that %$ax + by = a(x - b) + b(y + a)$%.
c9110b35 167 *
168 * This is currently bodgy. It needs sorting out at some time.
ef5f4810 169 */
170
d3409d5e 171 if (yy) {
a69a3efd 172 if (MP_NEGP(y)) {
c9110b35 173 do {
174 y = mp_add(y, y, a);
175 x = mp_sub(x, x, b);
a69a3efd 176 } while (MP_NEGP(y));
c9110b35 177 } else {
178 while (MP_CMP(y, >=, a)) {
179 y = mp_sub(y, y, a);
180 x = mp_add(x, x, b);
181 }
d3409d5e 182 }
c9110b35 183 } else {
a7b6288f
MW
184 if (MP_NEGP(x))
185 do x = mp_add(x, x, b); while (MP_NEGP(x));
186 else
187 while (MP_CMP(x, >=, b)) x = mp_sub(x, x, b);
c9110b35 188 }
d3409d5e 189
ef5f4810 190 /* --- Twiddle the signs --- */
191
9b8ff1cf
MW
192 if (f & f_aneg) { MP_SPLIT(x); x->f ^= MP_NEG; }
193 if (f & f_bneg) { MP_SPLIT(y); y->f ^= MP_NEG; }
ef5f4810 194
195 /* --- Store the results --- */
196
197 if (!xx)
198 MP_DROP(x);
199 else {
200 if (*xx) MP_DROP(*xx);
201 *xx = x;
202 }
203
204 if (!yy)
205 MP_DROP(y);
206 else {
207 if (*yy) MP_DROP(*yy);
208 *yy = y;
209 }
d3409d5e 210 }
211
b2253c70 212 MP_DROP(v);
d3409d5e 213 MP_DROP(X); MP_DROP(Y);
214 MP_DROP(a); MP_DROP(b);
215}
216
b817bfc6 217/* -- @mp_modinv@ --- *
218 *
219 * Arguments: @mp *d@ = destination
220 * @mp *x@ = argument
221 * @mp *p@ = modulus
222 *
223 * Returns: The inverse %$x^{-1} \bmod p$%.
224 *
225 * Use: Computes a modular inverse. An assertion fails if %$p$%
226 * has no inverse.
227 */
228
229mp *mp_modinv(mp *d, mp *x, mp *p)
230{
231 mp *g = MP_NEW;
232 mp_gcd(&g, 0, &d, p, x);
233 assert(MP_EQ(g, MP_ONE));
234 mp_drop(g);
235 return (d);
236}
237
d3409d5e 238/*----- Test rig ----------------------------------------------------------*/
239
240#ifdef TEST_RIG
241
7c404803
MW
242static int modinv(dstr *v)
243{
244 int ok = 1;
245 mp *x = *(mp **)v[0].buf;
246 mp *m = *(mp **)v[1].buf;
247 mp *r = *(mp **)v[2].buf;
248
249 mp *y = mp_modinv(MP_NEW, x, m);
250 if (!MP_EQ(y, r)) {
251 fputs("\n*** mp_modinv failed", stderr);
45c0fd36
MW
252 fputs("\nx = ", stderr); mp_writefile(x, stderr, 10);
253 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
7c404803
MW
254 fputs("\nexpect = ", stderr); mp_writefile(r, stderr, 10);
255 fputs("\nresult = ", stderr); mp_writefile(y, stderr, 10);
256 ok = 0;
257 }
258 MP_DROP(x); MP_DROP(m); MP_DROP(r); MP_DROP(y);
259 assert(mparena_count(MPARENA_GLOBAL) == 0);
45c0fd36 260 return (ok);
7c404803
MW
261}
262
d3409d5e 263static int gcd(dstr *v)
264{
265 int ok = 1;
266 mp *a = *(mp **)v[0].buf;
267 mp *b = *(mp **)v[1].buf;
268 mp *g = *(mp **)v[2].buf;
269 mp *x = *(mp **)v[3].buf;
270 mp *y = *(mp **)v[4].buf;
271
ef5f4810 272 mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW;
d3409d5e 273 mp_gcd(&gg, &xx, &yy, a, b);
b2253c70 274 if (!MP_EQ(x, xx)) {
d3409d5e 275 fputs("\n*** mp_gcd(x) failed", stderr);
45c0fd36
MW
276 fputs("\na = ", stderr); mp_writefile(a, stderr, 10);
277 fputs("\nb = ", stderr); mp_writefile(b, stderr, 10);
d3409d5e 278 fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 10);
279 fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 10);
280 fputc('\n', stderr);
281 ok = 0;
282 }
b2253c70 283 if (!MP_EQ(y, yy)) {
d3409d5e 284 fputs("\n*** mp_gcd(y) failed", stderr);
45c0fd36
MW
285 fputs("\na = ", stderr); mp_writefile(a, stderr, 10);
286 fputs("\nb = ", stderr); mp_writefile(b, stderr, 10);
d3409d5e 287 fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 10);
288 fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 10);
289 fputc('\n', stderr);
290 ok = 0;
291 }
292
293 if (!ok) {
294 mp *ax = mp_mul(MP_NEW, a, xx);
295 mp *by = mp_mul(MP_NEW, b, yy);
296 ax = mp_add(ax, ax, by);
b2253c70 297 if (MP_EQ(ax, gg))
d3409d5e 298 fputs("\n*** (Alternative result found.)\n", stderr);
299 MP_DROP(ax);
300 MP_DROP(by);
301 }
302
b2253c70 303 if (!MP_EQ(g, gg)) {
d3409d5e 304 fputs("\n*** mp_gcd(gcd) failed", stderr);
45c0fd36
MW
305 fputs("\na = ", stderr); mp_writefile(a, stderr, 10);
306 fputs("\nb = ", stderr); mp_writefile(b, stderr, 10);
d3409d5e 307 fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 10);
308 fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 10);
309 fputc('\n', stderr);
310 ok = 0;
311 }
312 MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y);
313 MP_DROP(gg); MP_DROP(xx); MP_DROP(yy);
ef5f4810 314 assert(mparena_count(MPARENA_GLOBAL) == 0);
d3409d5e 315 return (ok);
316}
317
318static test_chunk tests[] = {
319 { "gcd", gcd, { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
7c404803 320 { "modinv", modinv, { &type_mp, &type_mp, &type_mp, 0 } },
d3409d5e 321 { 0, 0, { 0 } }
322};
323
324int main(int argc, char *argv[])
325{
326 sub_init();
0f00dc4c 327 test_run(argc, argv, tests, SRCDIR "/t/mp");
d3409d5e 328 return (0);
329}
330
331#endif
332
333/*----- That's all, folks -------------------------------------------------*/