X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/343509982ee8c88ddafd0129b4dcf97e3c7a672d..ceb3f0c0a3b7bb3fa3250d31b04c382894095e52:/gf-gcd.c diff --git a/gf-gcd.c b/gf-gcd.c new file mode 100644 index 0000000..64d61d4 --- /dev/null +++ b/gf-gcd.c @@ -0,0 +1,259 @@ +/* -*-c-*- + * + * $Id: gf-gcd.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $ + * + * Euclidian algorithm on binary polynomials + * + * (c) 2004 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: gf-gcd.c,v $ + * Revision 1.1.2.1 2004/03/21 22:39:46 mdw + * Elliptic curves on binary fields work. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "gf.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @gf_gcd@ --- * + * + * Arguments: @mp **gcd, **xx, **yy@ = where to write the results + * @mp *a, *b@ = sources (must be nonzero) + * + * + * Returns: --- + * + * Use: Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that + * @ax + by = gcd(a, b)@. This is useful for computing modular + * inverses. + */ + +void gf_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b) +{ + mp *x = MP_ONE, *X = MP_ZERO; + mp *y = MP_ZERO, *Y = MP_ONE; + mp *u, *v; + mp *q = MP_NEW; + unsigned f = 0; + +#define f_swap 1u +#define f_ext 2u + + /* --- Sort out some initial flags --- */ + + if (xx || yy) + f |= f_ext; + + /* --- Ensure that @a@ is larger than @b@ --- * + * + * If they're the same length we don't care which order they're in, so this + * unsigned comparison is fine. + */ + + if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) { + { mp *t = a; a = b; b = t; } + f |= f_swap; + } + + /* --- Check for zeroness --- */ + + if (MP_EQ(b, MP_ZERO)) { + + /* --- Store %$|a|$% as the GCD --- */ + + if (gcd) { + if (*gcd) MP_DROP(*gcd); + a = MP_COPY(a); + *gcd = a; + } + + /* --- Store %$1$% and %$0$% in the appropriate bins --- */ + + if (f & f_ext) { + if (f & f_swap) { + mp **t = xx; xx = yy; yy = t; + } + if (xx) { + if (*xx) MP_DROP(*xx); + if (MP_EQ(a, MP_ZERO)) + *xx = MP_ZERO; + else + *xx = MP_ONE; + } + if (yy) { + if (*yy) MP_DROP(*yy); + *yy = MP_ZERO; + } + } + return; + } + + /* --- Take a reference to the arguments --- */ + + a = MP_COPY(a); + b = MP_COPY(b); + + /* --- Make sure @a@ and @b@ are not both even --- */ + + MP_SPLIT(a); a->f &= ~MP_NEG; + MP_SPLIT(b); b->f &= ~MP_NEG; + + u = MP_COPY(a); + v = MP_COPY(b); + + while (MP_LEN(v)) { + mp *t; + gf_div(&q, &u, u, v); + if (f & f_ext) { + t = gf_mul(MP_NEW, X, q); + t = gf_add(t, x, t); + MP_DROP(x); x = X; X = t; + t = gf_mul(MP_NEW, Y, q); + t = gf_add(t, y, t); + MP_DROP(y); y = Y; Y = t; + } + t = u; u = v; v = t; + } + + MP_DROP(q); + if (!gcd) + MP_DROP(u); + else { + if (*gcd) MP_DROP(*gcd); + u->f &= ~MP_NEG; + *gcd = u; + } + + /* --- Perform a little normalization --- */ + + if (f & f_ext) { + + /* --- If @a@ and @b@ got swapped, swap the coefficients back --- */ + + if (f & f_swap) { + mp *t = x; x = y; y = t; + t = a; a = b; b = t; + } + + /* --- Store the results --- */ + + if (!xx) + MP_DROP(x); + else { + if (*xx) MP_DROP(*xx); + *xx = x; + } + + if (!yy) + MP_DROP(y); + else { + if (*yy) MP_DROP(*yy); + *yy = y; + } + } + + MP_DROP(v); + MP_DROP(X); MP_DROP(Y); + MP_DROP(a); MP_DROP(b); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +static int gcd(dstr *v) +{ + int ok = 1; + mp *a = *(mp **)v[0].buf; + mp *b = *(mp **)v[1].buf; + mp *g = *(mp **)v[2].buf; + mp *x = *(mp **)v[3].buf; + mp *y = *(mp **)v[4].buf; + + mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW; + gf_gcd(&gg, &xx, &yy, a, b); + if (!MP_EQ(x, xx)) { + fputs("\n*** mp_gcd(x) failed", stderr); + fputs("\na = ", stderr); mp_writefile(a, stderr, 16); + fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); + fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 16); + fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 16); + fputc('\n', stderr); + ok = 0; + } + if (!MP_EQ(y, yy)) { + fputs("\n*** mp_gcd(y) failed", stderr); + fputs("\na = ", stderr); mp_writefile(a, stderr, 16); + fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); + fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 16); + fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 16); + fputc('\n', stderr); + ok = 0; + } + + if (!ok) { + mp *ax = gf_mul(MP_NEW, a, xx); + mp *by = gf_mul(MP_NEW, b, yy); + ax = gf_add(ax, ax, by); + if (MP_EQ(ax, gg)) + fputs("\n*** (Alternative result found.)\n", stderr); + MP_DROP(ax); + MP_DROP(by); + } + + if (!MP_EQ(g, gg)) { + fputs("\n*** mp_gcd(gcd) failed", stderr); + fputs("\na = ", stderr); mp_writefile(a, stderr, 16); + fputs("\nb = ", stderr); mp_writefile(b, stderr, 16); + fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 16); + fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 16); + fputc('\n', stderr); + ok = 0; + } + MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y); + MP_DROP(gg); MP_DROP(xx); MP_DROP(yy); + assert(mparena_count(MPARENA_GLOBAL) == 0); + return (ok); +} + +static test_chunk tests[] = { + { "gcd", gcd, { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + sub_init(); + test_run(argc, argv, tests, SRCDIR "/tests/gf"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/