git.distorted.org.uk Git - u/mdw/catacomb/blob - mp-gcd.c

   1 /* -*-c-*-
   2  *
   3  * $Id: mp-gcd.c,v 1.5.4.1 2004/03/21 22:39:46 mdw Exp $
   4  *
   5  * Extended GCD calculation
   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: mp-gcd.c,v $
  33  * Revision 1.5.4.1  2004/03/21 22:39:46  mdw
  34  * Elliptic curves on binary fields work.
  35  *
  36  * Revision 1.5  2000/10/08 12:02:41  mdw
  37  * Use Euclid's algorithm rather than the binary one.
  38  *
  39  * Revision 1.4  2000/06/17 11:34:46  mdw
  40  * More hacking for the signs of the outputs.
  41  *
  42  * Revision 1.3  1999/12/10 23:18:39  mdw
  43  * Change interface for suggested destinations.
  44  *
  45  * Revision 1.2  1999/11/22 20:49:56  mdw
  46  * Fix bug which failed to favour `x' when `y' wasn't wanted and the two
  47  * arguments needed swapping.
  48  *
  49  * Revision 1.1  1999/11/17 18:02:16  mdw
  50  * New multiprecision integer arithmetic suite.
  51  *
  52  */
  53
  54 /*----- Header files ------------------------------------------------------*/
  55
  56 #include "mp.h"
  57
  58 /*----- Main code ---------------------------------------------------------*/
  59
  60 /* --- @mp_gcd@ --- *
  61  *
  62  * Arguments:   @mp **gcd, **xx, **yy@ = where to write the results
  63  *              @mp *a, *b@ = sources (must be nonzero)
  64  *
  65  * Returns:     ---
  66  *
  67  * Use:         Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
  68  *              @ax + by = gcd(a, b)@.  This is useful for computing modular
  69  *              inverses.
  70  */
  71
  72 void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
  73 {
  74   mp *x = MP_ONE, *X = MP_ZERO;
  75   mp *y = MP_ZERO, *Y = MP_ONE;
  76   mp *u, *v;
  77   mp *q = MP_NEW;
  78   unsigned f = 0;
  79
  80 #define f_swap 1u
  81 #define f_aneg 2u
  82 #define f_bneg 4u
  83 #define f_ext 8u
  84
  85   /* --- Sort out some initial flags --- */
  86
  87   if (xx || yy)
  88     f |= f_ext;
  89
  90   if (a->f & MP_NEG)
  91     f |= f_aneg;
  92   if (b->f & MP_NEG)
  93     f |= f_bneg;
  94
  95   /* --- Ensure that @a@ is larger than @b@ --- *
  96    *
  97    * Use absolute values here!
  98    */
  99
 100   if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
 101     { mp *t = a; a = b; b = t; }
 102     f |= f_swap;
 103   }
 104
 105   /* --- Check for zeroness --- */
 106
 107   if (MP_EQ(b, MP_ZERO)) {
 108
 109     /* --- Store %$|a|$% as the GCD --- */
 110
 111     if (gcd) {
 112       if (*gcd) MP_DROP(*gcd);
 113       a = MP_COPY(a);
 114       if (a->f & MP_NEG) {
 115         MP_SPLIT(a);
 116         a->f &= ~MP_NEG;
 117         f |= f_aneg;
 118       }
 119       *gcd = a;
 120     }
 121
 122     /* --- Store %$1$% and %$0$% in the appropriate bins --- */
 123
 124     if (f & f_ext) {
 125       if (f & f_swap) {
 126         mp **t = xx; xx = yy; yy = t;
 127       }
 128       if (xx) {
 129         if (*xx) MP_DROP(*xx);
 130         if (MP_EQ(a, MP_ZERO))
 131           *xx = MP_ZERO;
 132         else if (f & f_aneg)
 133           *xx = MP_MONE;
 134         else
 135           *xx = MP_ONE;
 136       }
 137       if (yy) {
 138         if (*yy) MP_DROP(*yy);
 139         *yy = MP_ZERO;
 140       }
 141     }
 142     return;
 143   }
 144
 145   /* --- Take a reference to the arguments --- */
 146
 147   a = MP_COPY(a);
 148   b = MP_COPY(b);
 149
 150   /* --- Make sure @a@ and @b@ are not both even --- */
 151
 152   MP_SPLIT(a); a->f &= ~MP_NEG;
 153   MP_SPLIT(b); b->f &= ~MP_NEG;
 154
 155   u = MP_COPY(a);
 156   v = MP_COPY(b);
 157
 158   while (MP_LEN(v)) {
 159     mp *t;
 160     mp_div(&q, &u, u, v);
 161     if (f & f_ext) {
 162       t = mp_mul(MP_NEW, X, q);
 163       t = mp_sub(t, x, t);
 164       MP_DROP(x); x = X; X = t;
 165       t = mp_mul(MP_NEW, Y, q);
 166       t = mp_sub(t, y, t);
 167       MP_DROP(y); y = Y; Y = t;
 168     }
 169     t = u; u = v; v = t;
 170   }
 171
 172   MP_DROP(q);
 173   if (!gcd)
 174     MP_DROP(u);
 175   else {
 176     if (*gcd) MP_DROP(*gcd);
 177     u->f &= ~MP_NEG;
 178     *gcd = u;
 179   }
 180
 181   /* --- Perform a little normalization --- *
 182    *
 183    * Ensure that the coefficient returned is positive, if there is only one.
 184    * If there are two, favour @y@.  Of course, if the original arguments were
 185    * negative then I'll need to twiddle their signs as well.
 186    */
 187
 188   if (f & f_ext) {
 189
 190     /* --- If @a@ and @b@ got swapped, swap the coefficients back --- */
 191
 192     if (f & f_swap) {
 193       mp *t = x; x = y; y = t;
 194       t = a; a = b; b = t;
 195     }
 196
 197     /* --- Sort out the signs --- *
 198      *
 199      * Note that %$ax + by = a(x - b) + b(y + a)$%.
 200      *
 201      * This is currently bodgy.  It needs sorting out at some time.
 202      */
 203
 204     if (yy) {
 205       if (y->f & MP_NEG) {
 206         do {
 207           y = mp_add(y, y, a);
 208           x = mp_sub(x, x, b);
 209         } while (y->f & MP_NEG);
 210       } else {
 211         while (MP_CMP(y, >=, a)) {
 212           y = mp_sub(y, y, a);
 213           x = mp_add(x, x, b);
 214         }
 215       }
 216     } else {
 217       if (x->f & MP_NEG) {
 218         do
 219           x = mp_add(x, x, b);
 220         while (x->f & MP_NEG);
 221       } else {
 222         while (MP_CMP(x, >=, b))
 223           x = mp_sub(x, x, b);
 224       }
 225     }
 226
 227     /* --- Twiddle the signs --- */
 228
 229     if (f & f_aneg)
 230       x->f ^= MP_NEG;
 231     if (f & f_bneg)
 232       y->f ^= MP_NEG;
 233
 234     /* --- Store the results --- */
 235
 236     if (!xx)
 237       MP_DROP(x);
 238     else {
 239       if (*xx) MP_DROP(*xx);
 240       *xx = x;
 241     }
 242
 243     if (!yy)
 244       MP_DROP(y);
 245     else {
 246       if (*yy) MP_DROP(*yy);
 247       *yy = y;
 248     }
 249   }
 250
 251   MP_DROP(v);
 252   MP_DROP(X); MP_DROP(Y);
 253   MP_DROP(a); MP_DROP(b);
 254 }
 255
 256 /*----- Test rig ----------------------------------------------------------*/
 257
 258 #ifdef TEST_RIG
 259
 260 static int gcd(dstr *v)
 261 {
 262   int ok = 1;
 263   mp *a = *(mp **)v[0].buf;
 264   mp *b = *(mp **)v[1].buf;
 265   mp *g = *(mp **)v[2].buf;
 266   mp *x = *(mp **)v[3].buf;
 267   mp *y = *(mp **)v[4].buf;
 268
 269   mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW;
 270   mp_gcd(&gg, &xx, &yy, a, b);
 271   if (!MP_EQ(x, xx)) {
 272     fputs("\n*** mp_gcd(x) failed", stderr);
 273     fputs("\na      = ", stderr); mp_writefile(a, stderr, 10);
 274     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 10);
 275     fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 10);
 276     fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 10);
 277     fputc('\n', stderr);
 278     ok = 0;
 279   }
 280   if (!MP_EQ(y, yy)) {
 281     fputs("\n*** mp_gcd(y) failed", stderr);
 282     fputs("\na      = ", stderr); mp_writefile(a, stderr, 10);
 283     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 10);
 284     fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 10);
 285     fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 10);
 286     fputc('\n', stderr);
 287     ok = 0;
 288   }
 289
 290   if (!ok) {
 291     mp *ax = mp_mul(MP_NEW, a, xx);
 292     mp *by = mp_mul(MP_NEW, b, yy);
 293     ax = mp_add(ax, ax, by);
 294     if (MP_EQ(ax, gg))
 295       fputs("\n*** (Alternative result found.)\n", stderr);
 296     MP_DROP(ax);
 297     MP_DROP(by);
 298   }
 299
 300   if (!MP_EQ(g, gg)) {
 301     fputs("\n*** mp_gcd(gcd) failed", stderr);
 302     fputs("\na      = ", stderr); mp_writefile(a, stderr, 10);
 303     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 10);
 304     fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 10);
 305     fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 10);
 306     fputc('\n', stderr);
 307     ok = 0;
 308   }
 309   MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y);
 310   MP_DROP(gg); MP_DROP(xx); MP_DROP(yy);
 311   assert(mparena_count(MPARENA_GLOBAL) == 0);
 312   return (ok);
 313 }
 314
 315 static test_chunk tests[] = {
 316   { "gcd", gcd, { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
 317   { 0, 0, { 0 } }
 318 };
 319
 320 int main(int argc, char *argv[])
 321 {
 322   sub_init();
 323   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 324   return (0);
 325 }
 326
 327 #endif
 328
 329 /*----- That's all, folks -------------------------------------------------*/