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

   1 /* -*-c-*-
   2  *
   3  * $Id: mp-gcd.c,v 1.4 2000/06/17 11:34: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.4  2000/06/17 11:34:46  mdw
  34  * More hacking for the signs of the outputs.
  35  *
  36  * Revision 1.3  1999/12/10 23:18:39  mdw
  37  * Change interface for suggested destinations.
  38  *
  39  * Revision 1.2  1999/11/22 20:49:56  mdw
  40  * Fix bug which failed to favour `x' when `y' wasn't wanted and the two
  41  * arguments needed swapping.
  42  *
  43  * Revision 1.1  1999/11/17 18:02:16  mdw
  44  * New multiprecision integer arithmetic suite.
  45  *
  46  */
  47
  48 /*----- Header files ------------------------------------------------------*/
  49
  50 #include "mp.h"
  51
  52 /*----- Main code ---------------------------------------------------------*/
  53
  54 /* --- @mp_gcd@ --- *
  55  *
  56  * Arguments:   @mp **gcd, **xx, **yy@ = where to write the results
  57  *              @mp *a, *b@ = sources (must be nonzero)
  58  *
  59  * Returns:     ---
  60  *
  61  * Use:         Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
  62  *              @ax + by = gcd(a, b)@.  This is useful for computing modular
  63  *              inverses.
  64  */
  65
  66 void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
  67 {
  68   mp *X = MP_ONE, *Y = MP_ZERO;
  69   mp *x = MP_ZERO, *y = MP_ONE;
  70   mp *u, *v;
  71   size_t shift = 0;
  72   unsigned f = 0;
  73
  74   enum {
  75     f_swap = 1u,
  76     f_aneg = 2u,
  77     f_bneg = 4u,
  78     f_ext = 8u
  79   };
  80
  81   /* --- Sort out some initial flags --- */
  82
  83   if (xx || yy)
  84     f |= f_ext;
  85
  86   if (a->f & MP_NEG)
  87     f |= f_aneg;
  88   if (b->f & MP_NEG)
  89     f |= f_bneg;
  90
  91   /* --- Ensure that @a@ is larger than @b@ --- *
  92    *
  93    * Use absolute values here!
  94    */
  95
  96   if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
  97     { mp *t = a; a = b; b = t; }
  98     f |= f_swap;
  99   }
 100
 101   /* --- Check for zeroness --- */
 102
 103   if (MP_CMP(b, ==, MP_ZERO)) {
 104
 105     /* --- Store %$|a|$% as the GCD --- */
 106
 107     if (gcd) {
 108       if (*gcd) MP_DROP(*gcd);
 109       a = MP_COPY(a);
 110       if (a->f & MP_NEG) {
 111         MP_SPLIT(a);
 112         a->f &= ~MP_NEG;
 113         f |= f_aneg;
 114       }
 115       *gcd = a;
 116     }
 117
 118     /* --- Store %$1$% and %$0$% in the appropriate bins --- */
 119
 120     if (f & f_ext) {
 121       if (f & f_swap) {
 122         mp **t = xx; xx = yy; yy = t;
 123       }
 124       if (xx) {
 125         if (*xx) MP_DROP(*xx);
 126         if (MP_CMP(a, ==, MP_ZERO))
 127           *xx = MP_ZERO;
 128         else if (f & f_aneg)
 129           *xx = MP_MONE;
 130         else
 131           *xx = MP_ONE;
 132       }
 133       if (yy) {
 134         if (*yy) MP_DROP(*yy);
 135         *yy = MP_ZERO;
 136       }
 137     }
 138     return;
 139   }
 140
 141   /* --- Take a reference to the arguments --- */
 142
 143   a = MP_COPY(a);
 144   b = MP_COPY(b);
 145
 146   /* --- Make sure @a@ and @b@ are not both even --- */
 147
 148   MP_SPLIT(a); a->f &= ~MP_NEG;
 149   MP_SPLIT(b); b->f &= ~MP_NEG;
 150
 151   if (((a->v[0] | b->v[0]) & 1) == 0) {
 152     mpscan asc, bsc;
 153
 154     /* --- Break off my copies --- */
 155
 156     MP_SCAN(&asc, a);
 157     MP_SCAN(&bsc, b);
 158
 159     /* --- Start scanning --- */
 160
 161     for (;;) {
 162       if (!MP_STEP(&asc) || !MP_STEP(&bsc))
 163         assert(((void)"zero argument passed to mp_gcd", 0));
 164       if (MP_BIT(&asc) || MP_BIT(&bsc))
 165         break;
 166       shift++;
 167     }
 168
 169     /* --- Shift @a@ and @b@ down --- */
 170
 171     a = mp_lsr(a, a, shift);
 172     b = mp_lsr(b, b, shift);
 173   }
 174
 175   /* --- Set up @u@ and @v@ --- */
 176
 177   u = MP_COPY(a);
 178   v = MP_COPY(b);
 179
 180   /* --- Start the main loop --- */
 181
 182   for (;;) {
 183
 184     /* --- While @u@ is even --- */
 185
 186     {
 187       mpscan sc, xsc, ysc;
 188       size_t n = 0, nn = 0;
 189
 190       MP_SCAN(&sc, u);
 191       MP_SCAN(&xsc, X); MP_SCAN(&ysc, Y);
 192       for (;;) {
 193         MP_STEP(&sc);
 194         MP_STEP(&xsc); MP_STEP(&ysc);
 195         if (MP_BIT(&sc))
 196           break;
 197         if ((f & f_ext) && (MP_BIT(&xsc) | MP_BIT(&ysc))) {
 198           if (n) {
 199             X = mp_lsr(X, X, n);
 200             Y = mp_lsr(Y, Y, n);
 201             n = 0;
 202           }
 203           X = mp_add(X, X, b);
 204           Y = mp_sub(Y, Y, a);
 205           MP_SCAN(&xsc, X);
 206           MP_SCAN(&ysc, Y);
 207           MP_STEP(&xsc); MP_STEP(&ysc);
 208         }
 209         n++; nn++;
 210       }
 211
 212       if (nn) {
 213         u = mp_lsr(u, u, nn);
 214         if ((f & f_ext) && n) {
 215           X = mp_lsr(X, X, n);
 216           Y = mp_lsr(Y, Y, n);
 217         }
 218       }
 219     }
 220
 221     /* --- While @v@ is even --- */
 222
 223     {
 224       mpscan sc, xsc, ysc;
 225       size_t n = 0, nn = 0;
 226
 227       MP_SCAN(&sc, v);
 228       MP_SCAN(&xsc, x); MP_SCAN(&ysc, y);
 229       for (;;) {
 230         MP_STEP(&sc);
 231         MP_STEP(&xsc); MP_STEP(&ysc);
 232         if (MP_BIT(&sc))
 233           break;
 234         if ((f & f_ext) && (MP_BIT(&xsc) | MP_BIT(&ysc))) {
 235           if (n) {
 236             x = mp_lsr(x, x, n);
 237             y = mp_lsr(y, y, n);
 238             n = 0;
 239           }
 240           x = mp_add(x, x, b);
 241           y = mp_sub(y, y, a);
 242           MP_SCAN(&xsc, x);
 243           MP_SCAN(&ysc, y);
 244           MP_STEP(&xsc); MP_STEP(&ysc);
 245         }
 246         n++; nn++;
 247       }
 248
 249       if (nn) {
 250         v = mp_lsr(v, v, nn);
 251         if ((f & f_ext) && n) {
 252           x = mp_lsr(x, x, n);
 253           y = mp_lsr(y, y, n);
 254         }
 255       }
 256     }
 257
 258     /* --- End-of-loop fiddling --- */
 259
 260     if (MP_CMP(u, >=, v)) {
 261       u = mp_sub(u, u, v);
 262       if (f & f_ext) {
 263         X = mp_sub(X, X, x);
 264         Y = mp_sub(Y, Y, y);
 265       }
 266     } else {
 267       v = mp_sub(v, v, u);
 268       if (f & f_ext) {
 269         x = mp_sub(x, x, X);
 270         y = mp_sub(y, y, Y);
 271       }
 272     }
 273
 274     if (MP_CMP(u, ==, MP_ZERO))
 275       break;
 276   }
 277
 278   /* --- Write the results out --- */
 279
 280   if (!gcd)
 281     MP_DROP(v);
 282   else {
 283     if (*gcd) MP_DROP(*gcd);
 284     *gcd = mp_lsl(v, v, shift);
 285   }
 286
 287   /* --- Perform a little normalization --- *
 288    *
 289    * Ensure that the coefficient returned is positive, if there is only one.
 290    * If there are two, favour @y@.  Of course, if the original arguments were
 291    * negative then I'll need to twiddle their signs as well.
 292    */
 293
 294   if (f & f_ext) {
 295
 296     /* --- If @a@ and @b@ got swapped, swap the coefficients back --- */
 297
 298     if (f & f_swap) {
 299       mp *t = x; x = y; y = t;
 300       t = a; a = b; b = t;
 301     }
 302
 303     /* --- Sort out the signs --- *
 304      *
 305      * Note that %$ax + by = a(x - b) + b(y + a)$%.
 306      *
 307      * This is currently bodgy.  It needs sorting out at some time.
 308      */
 309
 310     if (yy) {
 311       if (y->f & MP_NEG) {
 312         do {
 313           y = mp_add(y, y, a);
 314           x = mp_sub(x, x, b);
 315         } while (y->f & MP_NEG);
 316       } else {
 317         while (MP_CMP(y, >=, a)) {
 318           y = mp_sub(y, y, a);
 319           x = mp_add(x, x, b);
 320         }
 321       }
 322     } else {
 323       if (x->f & MP_NEG) {
 324         do
 325           x = mp_add(x, x, b);
 326         while (x->f & MP_NEG);
 327       } else {
 328         while (MP_CMP(x, >=, b))
 329           x = mp_sub(x, x, b);
 330       }
 331     }
 332
 333     /* --- Twiddle the signs --- */
 334
 335     if (f & f_aneg)
 336       x->f ^= MP_NEG;
 337     if (f & f_bneg)
 338       y->f ^= MP_NEG;
 339
 340     /* --- Store the results --- */
 341
 342     if (!xx)
 343       MP_DROP(x);
 344     else {
 345       if (*xx) MP_DROP(*xx);
 346       *xx = x;
 347     }
 348
 349     if (!yy)
 350       MP_DROP(y);
 351     else {
 352       if (*yy) MP_DROP(*yy);
 353       *yy = y;
 354     }
 355   }
 356
 357   MP_DROP(u);
 358   MP_DROP(X); MP_DROP(Y);
 359   MP_DROP(a); MP_DROP(b);
 360 }
 361
 362 /*----- Test rig ----------------------------------------------------------*/
 363
 364 #ifdef TEST_RIG
 365
 366 static int gcd(dstr *v)
 367 {
 368   int ok = 1;
 369   mp *a = *(mp **)v[0].buf;
 370   mp *b = *(mp **)v[1].buf;
 371   mp *g = *(mp **)v[2].buf;
 372   mp *x = *(mp **)v[3].buf;
 373   mp *y = *(mp **)v[4].buf;
 374
 375   mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW;
 376   mp_gcd(&gg, &xx, &yy, a, b);
 377   if (MP_CMP(x, !=, xx)) {
 378     fputs("\n*** mp_gcd(x) failed", stderr);
 379     fputs("\na      = ", stderr); mp_writefile(a, stderr, 10);
 380     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 10);
 381     fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 10);
 382     fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 10);
 383     fputc('\n', stderr);
 384     ok = 0;
 385   }
 386   if (MP_CMP(y, !=, yy)) {
 387     fputs("\n*** mp_gcd(y) failed", stderr);
 388     fputs("\na      = ", stderr); mp_writefile(a, stderr, 10);
 389     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 10);
 390     fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 10);
 391     fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 10);
 392     fputc('\n', stderr);
 393     ok = 0;
 394   }
 395
 396   if (!ok) {
 397     mp *ax = mp_mul(MP_NEW, a, xx);
 398     mp *by = mp_mul(MP_NEW, b, yy);
 399     ax = mp_add(ax, ax, by);
 400     if (MP_CMP(ax, ==, gg))
 401       fputs("\n*** (Alternative result found.)\n", stderr);
 402     MP_DROP(ax);
 403     MP_DROP(by);
 404   }
 405
 406   if (MP_CMP(g, !=, gg)) {
 407     fputs("\n*** mp_gcd(gcd) failed", stderr);
 408     fputs("\na      = ", stderr); mp_writefile(a, stderr, 10);
 409     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 10);
 410     fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 10);
 411     fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 10);
 412     fputc('\n', stderr);
 413     ok = 0;
 414   }
 415   MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y);
 416   MP_DROP(gg); MP_DROP(xx); MP_DROP(yy);
 417   assert(mparena_count(MPARENA_GLOBAL) == 0);
 418   return (ok);
 419 }
 420
 421 static test_chunk tests[] = {
 422   { "gcd", gcd, { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
 423   { 0, 0, { 0 } }
 424 };
 425
 426 int main(int argc, char *argv[])
 427 {
 428   sub_init();
 429   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 430   return (0);
 431 }
 432
 433 #endif
 434
 435 /*----- That's all, folks -------------------------------------------------*/