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

   1 /* -*-c-*-
   2  *
   3  * $Id: mp-gcd.c,v 1.1 1999/11/17 18:02:16 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.1  1999/11/17 18:02:16  mdw
  34  * New multiprecision integer arithmetic suite.
  35  *
  36  */
  37
  38 /*----- Header files ------------------------------------------------------*/
  39
  40 #include "mp.h"
  41
  42 /*----- Main code ---------------------------------------------------------*/
  43
  44 /* --- @mp_gcd@ --- *
  45  *
  46  * Arguments:   @mp **gcd, **xx, **yy@ = where to write the results
  47  *              @mp *a, *b@ = sources (must be nonzero)
  48  *
  49  * Returns:     ---
  50  *
  51  * Use:         Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
  52  *              @ax + by = gcd(a, b)@.  This is useful for computing modular
  53  *              inverses.  Neither @a@ nor @b@ may be zero.  Note that,
  54  *              unlike @mp_div@ for example, it is not possible to specify
  55  *              explicit destinations -- new MPs are always allocated.
  56  */
  57
  58 void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
  59 {
  60   mp *X = MP_ONE, *Y = MP_ZERO;
  61   mp *x = MP_ZERO, *y = MP_ONE;
  62   mp *u, *v;
  63   size_t shift = 0;
  64   int ext = xx || yy;
  65   int swap = 0;
  66
  67   /* --- Ensure that @a@ is larger than @b@ --- */
  68
  69   if (MP_CMP(a, <, b)) {
  70     { mp *t = a; a = b; b = t; }
  71     swap = 1;
  72   }
  73
  74   /* --- Take a reference to the arguments --- */
  75
  76   a = MP_COPY(a);
  77   b = MP_COPY(b);
  78
  79   /* --- Make sure @a@ and @b@ are not both even --- */
  80
  81   if (((a->v[0] | b->v[0]) & 1) == 0) {
  82     mpscan asc, bsc;
  83
  84     /* --- Break off my copies --- */
  85
  86     MP_SPLIT(a);
  87     MP_SPLIT(b);
  88     MP_SCAN(&asc, a);
  89     MP_SCAN(&bsc, b);
  90
  91     /* --- Start scanning --- */
  92
  93     for (;;) {
  94       if (!MP_STEP(&asc) || !MP_STEP(&bsc))
  95         assert(((void)"zero argument passed to mp_gcd", 0));
  96       if (MP_BIT(&asc) || MP_BIT(&bsc))
  97         break;
  98       shift++;
  99     }
 100
 101     /* --- Shift @a@ and @b@ down --- */
 102
 103     a = mp_lsr(a, a, shift);
 104     b = mp_lsr(b, b, shift);
 105   }
 106
 107   /* --- Set up @u@ and @v@ --- */
 108
 109   u = MP_COPY(a);
 110   v = MP_COPY(b);
 111
 112   /* --- Start the main loop --- */
 113
 114   for (;;) {
 115
 116     /* --- While @u@ is even --- */
 117
 118     {
 119       mpscan sc, xsc, ysc;
 120       size_t n = 0, nn = 0;
 121
 122       MP_SCAN(&sc, u);
 123       MP_SCAN(&xsc, X); MP_SCAN(&ysc, Y);
 124       for (;;) {
 125         MP_STEP(&sc);
 126         MP_STEP(&xsc); MP_STEP(&ysc);
 127         if (MP_BIT(&sc))
 128           break;
 129         if (ext && (MP_BIT(&xsc) | MP_BIT(&ysc))) {
 130           if (n) {
 131             X = mp_lsr(X, X, n);
 132             Y = mp_lsr(Y, Y, n);
 133             n = 0;
 134           }
 135           X = mp_add(X, X, b);
 136           Y = mp_sub(Y, Y, a);
 137           MP_SCAN(&xsc, X);
 138           MP_SCAN(&ysc, Y);
 139           MP_STEP(&xsc); MP_STEP(&ysc);
 140         }
 141         n++; nn++;
 142       }
 143
 144       if (nn) {
 145         u = mp_lsr(u, u, nn);
 146         if (ext && n) {
 147           X = mp_lsr(X, X, n);
 148           Y = mp_lsr(Y, Y, n);
 149         }
 150       }
 151     }
 152
 153     /* --- While @v@ is even --- */
 154
 155     {
 156       mpscan sc, xsc, ysc;
 157       size_t n = 0, nn = 0;
 158
 159       MP_SCAN(&sc, v);
 160       MP_SCAN(&xsc, x); MP_SCAN(&ysc, y);
 161       for (;;) {
 162         MP_STEP(&sc);
 163         MP_STEP(&xsc); MP_STEP(&ysc);
 164         if (MP_BIT(&sc))
 165           break;
 166         if (ext && (MP_BIT(&xsc) | MP_BIT(&ysc))) {
 167           if (n) {
 168             x = mp_lsr(x, x, n);
 169             y = mp_lsr(y, y, n);
 170             n = 0;
 171           }
 172           x = mp_add(x, x, b);
 173           y = mp_sub(y, y, a);
 174           MP_SCAN(&xsc, x);
 175           MP_SCAN(&ysc, y);
 176           MP_STEP(&xsc); MP_STEP(&ysc);
 177         }
 178         n++; nn++;
 179       }
 180
 181       if (nn) {
 182         v = mp_lsr(v, v, nn);
 183         if (ext && n) {
 184           x = mp_lsr(x, x, n);
 185           y = mp_lsr(y, y, n);
 186         }
 187       }
 188     }
 189
 190     /* --- End-of-loop fiddling --- */
 191
 192     if (MP_CMP(u, >=, v)) {
 193       u = mp_sub(u, u, v);
 194       if (ext) {
 195         X = mp_sub(X, X, x);
 196         Y = mp_sub(Y, Y, y);
 197       }
 198     } else {
 199       v = mp_sub(v, v, u);
 200       if (ext) {
 201         x = mp_sub(x, x, X);
 202         y = mp_sub(y, y, Y);
 203       }
 204     }
 205
 206     if (MP_CMP(u, ==, MP_ZERO))
 207       break;
 208   }
 209
 210   /* --- Write the results out --- */
 211
 212   if (gcd)
 213     *gcd = mp_lsl(v, v, shift);
 214   else
 215     MP_DROP(v);
 216
 217   /* --- Perform a little normalization --- *
 218    *
 219    * Ensure that the coefficient returned is positive, if there is only one.
 220    * If there are two, favour @y@.
 221    */
 222
 223   if (ext) {
 224     if (swap) {
 225       mp *t = x; x = y; y = t;
 226     }
 227     if (yy) {
 228       if (y->f & MP_NEG) {
 229         y = mp_add(y, y, a);
 230         x = mp_sub(x, x, b);
 231       }
 232     } else if (x->f & MP_NEG)
 233       x = mp_add(x, x, b);
 234
 235     if (xx) *xx = x; else MP_DROP(x);
 236     if (yy) *yy = y; else MP_DROP(y);
 237   }
 238
 239   MP_DROP(u);
 240   MP_DROP(X); MP_DROP(Y);
 241   MP_DROP(a); MP_DROP(b);
 242 }
 243
 244 /*----- Test rig ----------------------------------------------------------*/
 245
 246 #ifdef TEST_RIG
 247
 248 static int gcd(dstr *v)
 249 {
 250   int ok = 1;
 251   mp *a = *(mp **)v[0].buf;
 252   mp *b = *(mp **)v[1].buf;
 253   mp *g = *(mp **)v[2].buf;
 254   mp *x = *(mp **)v[3].buf;
 255   mp *y = *(mp **)v[4].buf;
 256
 257   mp *gg, *xx, *yy;
 258   mp_gcd(&gg, &xx, &yy, a, b);
 259   if (MP_CMP(x, !=, xx)) {
 260     fputs("\n*** mp_gcd(x) failed", stderr);
 261     fputs("\na      = ", stderr); mp_writefile(a, stderr, 10);
 262     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 10);
 263     fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 10);
 264     fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 10);
 265     fputc('\n', stderr);
 266     ok = 0;
 267   }
 268   if (MP_CMP(y, !=, yy)) {
 269     fputs("\n*** mp_gcd(y) failed", stderr);
 270     fputs("\na      = ", stderr); mp_writefile(a, stderr, 10);
 271     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 10);
 272     fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 10);
 273     fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 10);
 274     fputc('\n', stderr);
 275     ok = 0;
 276   }
 277
 278   if (!ok) {
 279     mp *ax = mp_mul(MP_NEW, a, xx);
 280     mp *by = mp_mul(MP_NEW, b, yy);
 281     ax = mp_add(ax, ax, by);
 282     if (MP_CMP(ax, ==, gg))
 283       fputs("\n*** (Alternative result found.)\n", stderr);
 284     MP_DROP(ax);
 285     MP_DROP(by);
 286   }
 287
 288   if (MP_CMP(g, !=, gg)) {
 289     fputs("\n*** mp_gcd(gcd) failed", stderr);
 290     fputs("\na      = ", stderr); mp_writefile(a, stderr, 10);
 291     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 10);
 292     fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 10);
 293     fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 10);
 294     fputc('\n', stderr);
 295     ok = 0;
 296   }
 297   MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y);
 298   MP_DROP(gg); MP_DROP(xx); MP_DROP(yy);
 299   return (ok);
 300 }
 301
 302 static test_chunk tests[] = {
 303   { "gcd", gcd, { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
 304   { 0, 0, { 0 } }
 305 };
 306
 307 int main(int argc, char *argv[])
 308 {
 309   sub_init();
 310   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 311   return (0);
 312 }
 313
 314 #endif
 315
 316 /*----- That's all, folks -------------------------------------------------*/