mgcd.cpp

   1 //-----------------------------------------------------------------------------
   2 //
   3 //      Bignum GCD
   4 //
   5 //      Uses the binary GCD algorithm.
   6 //
   7 //      See "The Art of Computer Programming" p. 338.
   8 //
   9 //      mgcd always returns a positive value
  10 //
  11 //      mgcd(0, 0) = 0
  12 //
  13 //      mgcd(u, 0) = |u|
  14 //
  15 //      mgcd(0, v) = |v|
  16 //
  17 //-----------------------------------------------------------------------------
  18
  19 #include "stdafx.h"
  20 #include "defs.h"
  21
  22 unsigned int *
  23 mgcd(unsigned int *u, unsigned int *v)
  24 {
  25         int i, k, n;
  26         unsigned int *t;
  27
  28         if (MZERO(u)) {
  29                 t = mcopy(v);
  30                 MSIGN(t) = 1;
  31                 return t;
  32         }
  33
  34         if (MZERO(v)) {
  35                 t = mcopy(u);
  36                 MSIGN(t) = 1;
  37                 return t;
  38         }
  39
  40         u = mcopy(u);
  41         v = mcopy(v);
  42
  43         MSIGN(u) = 1;
  44         MSIGN(v) = 1;
  45
  46         k = 0;
  47
  48         while ((u[0] & 1) == 0 && (v[0] & 1) == 0) {
  49                 mshiftright(u);
  50                 mshiftright(v);
  51                 k++;
  52         }
  53
  54         if (u[0] & 1) {
  55                 t = mcopy(v);
  56                 MSIGN(t) *= -1;
  57         } else
  58                 t = mcopy(u);
  59
  60         while (1) {
  61
  62                 while ((t[0] & 1) == 0)
  63                         mshiftright(t);
  64
  65                 if (MSIGN(t) == 1) {
  66                         mfree(u);
  67                         u = mcopy(t);
  68                 } else {
  69                         mfree(v);
  70                         v = mcopy(t);
  71                         MSIGN(v) *= -1;
  72                 }
  73
  74                 mfree(t);
  75
  76                 t = msub(u, v);
  77
  78                 if (MZERO(t)) {
  79                         mfree(t);
  80                         mfree(v);
  81                         n = (k / 32) + 1;
  82                         v = mnew(n);
  83                         MSIGN(v) = 1;
  84                         MLENGTH(v) = n;
  85                         for (i = 0; i < n; i++)
  86                                 v[i] = 0;
  87                         mp_set_bit(v, k);
  88                         t = mmul(u, v);
  89                         mfree(u);
  90                         mfree(v);
  91                         return t;
  92                 }
  93         }
  94 }