zfunc.c

早期freebsd实现

 */
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 (!zispos(z))
		math_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) {
				zfree(val);
				val = temp;
				power += worth;
			} else
				zfree(temp);
		}
	}
	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);
	zfree(tmp);
	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 (zisleone(z2))
		math_error("Bad argument for facrem");
	/*
	 * Reject trivial cases.
	 */
	if ((z1.len < z2.len) || (zisodd(z1) && ziseven(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;
		zcopyval(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;
		zcopyval(z1, *rem);
		zshiftr(*rem, count);
		ztrim(rem);
		return count;
	}
	/*
	 * See if the factor goes in even once.
	 */
	zdiv(z1, z2, &temp1, &temp2);
	if (!ziszero(temp2)) {
		zfree(temp1);
		zfree(temp2);
		rem->v = alloc(z1.len);
		rem->len = z1.len;
		rem->sign = 0;
		zcopyval(z1, *rem);
		return 0;
	}
	zfree(temp2);
	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 (!ziszero(temp3)) {
			zfree(temp1);
			zfree(temp2);
			zfree(temp3);
			break;
		}
		zfree(temp3);
		zfree(z1);
		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 (ziszero(temp2)) {
				temp3 = z1;
				z1 = temp1;
				temp1 = temp3;
				count += worth;
			}
			zfree(temp1);
			zfree(temp2);
		}
		if (zp != squares)
			zfree(*zp);
	}
	*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 (zisunit(tmp1) || ziszero(tmp1)) {
		res->len = z1.len;
		res->v = alloc(z1.len);
		res->sign = z1.sign;
		zcopyval(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 (!zisunit(z2)) {
		(void) zfacrem(z1, z2, &tmp1);
		zfree(z1);
		z1 = tmp1;
		zgcd(z1, z2, &tmp1);
		zfree(z2);
		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) || ziszero(z))
		return 1;
	if (ziseven(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 (ziszero(tmp))
			return p;
		zfree(tmp);
		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 (zistiny(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 (ziszero(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);
	zfree(tmp1);
	zfree(tmp2);
	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)
		math_error("Square root of negative number");
	if (ziszero(z1)) {
		*dest = _zero_;
		return TRUE;
	}
	if ((*z1.v < 4) && zistiny(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 = 1L<<(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);
	zclearval(try);
	try.v[try.len-1] = (HALF)iquo;
	old.sign = 0;
	old.v = alloc(try.len);
	old.len = 1;
	zclearval(old);
	/*
	 * Main divide and average loop
	 */
	for (;;) {
		zdiv(z1, try, &quo, &rem);
		i = zrel(try, quo);
		if ((i == 0) && ziszero(rem)) {	/* exact square root */
			zfree(rem);
			zfree(quo);
			zfree(old);
			*dest = try;
			return TRUE;
		}
		zfree(rem);
		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)) {
				zfree(quo);
				zfree(old);
				*dest = try;
				return FALSE;
			}
			old.len = try.len;
			zcopyval(try, old);
		}
		/* average current try and quotent for the new try */
		zadd(try, quo, &temp);
		zfree(quo);
		zfree(try);
		try = temp;
		zshiftr(try, 1L);
		zquicktrim(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 && ziseven(z2))
		math_error("Even root of negative number");
	if (ziszero(z2) || zisneg(z2))
		math_error("Non-positive root");
	if (ziszero(z1)) {	/* root of zero */
		*dest = _zero_;
		return;
	}
	if (zisunit(z2)) {	/* first root */
		zcopy(z1, dest);
		return;
	}
	if (zisbig(z2)) {	/* humongous root */
		*dest = _one_;
		dest->sign = (HALF)sign;
		return;
	}
	k = (zistiny(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);
	zclearval(try);
	try.v[try.len-1] = ((HALF)1 << (highbit % BASEB));
	try.sign = 0;
	old.v = alloc(try.len);
	old.len = 1;
	zclearval(old);
	old.sign = 0;
	/*
	 * Main divide and average loop
	 */
	for (;;) {
		zpowi(try, k1, &temp);
		zquo(z1, temp, &quo);
		zfree(temp);
		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)) {
				zfree(quo);
				zfree(old);
				try.sign = (HALF)sign;
				zquicktrim(try);
				*dest = try;
				return;
			}
			old.len = try.len;
			zcopyval(try, old);
		}
		/* average current try and quotent for the new try */
		zmul(try, k1, &temp);
		zfree(try);
		zadd(quo, temp, &temp2);
		zfree(temp);
		zfree(quo);
		zquo(temp2, z2, &try);
		zfree(temp2);
	}
}


/*
 * 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 (zisleone(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 */
	zfree(tmp);
	return n;
}


/*
 * Return a trivial hash value for an integer.
 */
HASH
zhash(z)
	ZVALUE z;
{
	HASH hash;
	int i;

	hash = z.len * 1000003;
	if (z.sign)
		hash += 10000019;
	for (i = z.len - 1; i >= 0; i--)
		/* ignore Saber-C warning about Over/underflow */
		hash = hash * 79372817 + z.v[i] + 10000079;
	return hash;
}

/* END CODE */