156 lines
2.2 KiB
C++
156 lines
2.2 KiB
C++
//-----------------------------------------------------------------------------
|
|
//
|
|
// Bignum GCD
|
|
//
|
|
// Uses the binary GCD algorithm.
|
|
//
|
|
// See "The Art of Computer Programming" p. 338.
|
|
//
|
|
// mgcd always returns a positive value
|
|
//
|
|
// mgcd(0, 0) = 0
|
|
//
|
|
// mgcd(u, 0) = |u|
|
|
//
|
|
// mgcd(0, v) = |v|
|
|
//
|
|
//-----------------------------------------------------------------------------
|
|
|
|
#include "stdafx.h"
|
|
#include "defs.h"
|
|
|
|
unsigned int *
|
|
mgcd(unsigned int *u, unsigned int *v)
|
|
{
|
|
int i, k, n;
|
|
unsigned int *t;
|
|
|
|
if (MZERO(u)) {
|
|
t = mcopy(v);
|
|
MSIGN(t) = 1;
|
|
return t;
|
|
}
|
|
|
|
if (MZERO(v)) {
|
|
t = mcopy(u);
|
|
MSIGN(t) = 1;
|
|
return t;
|
|
}
|
|
|
|
u = mcopy(u);
|
|
v = mcopy(v);
|
|
|
|
MSIGN(u) = 1;
|
|
MSIGN(v) = 1;
|
|
|
|
k = 0;
|
|
|
|
while ((u[0] & 1) == 0 && (v[0] & 1) == 0) {
|
|
mshiftright(u);
|
|
mshiftright(v);
|
|
k++;
|
|
}
|
|
|
|
if (u[0] & 1) {
|
|
t = mcopy(v);
|
|
MSIGN(t) *= -1;
|
|
} else
|
|
t = mcopy(u);
|
|
|
|
while (1) {
|
|
|
|
while ((t[0] & 1) == 0)
|
|
mshiftright(t);
|
|
|
|
if (MSIGN(t) == 1) {
|
|
mfree(u);
|
|
u = mcopy(t);
|
|
} else {
|
|
mfree(v);
|
|
v = mcopy(t);
|
|
MSIGN(v) *= -1;
|
|
}
|
|
|
|
mfree(t);
|
|
|
|
t = msub(u, v);
|
|
|
|
if (MZERO(t)) {
|
|
mfree(t);
|
|
mfree(v);
|
|
n = (k / 32) + 1;
|
|
v = mnew(n);
|
|
MSIGN(v) = 1;
|
|
MLENGTH(v) = n;
|
|
for (i = 0; i < n; i++)
|
|
v[i] = 0;
|
|
mp_set_bit(v, k);
|
|
t = mmul(u, v);
|
|
mfree(u);
|
|
mfree(v);
|
|
return t;
|
|
}
|
|
}
|
|
}
|
|
|
|
#if SELFTEST
|
|
|
|
static unsigned int *egcd(unsigned int *, unsigned int *);
|
|
|
|
void
|
|
test_mgcd(void)
|
|
{
|
|
int i, j, n;
|
|
unsigned int *a, *b, *c, *d;
|
|
logout("testing mgcd\n");
|
|
n = mtotal;
|
|
for (i = 1; i < 100; i++) {
|
|
a = mint(i);
|
|
for (j = 1; j < 100; j++) {
|
|
b = mint(j);
|
|
c = mgcd(a, b);
|
|
d = egcd(a, b);
|
|
if (mcmp(c, d) != 0) {
|
|
logout("failed\n");
|
|
errout();
|
|
}
|
|
mfree(b);
|
|
mfree(c);
|
|
mfree(d);
|
|
}
|
|
mfree(a);
|
|
}
|
|
if (n != mtotal) {
|
|
logout("memory leak\n");
|
|
errout();
|
|
}
|
|
logout("ok\n");
|
|
}
|
|
|
|
// Euclid's algorithm
|
|
|
|
static unsigned int *
|
|
egcd(unsigned int *a, unsigned int *b)
|
|
{
|
|
int sign;
|
|
unsigned int *c;
|
|
if (MZERO(b))
|
|
stop("divide by zero");
|
|
b = mcopy(b);
|
|
if (MZERO(a))
|
|
return b;
|
|
sign = MSIGN(b);
|
|
a = mcopy(a);
|
|
while (!MZERO(b)) {
|
|
c = mmod(a, b);
|
|
mfree(a);
|
|
a = b;
|
|
b = c;
|
|
}
|
|
mfree(b);
|
|
MSIGN(a) = sign;
|
|
return a;
|
|
}
|
|
|
|
#endif
|