zfunc.c

Calc Software Package for Number Calc

	}
	power = ztoulong(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 >>= 1;
	bit >>= 1;
	zsquare(z1, &ans);
	if (bit & power) {
		zmul(ans, z1, &temp);
		zfree(ans);
		ans = temp;
	}
	bit >>= 1;
	while (bit) {
		zsquare(ans, &temp);
		zfree(ans);
		ans = temp;
		if (bit & power) {
			zmul(ans, z1, &temp);
			zfree(ans);
			ans = temp;
		}
		bit >>= 1;
	}
	/*
	 * 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(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) {
			if (i <= TEN_MAX) {
				zsquare(_tenpowers_[i-1], &_tenpowers_[i]);
			} else {
				math_error("cannot compute 10^2^(TEN_MAX+1)");
				/*NOTREACHED*/
			}
		}
		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 Euclidean 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(ZVALUE u, ZVALUE v, ZVALUE *res)
{
	FULL	q1, q2, ui3, vi3, uh, vh, A, B, C, D, T;
	ZVALUE	u2, u3, v2, v3, qz, tmp1, tmp2, tmp3;

	v.sign = 0;
	if (zisneg(u) || (zrel(u, v) >= 0))
		zmod(u, v, &v3, 0);
	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 HALF.	Throughout this loop u3 >= v3.
	 */
	while ((u3.len > 1) && !ziszero(v3)) {
		vh = 0;
#if LONG_BITS == BASEB
		uh =  u3.v[u3.len - 1];
		if (v3.len == u3.len)
			vh = v3.v[v3.len - 1];
#else
		uh = (((FULL) u3.v[u3.len - 1]) << BASEB) + u3.v[u3.len - 2];
		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];
#endif
		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, 0);
			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 = ztofull(u3);
	vi3 = ztofull(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;
}


/*
 * Compute the greatest common divisor of a pair of integers.
 */
void
zgcd(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	int h, i, j, k;
	LEN len, l, m, n, o, p, q;
	HALF u, v, w, x;
	HALF *a, *a0, *A, *b, *b0, *B, *c, *d;
	FULL f, g;
	ZVALUE gcd;
	BOOL needw;

	if (zisunit(z1) || zisunit(z2)) {
		*res = _one_;
		return;
	}
	z1.sign = 0;
	z2.sign = 0;
	if (ziszero(z1) || !zcmp(z1, z2)) {
		zcopy(z2, res);
		return;
	}
	if (ziszero(z2)) {
		zcopy(z1, res);
		return;
	}

	o = 0;
	while (!(z1.v[o] | z2.v[o])) o++;	/* Count common zero digits */

	c = z1.v + o;
	d = z2.v + o;

	m = z1.len - o;
	n = z2.len - o;
	u = *c | *d;			/* Count common zero bits */
	v = 1;
	p = 0;
	while (!(u & v)) {
		v <<= 1;
		p++;
	}

	while (!*c) {				/* Removing zero digits */
		c++;
		m--;
	}

	while (!*d) {
		d++;
		n--;
	}


	u = *d;			/* Count zero bits for *d */
	v = 1;
	q = 0;
	while (!(u & v)) {
		v <<= 1;
		q++;
	}

	a0 = A = alloc(m);
	b0 = B = alloc(n);

	memcpy(A, c, m * sizeof(HALF));			/* Copy c[] to A[] */

				/* Copy d[] to B[], shifting if necessary */
	if (q) {
		i = n;
		b = B + n;
		d += n;
		f = 0;
		while (i--) {
			f = f << BASEB | *--d;
			*--b = (HALF) (f >> q);
		}
		if (B[n-1] == 0) n--;
	}
	else memcpy(B, d, n * sizeof(HALF));

	if (n == 1) {		/* One digit case; use Euclid's algorithm */
		n = m;
		b0 = A;
		m = 1;
		a0 = B;
		if (m == 1) {				/* a has one digit */
			v = *a0;
			if (v > 1) {			/* Euclid's algorithm */
				b = b0 + n;
				i = n;
				u = 0;
				while (i--) {
					f = (FULL) u << BASEB | *--b;
					u = (HALF) (f % v);
				}
				while (u) { w = v % u; v = u; u = w; }
			}
			*b0 = v;
			n = 1;
		}
		len = n + o;
		gcd.v = alloc(len + 1);
		/* Common zero digits */
		if (o) memset(gcd.v, 0, o * sizeof(HALF));
		/* Left shift for common zero bits */
		if (p) {
			i = n;
			f = 0;
			b = b0;
			a = gcd.v + o;
			while (i--) {
				f = f >> BASEB | (FULL) *b++ << p;
				*a++ = (HALF) f;
			}
			if (f >>= BASEB) {len++; *a = (HALF) f;}
		} else {
			memcpy(gcd.v + o, b0, n * sizeof(HALF));
		}
		gcd.len = len;
		gcd.sign = 0;
		freeh(A);
		freeh(B);
		*res = gcd;
		return;
	}

	u = B[n-1];				/* Bit count for b */
	k = (n	- 1) * BASEB;
	while (u >>= 1) k++;

	needw = TRUE;

	w = 0;
	j = 0;
	while (m) {				/* START OF MAIN LOOP */
	   if (m - n < 2 || needw) {
		q = 0;
		u = *a0;
		v = 1;
		while (!(u & v)) {		/* count zero bits for *a0 */
			q++;
			v <<= 1;
		}

		if (q) {			/* right-justify a */
			a = a0 + m;
			i = m;
			f = 0;
			while (i--) {
				f = f << BASEB | *--a;
				*a = (HALF) (f >> q);
			}
			if (!a0[m-1]) m--;	/* top digit vanishes */
		}

		if (m == 1) break;

		u = a0[m-1];
		j = (m - 1) * BASEB;
		while (u >>= 1) j++;		/* counting bits for a */
		h = j - k;
		if (h < 0) {			/* swapping to get h > 0 */
			l = m;
			m = n;
			n = l;
			a = a0;
			a0 = b0;
			b0 = a;
			k = j;
			h = -h;
			needw = TRUE;
		}
		if (h > 1) {
			if (needw) {		/* find w = minv(*b0, h0) */
				u = 1;
				v = *b0;
				w = 0;
				x = 1;
				i = h;
				while (i-- && x) {
					if (u & x) { u -= v * x; w |= x;}
					x <<= 1;
				}
				needw = FALSE;
			}
			g = *a0 * w;
			if (h < BASEB) g &= (1 << h) - 1;
			else g &= BASE1;
		}
		else g = 1;
	   } else
		g = (HALF) *a0 * w;
		a = a0;
		b = b0;
		i = n;
		if (g > 1) {			/* a - g * b case */
			f = 0;
			while (i--) {
				f = (FULL) *a - g * *b++ - f;
				*a++ = (HALF) f;
				f >>= BASEB;
				f = -f & BASE1;
			}
			if (f) {
				i = m - n;
				while (i-- && f) {
					f = *a - f;
					*a++ = (HALF) f;
					f >>= BASEB;
					f = -f & BASE1;
				}
			}
			while (m && !*a0) { /* Removing trailing zeros */
				m--;
				a0++;
			}
			if (f) {		/* a - g * b < 0 */
				while (m > 1 && a0[m-1] == BASE1) m--;
				*a0 = - *a0;
				a = a0;
				i = m;
				while (--i) {
					a++;
					*a = ~*a;
				}
			}
		} else {				/* abs(a - b) case */
			while (i && *a++ == *b++) i--;
			q = n - i;