eigenmath/mgcd.cpp

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