3 * Extended GCD calculation
5 * (c) 1999 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
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.
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.
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,
28 /*----- Header files ------------------------------------------------------*/
32 /*----- Main code ---------------------------------------------------------*/
36 * Arguments: @mp **gcd, **xx, **yy@ = where to write the results
37 * @mp *a, *b@ = sources (must be nonzero)
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
46 void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
48 mp *x = MP_ONE, *X = MP_ZERO;
49 mp *y = MP_ZERO, *Y = MP_ONE;
51 mp *q = MP_NEW, *t, *spare = MP_NEW;
59 /* --- Sort out some initial flags --- */
61 if (xx || yy) f |= f_ext;
63 if (MP_NEGP(a)) f |= f_aneg;
64 if (MP_NEGP(b)) f |= f_bneg;
66 /* --- Ensure that @a@ is larger than @b@ --- *
68 * Use absolute values here!
71 if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
76 /* --- Check for zeroness --- */
80 /* --- Store %$|a|$% as the GCD --- */
83 if (*gcd) MP_DROP(*gcd);
93 /* --- Store %$1$% and %$0$% in the appropriate bins --- */
97 mp **tt = xx; xx = yy; yy = tt;
100 if (*xx) MP_DROP(*xx);
101 if (MP_EQ(a, MP_ZERO)) *xx = MP_ZERO;
102 else if (f & f_aneg) *xx = MP_MONE;
106 if (*yy) MP_DROP(*yy);
113 /* --- Force the signs on the arguments and take copies --- */
118 MP_SPLIT(a); a->f &= ~MP_NEG;
119 MP_SPLIT(b); b->f &= ~MP_NEG;
124 /* --- Main extended Euclidean algorithm --- */
126 while (!MP_ZEROP(v)) {
127 mp_div(&q, &u, u, v);
129 t = mp_mul(spare, X, q);
131 spare = x; x = X; X = t;
132 t = mp_mul(spare, Y, q);
134 spare = y; y = Y; Y = t;
139 MP_DROP(q); if (spare) MP_DROP(spare);
143 if (*gcd) MP_DROP(*gcd);
148 /* --- Perform a little normalization --- *
150 * Ensure that the coefficient returned is positive, if there is only one.
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.
157 /* --- If @a@ and @b@ got swapped, swap the coefficients back --- */
164 /* --- Sort out the signs --- *
166 * Note that %$ax + by = a(x - b) + b(y + a)$%.
168 * This is currently bodgy. It needs sorting out at some time.
176 } while (MP_NEGP(y));
178 while (MP_CMP(y, >=, a)) {
185 do x = mp_add(x, x, b); while (MP_NEGP(x));
187 while (MP_CMP(x, >=, b)) x = mp_sub(x, x, b);
190 /* --- Twiddle the signs --- */
192 if (f & f_aneg) { MP_SPLIT(x); x->f ^= MP_NEG; }
193 if (f & f_bneg) { MP_SPLIT(y); y->f ^= MP_NEG; }
195 /* --- Store the results --- */
200 if (*xx) MP_DROP(*xx);
207 if (*yy) MP_DROP(*yy);
213 MP_DROP(X); MP_DROP(Y);
214 MP_DROP(a); MP_DROP(b);
217 /* -- @mp_modinv@ --- *
219 * Arguments: @mp *d@ = destination
223 * Returns: The inverse %$x^{-1} \bmod p$%.
225 * Use: Computes a modular inverse. An assertion fails if %$p$%
229 mp *mp_modinv(mp *d, mp *x, mp *p)
232 mp_gcd(&g, 0, &d, p, x);
233 assert(MP_EQ(g, MP_ONE));
238 /*----- Test rig ----------------------------------------------------------*/
242 static int modinv(dstr *v)
245 mp *x = *(mp **)v[0].buf;
246 mp *m = *(mp **)v[1].buf;
247 mp *r = *(mp **)v[2].buf;
249 mp *y = mp_modinv(MP_NEW, x, m);
251 fputs("\n*** mp_modinv failed", stderr);
252 fputs("\nx = ", stderr); mp_writefile(x, stderr, 10);
253 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
254 fputs("\nexpect = ", stderr); mp_writefile(r, stderr, 10);
255 fputs("\nresult = ", stderr); mp_writefile(y, stderr, 10);
258 MP_DROP(x); MP_DROP(m); MP_DROP(r); MP_DROP(y);
259 assert(mparena_count(MPARENA_GLOBAL) == 0);
263 static int gcd(dstr *v)
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;
272 mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW;
273 mp_gcd(&gg, &xx, &yy, a, b);
275 fputs("\n*** mp_gcd(x) failed", stderr);
276 fputs("\na = ", stderr); mp_writefile(a, stderr, 10);
277 fputs("\nb = ", stderr); mp_writefile(b, stderr, 10);
278 fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 10);
279 fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 10);
284 fputs("\n*** mp_gcd(y) failed", stderr);
285 fputs("\na = ", stderr); mp_writefile(a, stderr, 10);
286 fputs("\nb = ", stderr); mp_writefile(b, stderr, 10);
287 fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 10);
288 fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 10);
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);
298 fputs("\n*** (Alternative result found.)\n", stderr);
304 fputs("\n*** mp_gcd(gcd) failed", stderr);
305 fputs("\na = ", stderr); mp_writefile(a, stderr, 10);
306 fputs("\nb = ", stderr); mp_writefile(b, stderr, 10);
307 fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 10);
308 fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 10);
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);
314 assert(mparena_count(MPARENA_GLOBAL) == 0);
318 static test_chunk tests[] = {
319 { "gcd", gcd, { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
320 { "modinv", modinv, { &type_mp, &type_mp, &type_mp, 0 } },
324 int main(int argc, char *argv[])
327 test_run(argc, argv, tests, SRCDIR "/t/mp");
333 /*----- That's all, folks -------------------------------------------------*/