zfunc.c

早期freebsd实现

/*
 * Copyright (c) 1994 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 "zmath.h"

static ZVALUE primeprod;		/* product of primes under 100 */
ZVALUE _tenpowers_[2 * BASEB];		/* table of 10^2^n */


/*
 * 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 (zisneg(z))
		math_error("Negative argument for factorial");
	if (zisbig(z))
		math_error("Very large factorial");
	n = (zistiny(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);
		zfree(res);
		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);
		zfree(res);
		res = temp;
	}
	zshift(res, ptwo, &temp);
	zfree(res);
	*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 (zisneg(z))
		math_error("Negative argument for factorial");
	if (zisbig(z))
		math_error("Very large factorial");
	n = (zistiny(z) ? z1tol(z) : z2tol(z));
	/*
	 * Multiply by the primes in order, collecting multiple numbers
	 * together until there is a chance 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);
		zfree(res);
		res = temp;
		mul = 1;
next: ;
	}
	/*
	 * Multiply by the final value if any.
	 */
	if (mul > 1) {
		zmuli(res, mul, &temp);
		zfree(res);
		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 = 0;		/* power of prime */
	long i;			/* test value */
	ZVALUE res, temp;

	if (zisneg(z) || ziszero(z))
		math_error("Non-positive argument for lcmfact");
	if (zisbig(z))
		math_error("Very large lcmfact");
	n = (zistiny(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);
		zfree(res);
		res = temp;
next: ;
	}
	/*
	 * Finish by scaling by the necessary power of two.
	 */
	zshift(res, zhighbit(z), dest);
	zfree(res);
}


/*
 * Compute the permutation function  M! / (M - N)!.
 */
void
zperm(z1, z2, res)
	ZVALUE z1, z2, *res;
{
	long count;
	ZVALUE cur, tmp, ans;

	if (zisneg(z1) || zisneg(z2))
		math_error("Negative argument for permutation");
	if (zrel(z1, z2) < 0)
		math_error("Second arg larger than first in permutation");
	if (zisbig(z2))
		math_error("Very large permutation");
	count = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
	zcopy(z1, &ans);
	zsub(z1, _one_, &cur);
	while (--count > 0) {
		zmul(ans, cur, &tmp);
		zfree(ans);
		ans = tmp;
		zsub(cur, _one_, &tmp);
		zfree(cur);
		cur = tmp;
	}
	zfree(cur);
	*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 (zisneg(z1) || zisneg(z2))
		math_error("Negative argument for combinatorial");
	zsub(z1, z2, &temp);
	if (zisneg(temp)) {
		zfree(temp);
		math_error("Second arg larger than first for combinatorial");
	}
	if (zisbig(z2) && zisbig(temp)) {
		zfree(temp);
		math_error("Very large combinatorial");
	}
	count = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
	i = (zistiny(temp) ? z1tol(temp) : z2tol(temp));
	if (zisbig(z2) || (!zisbig(temp) && (i < count)))
		count = i;
	zfree(temp);
	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);
		zfree(ans);
		zquo(temp, div, &ans);
		zfree(temp);
		zsub(mul, _one_, &temp);
		if (mul.v != z1.v)
			zfree(mul);
		mul = temp;
	}
	if (mul.v != z1.v)
		zfree(mul);
	*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 (ziseven(z))		/* if even, not prime if not 2 */
		return (zistwo(z) != 0);
	/*
	 * See if the number is small, and is either a small prime,
	 * or is divisible by a small prime.
	 */
	if (zistiny(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 (!zisunit(z1)) {
		zfree(z1);
		return FALSE;
	}
	zfree(z1);
	/*
	 * 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 = ((zistiny(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 (zisone(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);
			zfree(z3);
			zmod(z2, z, &z3);
			zfree(z2);
		}
		zfree(z3);
	}
	zfree(zm1);
	zfree(z1);
	return TRUE;	/* number might be prime */

notprime:
	zfree(z3);
	zfree(zm1);
	zfree(z1);
	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 (ziseven(z2) || zisneg(z2))
		return 0;
	val = 1;
	if (ziszero(z1) || zisone(z1))
		return val;
	if (zisunit(z1)) {
		if ((*z2.v - 1) & 0x2)
			val = -val;
		return val;
	}
	zcopy(z1, &p);
	zcopy(z2, &q);
	for (;;) {
		zmod(p, q, &tmp);
		zfree(p);
		p = tmp;
		if (ziszero(p)) {
			zfree(p);
			p = _one_;
		}
		if (ziseven(p)) {
			lowbit = zlowbit(p);
			zshift(p, -lowbit, &tmp);
			zfree(p);
			p = tmp;
			if ((lowbit & 1) && (((*q.v & 0x7) == 3) || ((*q.v & 0x7) == 5)))
				val = -val;
		}
		if (zisunit(p)) {
			zfree(p);
			zfree(q);
			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 (zisbig(z))
		math_error("Very large Fibonacci number");
	n = (zistiny(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);
		zfree(fnm1);
		zfree(fn);
		zfree(fnp1);
		zadd(t2, t3, &fnp1);
		zsub(t3, t1, &fn);
		zfree(t1);
		zfree(t2);
		zfree(t3);
		if (i & n) {
			fnm1 = fn;
			fn = fnp1;
			zadd(fnm1, fn, &fnp1);
		} else
			zsub(fnp1, fn, &fnm1);
		i >>= 1L;
	}
	zfree(fnm1);
	zfree(fnp1);
	*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 && zisodd(z2));
	z1.sign = 0;
	z2.sign = 0;
	if (ziszero(z2) && !ziszero(z1)) {	/* number raised to power 0 */
		*res = _one_;
		return;
	}
	if (zisleone(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 (zisbig(z2))
		math_error("Raising to very large power");
	power = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
	if (zistwo(z1)) {	/* two raised to a power */
		zbitvalue((long) power, res);
		return;
	}
	/*
	 * See if this is a power of ten
	 */
	if (zistiny(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);
				zcopyval(z1, ans);
				ans.sign = (BOOL)sign;
				*res = ans;
				return;
			case 2:
				zsquare(z1, res);
				return;
			case 3:
				zsquare(z1, &temp);
				zmul(z1, temp, res);
				zfree(temp);
				res->sign = (BOOL)sign;
				return;
			case 4:
				zsquare(z1, &temp);
				zsquare(temp, res);
				zfree(temp);
				return;
		}
	}
	/*
	 * Shift out all powers of twos so the multiplies are smaller.
	 * We will shift back the right amount when done.
	 */
	twos = 0;
	if (ziseven(z1)) {
		twos = zlowbit(z1);
		ans.v = alloc(z1.len);
		ans.len = z1.len;
		ans.sign = z1.sign;
		zcopyval(z1, ans);
		zshiftr(ans, twos);
		ztrim(&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) {