chiark - git - mdw - catacomb/blob - mp-gcd.c

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