4.4BSD/usr/src/contrib/calc-1.26.4/zfunc.c
/*
* Copyright (c) 1993 David I. Bell
* Permission is granted to use, distribute, or modify this source,
* provided that this copyright notice remains intact.
*
* Extended precision integral arithmetic non-primitive routines
*/
#include "math.h"
static ZVALUE primeprod; /* product of primes under 100 */
ZVALUE _tenpowers_[32]; /* table of 10^2^n */
#if 0
static char *abortmsg = "Calculation aborted";
static char *memmsg = "Not enough memory";
#endif
/*
* Compute the factorial of a number.
*/
void
zfact(z, dest)
ZVALUE z, *dest;
{
long ptwo; /* count of powers of two */
long n; /* current multiplication value */
long m; /* reduced multiplication value */
long mul; /* collected value to multiply by */
ZVALUE res, temp;
if (isneg(z))
error("Negative argument for factorial");
if (isbig(z))
error("Very large factorial");
n = (istiny(z) ? z1tol(z) : z2tol(z));
ptwo = 0;
mul = 1;
res = _one_;
/*
* Multiply numbers together, but squeeze out all powers of two.
* We will put them back in at the end. Also collect multiple
* numbers together until there is a risk of overflow.
*/
for (; n > 1; n--) {
for (m = n; ((m & 0x1) == 0); m >>= 1)
ptwo++;
mul *= m;
if (mul < BASE1/2)
continue;
zmuli(res, mul, &temp);
freeh(res.v);
res = temp;
mul = 1;
}
/*
* Multiply by the remaining value, then scale result by
* the proper power of two.
*/
if (mul > 1) {
zmuli(res, mul, &temp);
freeh(res.v);
res = temp;
}
zshift(res, ptwo, &temp);
freeh(res.v);
*dest = temp;
}
/*
* Compute the product of the primes up to the specified number.
*/
void
zpfact(z, dest)
ZVALUE z, *dest;
{
long n; /* limiting number to multiply by */
long p; /* current prime */
long i; /* test value */
long mul; /* collected value to multiply by */
ZVALUE res, temp;
if (isneg(z))
error("Negative argument for factorial");
if (isbig(z))
error("Very large factorial");
n = (istiny(z) ? z1tol(z) : z2tol(z));
/*
* Multiply by the primes in order, collecting multiple numbers
* together until there is a change of overflow.
*/
mul = 1 + (n > 1);
res = _one_;
for (p = 3; p <= n; p += 2) {
for (i = 3; (i * i) <= p; i += 2) {
if ((p % i) == 0)
goto next;
}
mul *= p;
if (mul < BASE1/2)
continue;
zmuli(res, mul, &temp);
freeh(res.v);
res = temp;
mul = 1;
next: ;
}
/*
* Multiply by the final value if any.
*/
if (mul > 1) {
zmuli(res, mul, &temp);
freeh(res.v);
res = temp;
}
*dest = res;
}
/*
* Compute the least common multiple of all the numbers up to the
* specified number.
*/
void
zlcmfact(z, dest)
ZVALUE z, *dest;
{
long n; /* limiting number to multiply by */
long p; /* current prime */
long pp; /* power of prime */
long i; /* test value */
ZVALUE res, temp;
if (isneg(z) || iszero(z))
error("Non-positive argument for lcmfact");
if (isbig(z))
error("Very large lcmfact");
n = (istiny(z) ? z1tol(z) : z2tol(z));
/*
* Multiply by powers of the necessary odd primes in order.
* The power for each prime is the highest one which is not
* more than the specified number.
*/
res = _one_;
for (p = 3; p <= n; p += 2) {
for (i = 3; (i * i) <= p; i += 2) {
if ((p % i) == 0)
goto next;
}
i = p;
while (i <= n) {
pp = i;
i *= p;
}
zmuli(res, pp, &temp);
freeh(res.v);
res = temp;
next: ;
}
/*
* Finish by scaling by the necessary power of two.
*/
zshift(res, zhighbit(z), dest);
freeh(res.v);
}
/*
* Compute the permuation function M! / (M - N)!.
*/
void
zperm(z1, z2, res)
ZVALUE z1, z2, *res;
{
long count;
ZVALUE cur, tmp, ans;
if (isneg(z1) || isneg(z2))
error("Negative argument for permutation");
if (zrel(z1, z2) < 0)
error("Second arg larger than first in permutation");
if (isbig(z2))
error("Very large permutation");
count = (istiny(z2) ? z1tol(z2) : z2tol(z2));
zcopy(z1, &ans);
zsub(z1, _one_, &cur);
while (--count > 0) {
zmul(ans, cur, &tmp);
freeh(ans.v);
ans = tmp;
zsub(cur, _one_, &tmp);
freeh(cur.v);
cur = tmp;
}
freeh(cur.v);
*res = ans;
}
/*
* Compute the combinatorial function M! / ( N! * (M - N)! ).
*/
void
zcomb(z1, z2, res)
ZVALUE z1, z2, *res;
{
ZVALUE ans;
ZVALUE mul, div, temp;
FULL count, i;
HALF dh[2];
if (isneg(z1) || isneg(z2))
error("Negative argument for combinatorial");
zsub(z1, z2, &temp);
if (isneg(temp)) {
freeh(temp.v);
error("Second arg larger than first for combinatorial");
}
if (isbig(z2) && isbig(temp)) {
freeh(temp.v);
error("Very large combinatorial");
}
count = (istiny(z2) ? z1tol(z2) : z2tol(z2));
i = (istiny(temp) ? z1tol(temp) : z2tol(temp));
if (isbig(z2) || (!isbig(temp) && (i < count)))
count = i;
freeh(temp.v);
mul = z1;
div.sign = 0;
div.v = dh;
ans = _one_;
for (i = 1; i <= count; i++) {
dh[0] = i & BASE1;
dh[1] = i / BASE;
div.len = 1 + (dh[1] != 0);
zmul(ans, mul, &temp);
freeh(ans.v);
zquo(temp, div, &ans);
freeh(temp.v);
zsub(mul, _one_, &temp);
if (mul.v != z1.v)
freeh(mul.v);
mul = temp;
}
if (mul.v != z1.v)
freeh(mul.v);
*res = ans;
}
/*
* Perform a probabilistic primality test (algorithm P in Knuth).
* Returns FALSE if definitely not prime, or TRUE if probably prime.
* Count determines how many times to check for primality.
* The chance of a non-prime passing this test is less than (1/4)^count.
* For example, a count of 100 fails for only 1 in 10^60 numbers.
*/
BOOL
zprimetest(z, count)
ZVALUE z; /* number to test for primeness */
long count;
{
long ij, ik, ix;
ZVALUE zm1, z1, z2, z3, ztmp;
HALF val[2];
z.sign = 0;
if (iseven(z)) /* if even, not prime if not 2 */
return (istwo(z) != 0);
/*
* See if the number is small, and is either a small prime,
* or is divisible by a small prime.
*/
if (istiny(z) && (*z.v <= (HALF)(101*101-1))) {
ix = *z.v;
for (ik = 3; (ik <= 97) && ((ik * ik) <= ix); ik += 2)
if ((ix % ik) == 0)
return FALSE;
return TRUE;
}
/*
* See if the number is divisible by one of the primes 3, 5,
* 7, 11, or 13. This is a very easy check.
*/
ij = zmodi(z, 15015L);
if (!(ij % 3) || !(ij % 5) || !(ij % 7) || !(ij % 11) || !(ij % 13))
return FALSE;
/*
* Check the gcd of the number and the product of more of the first
* few odd primes. We must build the prime product on the first call.
*/
ztmp.sign = 0;
ztmp.len = 1;
ztmp.v = val;
if (primeprod.len == 0) {
val[0] = 101;
zpfact(ztmp, &primeprod);
}
zgcd(z, primeprod, &z1);
if (!isunit(z1)) {
freeh(z1.v);
return FALSE;
}
freeh(z1.v);
/*
* Not divisible by a small prime, so onward with the real test.
* Make sure the count is limited by the number of odd numbers between
* three and the number being tested.
*/
ix = ((istiny(z) ? z1tol(z) : z2tol(z) - 3) / 2);
if (count > ix) count = ix;
zsub(z, _one_, &zm1);
ik = zlowbit(zm1);
zshift(zm1, -ik, &z1);
/*
* Loop over various "random" numbers, testing each one.
* These numbers are the odd numbers starting from three.
*/
for (ix = 0; ix < count; ix++) {
val[0] = (ix * 2) + 3;
ij = 0;
zpowermod(ztmp, z1, z, &z3);
for (;;) {
if (isone(z3)) {
if (ij) /* number is definitely not prime */
goto notprime;
break;
}
if (zcmp(z3, zm1) == 0)
break;
if (++ij >= ik)
goto notprime; /* number is definitely not prime */
zsquare(z3, &z2);
freeh(z3.v);
zmod(z2, z, &z3);
freeh(z2.v);
}
freeh(z3.v);
}
freeh(zm1.v);
freeh(z1.v);
return TRUE; /* number might be prime */
notprime:
freeh(z3.v);
freeh(zm1.v);
freeh(z1.v);
return FALSE;
}
/*
* Compute the Jacobi function (p / q) for odd q.
* If q is prime then the result is:
* 1 if p == x^2 (mod q) for some x.
* -1 otherwise.
* If q is not prime, then the result is not meaningful if it is 1.
* This function returns 0 if q is even or q < 0.
*/
FLAG
zjacobi(z1, z2)
ZVALUE z1, z2;
{
ZVALUE p, q, tmp;
long lowbit;
int val;
if (iseven(z2) || isneg(z2))
return 0;
val = 1;
if (iszero(z1) || isone(z1))
return val;
if (isunit(z1)) {
if ((*z2.v - 1) & 0x2)
val = -val;
return val;
}
zcopy(z1, &p);
zcopy(z2, &q);
for (;;) {
zmod(p, q, &tmp);
freeh(p.v);
p = tmp;
if (iszero(p)) {
freeh(p.v);
p = _one_;
}
if (iseven(p)) {
lowbit = zlowbit(p);
zshift(p, -lowbit, &tmp);
freeh(p.v);
p = tmp;
if ((lowbit & 1) && (((*q.v & 0x7) == 3) || ((*q.v & 0x7) == 5)))
val = -val;
}
if (isunit(p)) {
freeh(p.v);
freeh(q.v);
return val;
}
if ((*p.v & *q.v & 0x3) == 3)
val = -val;
tmp = q;
q = p;
p = tmp;
}
}
/*
* Return the Fibonacci number F(n).
* This is evaluated by recursively using the formulas:
* F(2N+1) = F(N+1)^2 + F(N)^2
* and
* F(2N) = F(N+1)^2 - F(N-1)^2
*/
void
zfib(z, res)
ZVALUE z, *res;
{
unsigned long i;
long n;
int sign;
ZVALUE fnm1, fn, fnp1; /* consecutive fibonacci values */
ZVALUE t1, t2, t3;
if (isbig(z))
error("Very large Fibonacci number");
n = (istiny(z) ? z1tol(z) : z2tol(z));
if (n == 0) {
*res = _zero_;
return;
}
sign = z.sign && ((n & 0x1) == 0);
if (n <= 2) {
*res = _one_;
res->sign = (BOOL)sign;
return;
}
i = TOPFULL;
while ((i & n) == 0)
i >>= 1L;
i >>= 1L;
fnm1 = _zero_;
fn = _one_;
fnp1 = _one_;
while (i) {
zsquare(fnm1, &t1);
zsquare(fn, &t2);
zsquare(fnp1, &t3);
freeh(fnm1.v);
freeh(fn.v);
freeh(fnp1.v);
zadd(t2, t3, &fnp1);
zsub(t3, t1, &fn);
freeh(t1.v);
freeh(t2.v);
freeh(t3.v);
if (i & n) {
fnm1 = fn;
fn = fnp1;
zadd(fnm1, fn, &fnp1);
} else
zsub(fnp1, fn, &fnm1);
i >>= 1L;
}
freeh(fnm1.v);
freeh(fnp1.v);
*res = fn;
res->sign = (BOOL)sign;
}
/*
* Compute the result of raising one number to the power of another
* The second number is assumed to be non-negative.
* It cannot be too large except for trivial cases.
*/
void
zpowi(z1, z2, res)
ZVALUE z1, z2, *res;
{
int sign; /* final sign of number */
unsigned long power; /* power to raise to */
unsigned long bit; /* current bit value */
long twos; /* count of times 2 is in result */
ZVALUE ans, temp;
sign = (z1.sign && isodd(z2));
z1.sign = 0;
z2.sign = 0;
if (iszero(z2)) { /* number raised to power 0 */
if (iszero(z1))
error("Zero raised to zero power");
*res = _one_;
return;
}
if (isleone(z1)) { /* 0, 1, or -1 raised to a power */
ans = _one_;
ans.sign = (BOOL)sign;
if (*z1.v == 0)
ans = _zero_;
*res = ans;
return;
}
if (isbig(z2))
error("Raising to very large power");
power = (istiny(z2) ? z1tol(z2) : z2tol(z2));
if (istwo(z1)) { /* two raised to a power */
zbitvalue((long) power, res);
return;
}
/*
* See if this is a power of ten
*/
if (istiny(z1) && (*z1.v == 10)) {
ztenpow((long) power, res);
res->sign = (BOOL)sign;
return;
}
/*
* Handle low powers specially
*/
if (power <= 4) {
switch ((int) power) {
case 1:
ans.len = z1.len;
ans.v = alloc(ans.len);
copyval(z1, ans);
ans.sign = (BOOL)sign;
*res = ans;
return;
case 2:
zsquare(z1, res);
return;
case 3:
zsquare(z1, &temp);
zmul(z1, temp, res);
freeh(temp.v);
res->sign = (BOOL)sign;
return;
case 4:
zsquare(z1, &temp);
zsquare(temp, res);
freeh(temp.v);
return;
}
}
/*
* Shift out all powers of twos so the multiplies are smaller.
* We will shift back the right amount when done.
*/
twos = 0;
if (iseven(z1)) {
twos = zlowbit(z1);
ans.v = alloc(z1.len);
ans.len = z1.len;
copyval(z1, ans);
shiftr(ans, twos);
trim(&ans);
z1 = ans;
twos *= power;
}
/*
* Compute the power by squaring and multiplying.
* This uses the left to right method of power raising.
*/
bit = TOPFULL;
while ((bit & power) == 0)
bit >>= 1L;
bit >>= 1L;
zsquare(z1, &ans);
if (bit & power) {
zmul(ans, z1, &temp);
freeh(ans.v);
ans = temp;
}
bit >>= 1L;
while (bit) {
zsquare(ans, &temp);
freeh(ans.v);
ans = temp;
if (bit & power) {
zmul(ans, z1, &temp);
freeh(ans.v);
ans = temp;
}
bit >>= 1L;
}
/*
* Scale back up by proper power of two
*/
if (twos) {
zshift(ans, twos, &temp);
freeh(ans.v);
ans = temp;
freeh(z1.v);
}
ans.sign = (BOOL)sign;
*res = ans;
}
/*
* Compute ten to the specified power
* This saves some work since the squares of ten are saved.
*/
void
ztenpow(power, res)
long power;
ZVALUE *res;
{
long i;
ZVALUE ans;
ZVALUE temp;
if (power <= 0) {
*res = _one_;
return;
}
ans = _one_;
_tenpowers_[0] = _ten_;
for (i = 0; power; i++) {
if (_tenpowers_[i].len == 0)
zsquare(_tenpowers_[i-1], &_tenpowers_[i]);
if (power & 0x1) {
zmul(ans, _tenpowers_[i], &temp);
freeh(ans.v);
ans = temp;
}
power /= 2;
}
*res = ans;
}
/*
* Calculate modular inverse suppressing unnecessary divisions.
* This is based on the Euclidian algorithm for large numbers.
* (Algorithm X from Knuth Vol 2, section 4.5.2. and exercise 17)
* Returns TRUE if there is no solution because the numbers
* are not relatively prime.
*/
BOOL
zmodinv(u, v, res)
ZVALUE u, v;
ZVALUE *res;
{
FULL q1, q2, ui3, vi3, uh, vh, A, B, C, D, T;
ZVALUE u2, u3, v2, v3, qz, tmp1, tmp2, tmp3;
if (isneg(u) || isneg(v) || (zrel(u, v) >= 0))
zmod(u, v, &v3);
else
zcopy(u, &v3);
zcopy(v, &u3);
u2 = _zero_;
v2 = _one_;
/*
* Loop here while the size of the numbers remain above
* the size of a FULL. Throughout this loop u3 >= v3.
*/
while ((u3.len > 1) && !iszero(v3)) {
uh = (((FULL) u3.v[u3.len - 1]) << BASEB) + u3.v[u3.len - 2];
vh = 0;
if ((v3.len + 1) >= u3.len)
vh = v3.v[v3.len - 1];
if (v3.len == u3.len)
vh = (vh << BASEB) + v3.v[v3.len - 2];
A = 1;
B = 0;
C = 0;
D = 1;
/*
* Calculate successive quotients of the continued fraction
* expansion using only single precision arithmetic until
* greater precision is required.
*/
while ((vh + C) && (vh + D)) {
q1 = (uh + A) / (vh + C);
q2 = (uh + B) / (vh + D);
if (q1 != q2)
break;
T = A - q1 * C;
A = C;
C = T;
T = B - q1 * D;
B = D;
D = T;
T = uh - q1 * vh;
uh = vh;
vh = T;
}
/*
* If B is zero, then we made no progress because
* the calculation requires a very large quotient.
* So we must do this step of the calculation in
* full precision
*/
if (B == 0) {
zquo(u3, v3, &qz);
zmul(qz, v2, &tmp1);
zsub(u2, tmp1, &tmp2);
freeh(tmp1.v);
freeh(u2.v);
u2 = v2;
v2 = tmp2;
zmul(qz, v3, &tmp1);
zsub(u3, tmp1, &tmp2);
freeh(tmp1.v);
freeh(u3.v);
u3 = v3;
v3 = tmp2;
freeh(qz.v);
continue;
}
/*
* Apply the calculated A,B,C,D numbers to the current
* values to update them as if the full precision
* calculations had been carried out.
*/
zmuli(u2, (long) A, &tmp1);
zmuli(v2, (long) B, &tmp2);
zadd(tmp1, tmp2, &tmp3);
freeh(tmp1.v);
freeh(tmp2.v);
zmuli(u2, (long) C, &tmp1);
zmuli(v2, (long) D, &tmp2);
freeh(u2.v);
freeh(v2.v);
u2 = tmp3;
zadd(tmp1, tmp2, &v2);
freeh(tmp1.v);
freeh(tmp2.v);
zmuli(u3, (long) A, &tmp1);
zmuli(v3, (long) B, &tmp2);
zadd(tmp1, tmp2, &tmp3);
freeh(tmp1.v);
freeh(tmp2.v);
zmuli(u3, (long) C, &tmp1);
zmuli(v3, (long) D, &tmp2);
freeh(u3.v);
freeh(v3.v);
u3 = tmp3;
zadd(tmp1, tmp2, &v3);
freeh(tmp1.v);
freeh(tmp2.v);
}
/*
* Here when the remaining numbers become single precision in size.
* Finish the procedure using single precision calculations.
*/
if (iszero(v3) && !isone(u3)) {
freeh(u3.v);
freeh(v3.v);
freeh(u2.v);
freeh(v2.v);
return TRUE;
}
ui3 = (istiny(u3) ? z1tol(u3) : z2tol(u3));
vi3 = (istiny(v3) ? z1tol(v3) : z2tol(v3));
freeh(u3.v);
freeh(v3.v);
while (vi3) {
q1 = ui3 / vi3;
zmuli(v2, (long) q1, &tmp1);
zsub(u2, tmp1, &tmp2);
freeh(tmp1.v);
freeh(u2.v);
u2 = v2;
v2 = tmp2;
q2 = ui3 - q1 * vi3;
ui3 = vi3;
vi3 = q2;
}
freeh(v2.v);
if (ui3 != 1) {
freeh(u2.v);
return TRUE;
}
if (isneg(u2)) {
zadd(v, u2, res);
freeh(u2.v);
return FALSE;
}
*res = u2;
return FALSE;
}
#if 0
/*
* Approximate the quotient of two integers by another set of smaller
* integers. This uses continued fractions to determine the smaller set.
*/
void
zapprox(z1, z2, res1, res2)
ZVALUE z1, z2, *res1, *res2;
{
int sign;
ZVALUE u1, v1, u3, v3, q, t1, t2, t3;
sign = ((z1.sign != 0) ^ (z2.sign != 0));
z1.sign = 0;
z2.sign = 0;
v3 = z2;
u3 = z1;
u1 = _one_;
v1 = _zero_;
while (!iszero(v3)) {
zdiv(u3, v3, &q, &t1);
zmul(v1, q, &t2);
zsub(u1, t2, &t3);
freeh(q.v);
freeh(t2.v);
freeh(u1.v);
if ((u3.v != z1.v) && (u3.v != z2.v))
freeh(u3.v);
u1 = v1;
u3 = v3;
v1 = t3;
v3 = t1;
}
if (!isunit(u3))
error("Non-relativly prime numbers for approx");
if ((u3.v != z1.v) && (u3.v != z2.v))
freeh(u3.v);
if ((v3.v != z1.v) && (v3.v != z2.v))
freeh(v3.v);
freeh(v1.v);
zmul(u1, z1, &t1);
zsub(t1, _one_, &t2);
freeh(t1.v);
zquo(t2, z2, &t1);
freeh(t2.v);
u1.sign = (BOOL)sign;
t1.sign = 0;
*res1 = t1;
*res2 = u1;
}
#endif
/*
* Binary gcd algorithm
* This algorithm taken from Knuth
*/
void
zgcd(z1, z2, res)
ZVALUE z1, z2, *res;
{
ZVALUE u, v, t;
register long j, k, olen, mask;
register HALF h;
HALF *oldv1, *oldv2;
/*
* First see if one number is very much larger than the other.
* If so, then divide as necessary to get the numbers near each other.
*/
z1.sign = 0;
z2.sign = 0;
oldv1 = z1.v;
oldv2 = z2.v;
if (z1.len < z2.len) {
t = z1;
z1 = z2;
z2 = t;
}
while ((z1.len > (z2.len + 5)) && !iszero(z2)) {
zmod(z1, z2, &t);
if ((z1.v != oldv1) && (z1.v != oldv2))
freeh(z1.v);
z1 = z2;
z2 = t;
}
/*
* Ok, now do the binary method proper
*/
u.len = z1.len;
v.len = z2.len;
u.sign = 0;
v.sign = 0;
if (!ztest(z1)) {
v.v = alloc(v.len);
copyval(z2, v);
*res = v;
goto done;
}
if (!ztest(z2)) {
u.v = alloc(u.len);
copyval(z1, u);
*res = u;
goto done;
}
u.v = alloc(u.len);
v.v = alloc(v.len);
copyval(z1, u);
copyval(z2, v);
k = 0;
while (u.v[k] == 0 && v.v[k] == 0)
++k;
mask = 01;
h = u.v[k] | v.v[k];
k *= BASEB;
while (!(h & mask)) {
mask <<= 1;
k++;
}
shiftr(u, k);
shiftr(v, k);
trim(&u);
trim(&v);
if (isodd(u)) {
t.v = alloc(v.len);
t.len = v.len;
copyval(v, t);
t.sign = !v.sign;
} else {
t.v = alloc(u.len);
t.len = u.len;
copyval(u, t);
t.sign = u.sign;
}
while (ztest(t)) {
j = 0;
while (!t.v[j])
++j;
mask = 01;
h = t.v[j];
j *= BASEB;
while (!(h & mask)) {
mask <<= 1;
j++;
}
shiftr(t, j);
trim(&t);
if (ztest(t) > 0) {
freeh(u.v);
u = t;
} else {
freeh(v.v);
v = t;
v.sign = !t.sign;
}
zsub(u, v, &t);
}
freeh(t.v);
freeh(v.v);
if (k) {
olen = u.len;
u.len += k / BASEB + 1;
u.v = (HALF *) realloc(u.v, u.len * sizeof(HALF));
while (olen != u.len)
u.v[olen++] = 0;
shiftl(u, k);
}
trim(&u);
*res = u;
done:
if ((z1.v != oldv1) && (z1.v != oldv2))
freeh(z1.v);
if ((z2.v != oldv1) && (z2.v != oldv2))
freeh(z2.v);
}
/*
* Compute the lcm of two integers (least common multiple).
* This is done using the formula: gcd(a,b) * lcm(a,b) = a * b.
*/
void
zlcm(z1, z2, res)
ZVALUE z1, z2, *res;
{
ZVALUE temp1, temp2;
zgcd(z1, z2, &temp1);
zquo(z1, temp1, &temp2);
freeh(temp1.v);
zmul(temp2, z2, res);
freeh(temp2.v);
}
/*
* Return whether or not two numbers are relatively prime to each other.
*/
BOOL
zrelprime(z1, z2)
ZVALUE z1, z2; /* numbers to be tested */
{
FULL rem1, rem2; /* remainders */
ZVALUE rem;
BOOL result;
z1.sign = 0;
z2.sign = 0;
if (iseven(z1) && iseven(z2)) /* false if both even */
return FALSE;
if (isunit(z1) || isunit(z2)) /* true if either is a unit */
return TRUE;
if (iszero(z1) || iszero(z2)) /* false if either is zero */
return FALSE;
if (istwo(z1) || istwo(z2)) /* true if either is two */
return TRUE;
/*
* Try reducing each number by the product of the first few odd primes
* to see if any of them are a common factor.
*/
rem1 = zmodi(z1, 3L * 5 * 7 * 11 * 13);
rem2 = zmodi(z2, 3L * 5 * 7 * 11 * 13);
if (((rem1 % 3) == 0) && ((rem2 % 3) == 0))
return FALSE;
if (((rem1 % 5) == 0) && ((rem2 % 5) == 0))
return FALSE;
if (((rem1 % 7) == 0) && ((rem2 % 7) == 0))
return FALSE;
if (((rem1 % 11) == 0) && ((rem2 % 11) == 0))
return FALSE;
if (((rem1 % 13) == 0) && ((rem2 % 13) == 0))
return FALSE;
/*
* Try a new batch of primes now
*/
rem1 = zmodi(z1, 17L * 19 * 23);
rem2 = zmodi(z2, 17L * 19 * 23);
if (((rem1 % 17) == 0) && ((rem2 % 17) == 0))
return FALSE;
if (((rem1 % 19) == 0) && ((rem2 % 19) == 0))
return FALSE;
if (((rem1 % 23) == 0) && ((rem2 % 23) == 0))
return FALSE;
/*
* Yuk, we must actually compute the gcd to know the answer
*/
zgcd(z1, z2, &rem);
result = isunit(rem);
freeh(rem.v);
return result;
}
/*
* Compute the log of one number base another, to the closest integer.
* This is the largest integer which when the second number is raised to it,
* the resulting value is less than or equal to the first number.
* Example: zlog(123456, 10) = 5.
*/
long
zlog(z1, z2)
ZVALUE z1, z2;
{
register ZVALUE *zp; /* current square */
long power; /* current power */
long worth; /* worth of current square */
ZVALUE val; /* current value of power */
ZVALUE temp; /* temporary */
ZVALUE squares[32]; /* table of squares of base */
/*
* Make sure that the numbers are nonzero and the base is greater than one.
*/
if (isneg(z1) || iszero(z1) || isneg(z2) || isleone(z2))
error("Bad arguments for log");
/*
* Reject trivial cases.
*/
if (z1.len < z2.len)
return 0;
if ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))
return 0;
power = zrel(z1, z2);
if (power <= 0)
return (power + 1);
/*
* Handle any power of two special.
*/
if (zisonebit(z2))
return (zhighbit(z1) / zlowbit(z2));
/*
* Handle base 10 special
*/
if ((z2.len == 1) && (*z2.v == 10))
return zlog10(z1);
/*
* Now loop by squaring the base each time, and see whether or
* not each successive square is still smaller than the number.
*/
worth = 1;
zp = &squares[0];
*zp = z2;
while (((zp->len * 2) - 1) <= z1.len) { /* while square not too large */
zsquare(*zp, zp + 1);
zp++;
worth *= 2;
}
/*
* Now back down the squares, and multiply them together to see
* exactly how many times the base can be raised by.
*/
val = _one_;
power = 0;
for (; zp >= squares; zp--, worth /= 2) {
if ((val.len + zp->len - 1) <= z1.len) {
zmul(val, *zp, &temp);
if (zrel(z1, temp) >= 0) {
freeh(val.v);
val = temp;
power += worth;
} else
freeh(temp.v);
}
if (zp != squares)
freeh(zp->v);
}
return power;
}
/*
* Return the integral log base 10 of a number.
*/
long
zlog10(z)
ZVALUE z;
{
register ZVALUE *zp; /* current square */
long power; /* current power */
long worth; /* worth of current square */
ZVALUE val; /* current value of power */
ZVALUE temp; /* temporary */
if (!ispos(z))
error("Non-positive number for log10");
/*
* Loop by squaring the base each time, and see whether or
* not each successive square is still smaller than the number.
*/
worth = 1;
zp = &_tenpowers_[0];
*zp = _ten_;
while (((zp->len * 2) - 1) <= z.len) { /* while square not too large */
if (zp[1].len == 0)
zsquare(*zp, zp + 1);
zp++;
worth *= 2;
}
/*
* Now back down the squares, and multiply them together to see
* exactly how many times the base can be raised by.
*/
val = _one_;
power = 0;
for (; zp >= _tenpowers_; zp--, worth /= 2) {
if ((val.len + zp->len - 1) <= z.len) {
zmul(val, *zp, &temp);
if (zrel(z, temp) >= 0) {
freeh(val.v);
val = temp;
power += worth;
} else
freeh(temp.v);
}
}
return power;
}
/*
* Return the number of times that one number will divide another.
* This works similarly to zlog, except that divisions must be exact.
* For example, zdivcount(540, 3) = 3, since 3^3 divides 540 but 3^4 won't.
*/
long
zdivcount(z1, z2)
ZVALUE z1, z2;
{
long count; /* number of factors removed */
ZVALUE tmp; /* ignored return value */
count = zfacrem(z1, z2, &tmp);
freeh(tmp.v);
return count;
}
/*
* Remove all occurances of the specified factor from a number.
* Also returns the number of factors removed as a function return value.
* Example: zfacrem(540, 3, &x) returns 3 and sets x to 20.
*/
long
zfacrem(z1, z2, rem)
ZVALUE z1, z2, *rem;
{
register ZVALUE *zp; /* current square */
long count; /* total count of divisions */
long worth; /* worth of current square */
ZVALUE temp1, temp2, temp3; /* temporaries */
ZVALUE squares[32]; /* table of squares of factor */
z1.sign = 0;
z2.sign = 0;
/*
* Make sure factor isn't 0 or 1.
*/
if (isleone(z2))
error("Bad argument for facrem");
/*
* Reject trivial cases.
*/
if ((z1.len < z2.len) || (isodd(z1) && iseven(z2)) ||
((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))) {
rem->v = alloc(z1.len);
rem->len = z1.len;
rem->sign = 0;
copyval(z1, *rem);
return 0;
}
/*
* Handle any power of two special.
*/
if (zisonebit(z2)) {
count = zlowbit(z1) / zlowbit(z2);
rem->v = alloc(z1.len);
rem->len = z1.len;
rem->sign = 0;
copyval(z1, *rem);
shiftr(*rem, count);
trim(rem);
return count;
}
/*
* See if the factor goes in even once.
*/
zdiv(z1, z2, &temp1, &temp2);
if (!iszero(temp2)) {
freeh(temp1.v);
freeh(temp2.v);
rem->v = alloc(z1.len);
rem->len = z1.len;
rem->sign = 0;
copyval(z1, *rem);
return 0;
}
freeh(temp2.v);
z1 = temp1;
/*
* Now loop by squaring the factor each time, and see whether
* or not each successive square will still divide the number.
*/
count = 1;
worth = 1;
zp = &squares[0];
*zp = z2;
while (((zp->len * 2) - 1) <= z1.len) { /* while square not too large */
zsquare(*zp, &temp1);
zdiv(z1, temp1, &temp2, &temp3);
if (!iszero(temp3)) {
freeh(temp1.v);
freeh(temp2.v);
freeh(temp3.v);
break;
}
freeh(temp3.v);
freeh(z1.v);
z1 = temp2;
*++zp = temp1;
worth *= 2;
count += worth;
}
/*
* Now back down the list of squares, and see if the lower powers
* will divide any more times.
*/
for (; zp >= squares; zp--, worth /= 2) {
if (zp->len <= z1.len) {
zdiv(z1, *zp, &temp1, &temp2);
if (iszero(temp2)) {
temp3 = z1;
z1 = temp1;
temp1 = temp3;
count += worth;
}
freeh(temp1.v);
freeh(temp2.v);
}
if (zp != squares)
freeh(zp->v);
}
*rem = z1;
return count;
}
/*
* Keep dividing a number by the gcd of it with another number until the
* result is relatively prime to the second number.
*/
void
zgcdrem(z1, z2, res)
ZVALUE z1, z2, *res;
{
ZVALUE tmp1, tmp2;
/*
* Begin by taking the gcd for the first time.
* If the number is already relatively prime, then we are done.
*/
zgcd(z1, z2, &tmp1);
if (isunit(tmp1) || iszero(tmp1)) {
res->len = z1.len;
res->v = alloc(z1.len);
res->sign = z1.sign;
copyval(z1, *res);
return;
}
zquo(z1, tmp1, &tmp2);
z1 = tmp2;
z2 = tmp1;
/*
* Now keep alternately taking the gcd and removing factors until
* the gcd becomes one.
*/
while (!isunit(z2)) {
(void) zfacrem(z1, z2, &tmp1);
freeh(z1.v);
z1 = tmp1;
zgcd(z1, z2, &tmp1);
freeh(z2.v);
z2 = tmp1;
}
*res = z1;
}
/*
* Find the lowest prime factor of a number if one can be found.
* Search is conducted for the first count primes.
* Returns 1 if no factor was found.
*/
long
zlowfactor(z, count)
ZVALUE z;
long count;
{
FULL p, d;
ZVALUE div, tmp;
HALF divval[2];
if ((--count < 0) || iszero(z))
return 1;
if (iseven(z))
return 2;
div.sign = 0;
div.v = divval;
for (p = 3; (count > 0); p += 2) {
for (d = 3; (d * d) <= p; d += 2)
if ((p % d) == 0)
goto next;
divval[0] = (p & BASE1);
divval[1] = (p / BASE);
div.len = 1 + (p >= BASE);
zmod(z, div, &tmp);
if (iszero(tmp))
return p;
freeh(tmp.v);
count--;
next:;
}
return 1;
}
/*
* Return the number of digits (base 10) in a number, ignoring the sign.
*/
long
zdigits(z1)
ZVALUE z1;
{
long count, val;
z1.sign = 0;
if (istiny(z1)) { /* do small numbers ourself */
count = 1;
val = 10;
while (*z1.v >= (HALF)val) {
count++;
val *= 10;
}
return count;
}
return (zlog10(z1) + 1);
}
/*
* Return the single digit at the specified decimal place of a number,
* where 0 means the rightmost digit. Example: zdigit(1234, 1) = 3.
*/
FLAG
zdigit(z1, n)
ZVALUE z1;
long n;
{
ZVALUE tmp1, tmp2;
FLAG res;
z1.sign = 0;
if (iszero(z1) || (n < 0) || (n / BASEDIG >= z1.len))
return 0;
if (n == 0)
return zmodi(z1, 10L);
if (n == 1)
return zmodi(z1, 100L) / 10;
if (n == 2)
return zmodi(z1, 1000L) / 100;
if (n == 3)
return zmodi(z1, 10000L) / 1000;
ztenpow(n, &tmp1);
zquo(z1, tmp1, &tmp2);
res = zmodi(tmp2, 10L);
freeh(tmp1.v);
freeh(tmp2.v);
return res;
}
/*
* Find the square root of a number. This is the greatest integer whose
* square is less than or equal to the number. Returns TRUE if the
* square root is exact.
*/
BOOL
zsqrt(z1, dest)
ZVALUE z1, *dest;
{
ZVALUE try, quo, rem, old, temp;
FULL iquo, val;
long i,j;
if (z1.sign)
error("Square root of negative number");
if (iszero(z1)) {
*dest = _zero_;
return TRUE;
}
if ((*z1.v < 4) && istiny(z1)) {
*dest = _one_;
return (*z1.v == 1);
}
/*
* Pick the square root of the leading one or two digits as a first guess.
*/
val = z1.v[z1.len-1];
if ((z1.len & 0x1) == 0)
val = (val * BASE) + z1.v[z1.len-2];
/*
* Find the largest power of 2 that when squared
* is <= val > 0. Avoid multiply overflow by doing
* a careful check at the BASE boundary.
*/
j = 1<<(BASEB+BASEB-2);
if (val > j) {
iquo = BASE;
} else {
i = BASEB-1;
while (j > val) {
--i;
j >>= 2;
}
iquo = bitmask[i];
}
for (i = 8; i > 0; i--)
iquo = (iquo + (val / iquo)) / 2;
if (iquo > BASE1)
iquo = BASE1;
/*
* Allocate the numbers to use for the main loop.
* The size and high bits of the final result are correctly set here.
* Notice that the remainder of the test value is rubbish, but this
* is unimportant.
*/
try.sign = 0;
try.len = (z1.len + 1) / 2;
try.v = alloc(try.len);
try.v[try.len-1] = (HALF)iquo;
old.sign = 0;
old.v = alloc(try.len);
old.len = 1;
*old.v = 0;
/*
* Main divide and average loop
*/
for (;;) {
zdiv(z1, try, &quo, &rem);
i = zrel(try, quo);
if ((i == 0) && iszero(rem)) { /* exact square root */
freeh(rem.v);
freeh(quo.v);
freeh(old.v);
*dest = try;
return TRUE;
}
freeh(rem.v);
if (i <= 0) {
/*
* Current try is less than or equal to the square root since
* it is less than the quotient. If the quotient is equal to
* the try, we are all done. Also, if the try is equal to the
* old try value, we are done since no improvement occurred.
* If not, save the improved value and loop some more.
*/
if ((i == 0) || (zcmp(old, try) == 0)) {
freeh(quo.v);
freeh(old.v);
*dest = try;
return FALSE;
}
old.len = try.len;
copyval(try, old);
}
/* average current try and quotent for the new try */
zadd(try, quo, &temp);
freeh(quo.v);
freeh(try.v);
try = temp;
shiftr(try, 1L);
quicktrim(try);
}
}
/*
* Take an arbitrary root of a number (to the greatest integer).
* This uses the following iteration to get the Kth root of N:
* x = ((K-1) * x + N / x^(K-1)) / K
*/
void
zroot(z1, z2, dest)
ZVALUE z1, z2, *dest;
{
ZVALUE try, quo, old, temp, temp2;
ZVALUE k1; /* holds k - 1 */
int sign;
long i, k, highbit;
SIUNION sival;
sign = z1.sign;
if (sign && iseven(z2))
error("Even root of negative number");
if (iszero(z2) || isneg(z2))
error("Non-positive root");
if (iszero(z1)) { /* root of zero */
*dest = _zero_;
return;
}
if (isunit(z2)) { /* first root */
zcopy(z1, dest);
return;
}
if (isbig(z2)) { /* humongous root */
*dest = _one_;
dest->sign = (HALF)sign;
return;
}
k = (istiny(z2) ? z1tol(z2) : z2tol(z2));
highbit = zhighbit(z1);
if (highbit < k) { /* too high a root */
*dest = _one_;
dest->sign = (HALF)sign;
return;
}
sival.ivalue = k - 1;
k1.v = &sival.silow;
k1.len = 1 + (sival.sihigh != 0);
k1.sign = 0;
z1.sign = 0;
/*
* Allocate the numbers to use for the main loop.
* The size and high bits of the final result are correctly set here.
* Notice that the remainder of the test value is rubbish, but this
* is unimportant.
*/
highbit = (highbit + k - 1) / k;
try.len = (highbit / BASEB) + 1;
try.v = alloc(try.len);
try.v[try.len-1] = ((HALF)1 << (highbit % BASEB));
try.sign = 0;
old.v = alloc(try.len);
old.len = 1;
*old.v = 0;
old.sign = 0;
/*
* Main divide and average loop
*/
for (;;) {
zpowi(try, k1, &temp);
zquo(z1, temp, &quo);
freeh(temp.v);
i = zrel(try, quo);
if (i <= 0) {
/*
* Current try is less than or equal to the root since it is
* less than the quotient. If the quotient is equal to the try,
* we are all done. Also, if the try is equal to the old value,
* we are done since no improvement occurred.
* If not, save the improved value and loop some more.
*/
if ((i == 0) || (zcmp(old, try) == 0)) {
freeh(quo.v);
freeh(old.v);
try.sign = (HALF)sign;
quicktrim(try);
*dest = try;
return;
}
old.len = try.len;
copyval(try, old);
}
/* average current try and quotent for the new try */
zmul(try, k1, &temp);
freeh(try.v);
zadd(quo, temp, &temp2);
freeh(temp.v);
freeh(quo.v);
zquo(temp2, z2, &try);
freeh(temp2.v);
}
}
/*
* Test to see if a number is an exact square or not.
*/
BOOL
zissquare(z)
ZVALUE z; /* number to be tested */
{
long n, i;
ZVALUE tmp;
if (z.sign) /* negative */
return FALSE;
while ((z.len > 1) && (*z.v == 0)) { /* ignore trailing zero words */
z.len--;
z.v++;
}
if (isleone(z)) /* zero or one */
return TRUE;
n = *z.v & 0xf; /* check mod 16 values */
if ((n != 0) && (n != 1) && (n != 4) && (n != 9))
return FALSE;
n = *z.v & 0xff; /* check mod 256 values */
i = 0x80;
while (((i * i) & 0xff) != n)
if (--i <= 0)
return FALSE;
n = zsqrt(z, &tmp); /* must do full square root test now */
freeh(tmp.v);
return n;
}
/* END CODE */