zfunc.c

早期freebsd实现

		zmul(ans, z1, &temp);
		zfree(ans);
		ans = temp;
	}
	bit >>= 1L;
	while (bit) {
		zsquare(ans, &temp);
		zfree(ans);
		ans = temp;
		if (bit & power) {
			zmul(ans, z1, &temp);
			zfree(ans);
			ans = temp;
		}
		bit >>= 1L;
	}
	/*
	 * Scale back up by proper power of two
	 */
	if (twos) {
		zshift(ans, twos, &temp);
		zfree(ans);
		ans = temp;
		zfree(z1);
	}
	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);
			zfree(ans);
			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 (zisneg(u) || zisneg(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) && !ziszero(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);
			zfree(tmp1);
			zfree(u2);
			u2 = v2;
			v2 = tmp2;
			zmul(qz, v3, &tmp1);
			zsub(u3, tmp1, &tmp2);
			zfree(tmp1);
			zfree(u3);
			u3 = v3;
			v3 = tmp2;
			zfree(qz);
			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);
		zfree(tmp1);
		zfree(tmp2);
		zmuli(u2, (long) C, &tmp1);
		zmuli(v2, (long) D, &tmp2);
		zfree(u2);
		zfree(v2);
		u2 = tmp3;
		zadd(tmp1, tmp2, &v2);
		zfree(tmp1);
		zfree(tmp2);
		zmuli(u3, (long) A, &tmp1);
		zmuli(v3, (long) B, &tmp2);
		zadd(tmp1, tmp2, &tmp3);
		zfree(tmp1);
		zfree(tmp2);
		zmuli(u3, (long) C, &tmp1);
		zmuli(v3, (long) D, &tmp2);
		zfree(u3);
		zfree(v3);
		u3 = tmp3;
		zadd(tmp1, tmp2, &v3);
		zfree(tmp1);
		zfree(tmp2);
	}

	/*
	 * Here when the remaining numbers become single precision in size.
	 * Finish the procedure using single precision calculations.
	 */
	if (ziszero(v3) && !zisone(u3)) {
		zfree(u3);
		zfree(v3);
		zfree(u2);
		zfree(v2);
		return TRUE;
	}
	ui3 = (zistiny(u3) ? z1tol(u3) : z2tol(u3));
	vi3 = (zistiny(v3) ? z1tol(v3) : z2tol(v3));
	zfree(u3);
	zfree(v3);
	while (vi3) {
		q1 = ui3 / vi3;
		zmuli(v2, (long) q1, &tmp1);
		zsub(u2, tmp1, &tmp2);
		zfree(tmp1);
		zfree(u2);
		u2 = v2;
		v2 = tmp2;
		q2 = ui3 - q1 * vi3;
		ui3 = vi3;
		vi3 = q2;
	}
	zfree(v2);
	if (ui3 != 1) {
		zfree(u2);
		return TRUE;
	}
	if (zisneg(u2)) {
		zadd(v, u2, res);
		zfree(u2);
		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 (!ziszero(v3)) {
		zdiv(u3, v3, &q, &t1);
		zmul(v1, q, &t2);
		zsub(u1, t2, &t3);
		zfree(q);
		zfree(t2);
		zfree(u1);
		if ((u3.v != z1.v) && (u3.v != z2.v))
			zfree(u3);
		u1 = v1;
		u3 = v3;
		v1 = t3;
		v3 = t1;
	}
	if (!zisunit(u3))
		math_error("Non-relativly prime numbers for approx");
	if ((u3.v != z1.v) && (u3.v != z2.v))
		zfree(u3);
	if ((v3.v != z1.v) && (v3.v != z2.v))
		zfree(v3);
	zfree(v1);
	zmul(u1, z1, &t1);
	zsub(t1, _one_, &t2);
	zfree(t1);
	zquo(t2, z2, &t1);
	zfree(t2);
	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;

	z1.sign = 0;
	z2.sign = 0;
	if (ziszero(z1)) {
		zcopy(z2, res);
		return;
	}
	if (ziszero(z2)) {
		zcopy(z1, res);
		return;
	}
	if (zisunit(z1) || zisunit(z2)) {
		*res = _one_;
		return;
	}
	/*
	 * 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.
	 */
	oldv1 = z1.v;
	oldv2 = z2.v;
	if (z1.len < z2.len) {
		t = z1;
		z1 = z2;
		z2 = t;
	}
	while ((z1.len > (z2.len + 5)) && !ziszero(z2)) {
		zmod(z1, z2, &t);
		if ((z1.v != oldv1) && (z1.v != oldv2))
			zfree(z1);
		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);
		zcopyval(z2, v);
		*res = v;
		goto done;
	}
	if (!ztest(z2)) {
		u.v = alloc(u.len);
		zcopyval(z1, u);
		*res = u;
		goto done;
	}
	u.v = alloc(u.len);
	v.v = alloc(v.len);
	zcopyval(z1, u);
	zcopyval(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++;
	}
	zshiftr(u, k);
	zshiftr(v, k);
	ztrim(&u);
	ztrim(&v);
	if (zisodd(u)) {
		t.v = alloc(v.len);
		t.len = v.len;
		zcopyval(v, t);
		t.sign = !v.sign;
	} else {
		t.v = alloc(u.len);
		t.len = u.len;
		zcopyval(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++;
		}
		zshiftr(t, j);
		ztrim(&t);
		if (ztest(t) > 0) {
			zfree(u);
			u = t;
		} else {
			zfree(v);
			v = t;
			v.sign = !t.sign;
		}
		zsub(u, v, &t);
	}
	zfree(t);
	zfree(v);
	if (k) {
		olen = u.len;
		u.len += k / BASEB + 1;
		u.v = (HALF *) realloc(u.v, u.len * sizeof(HALF));
		if (u.v == NULL) {
		    math_error("Not enough memory to expand number");
		}
		while (olen != u.len)
			u.v[olen++] = 0;
		zshiftl(u, k);
	}
	ztrim(&u);
	*res = u;

done:
	if ((z1.v != oldv1) && (z1.v != oldv2))
		zfree(z1);
	if ((z2.v != oldv1) && (z2.v != oldv2))
		zfree(z2);
}


/*
 * 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);
	zfree(temp1);
	zmul(temp2, z2, res);
	zfree(temp2);
}


/*
 * 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 (ziseven(z1) && ziseven(z2))	/* false if both even */
		return FALSE;
	if (zisunit(z1) || zisunit(z2))	/* true if either is a unit */
		return TRUE;
	if (ziszero(z1) || ziszero(z2))	/* false if either is zero */
		return FALSE;
	if (zistwo(z1) || zistwo(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 = zisunit(rem);
	zfree(rem);
	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 (zisneg(z1) || ziszero(z1) || zisneg(z2) || zisleone(z2))
		math_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) {
				zfree(val);
				val = temp;
				power += worth;
			} else
				zfree(temp);
		}
		if (zp != squares)
			zfree(*zp);
	}
	return power;
}


/*
 * Return the integral log base 10 of a number.