zfunc.c

Calc Software Package for Number Calc

long
zgcdrem(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	ZVALUE tmp1, tmp2;
	long count, onecount;
	long sh;

	if (ziszero(z1) || ziszero(z2)) {
		math_error("Zero argument in call to zgcdrem!!!");
		/*NOTREACHED*/
	}
	/*
	 * Begin by taking the gcd for the first time.
	 * If the number is already relatively prime, then we are done.
	 */
	z1.sign = 0;
	z2.sign = 0;
	if (zisone(z2))
		return 0;
	if (zisonebit(z2)) {
		sh = zlowbit(z1);
		if (sh == 0)
			return 0;
		zshift(z1, -sh, res);
		return 1 + (sh - 1)/zlowbit(z2);
	}
	if (zisonebit(z1)) {
		if (zisodd(z2))
			return 0;
		*res = _one_;
		return zlowbit(z1);
	}

	zgcd(z1, z2, &tmp1);
	if (zisunit(tmp1) || ziszero(tmp1))
		return 0;
	zequo(z1, tmp1, &tmp2);
	count = 1;
	z1 = tmp2;
	z2 = tmp1;
	/*
	 * Now keep alternately taking the gcd and removing factors until
	 * the gcd becomes one.
	 */
	while (!zisunit(z2)) {
		onecount = zfacrem(z1, z2, &tmp1);
		if (onecount) {
			count += onecount;
			zfree(z1);
			z1 = tmp1;
		}
		zgcd(z1, z2, &tmp1);
		zfree(z2);
		z2 = tmp1;
	}
	*res = z1;
	return count;
}


/*
 * Return the number of digits (base 10) in a number, ignoring the sign.
 */
long
zdigits(ZVALUE z1)
{
	long count, val;

	z1.sign = 0;
	if (!zge16b(z1)) {	/* do small numbers ourself */
		count = 1;
		val = 10;
		while (*z1.v >= (HALF)val) {
			count++;
			val *= 10;
		}
		return count;
	}
	return (zlog10(z1, NULL) + 1);
}


/*
 * Return the single digit at the specified decimal place of a number,
 * where 0 means the rightmost digit.  Example:	 zdigit(1234, 1) = 3.
 */
long
zdigit(ZVALUE z1, long n)
{
	ZVALUE tmp1, tmp2;
	long 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, 0);
	res = zmodi(tmp2, 10L);
	zfree(tmp1);
	zfree(tmp2);
	return res;
}


/*
 * z is to be a nonnegative integer
 * If z is the square of a integer stores at dest the square root of z;
 *	otherwise stores at z an integer differing from the square root
 *	by less than 1.	 Returns the sign of the true square root minus
 *	the calculated integer.	 Type of rounding is determined by
 *	rnd as follows: rnd = 0 gives round down, rnd = 1
 *	rounds up, rnd = 8 rounds to even integer, rnd = 9 rounds to odd
 *	integer, rnd = 16 rounds to nearest integer.
 */
FLAG
zsqrt(ZVALUE z, ZVALUE *dest, long rnd)
{
	HALF *a, *A, *b, *a0, u;
	int i, j, j1, j2, k, k1, m, m0, m1, n, n0, o;
	FULL d, e, f, g, h, s, t, x, topbit;
	int remsign;
	BOOL up, onebit;
	ZVALUE sqrt;

	if (z.sign) {
		math_error("Square root of negative number");
		/*NOTREACHED*/
	}
	if (ziszero(z)) {
		*dest = _zero_;
		return 0;
	}
	m0 = z.len;
	o = m0 & 1;
	m = m0 + o;		/* m is smallest even number >= z.len */
	n0 = n = m / 2;
	f = z.v[z.len - 1];
	k = 1;
	while (f >>= 2)
		k++;
	if (!o)
		k += BASEB/2;
	j = BASEB - k;
	m1 = m;
	if (k == BASEB) {
		m1 += 2;
		n0++;
	}
	A = alloc(m1);
	A[m1] = 0;
	a0 = A + n0;
	memcpy(A, z.v, m0 * sizeof(HALF));
	if (o)
		A[m - 1] = 0;
	if (n == 1) {
		if (j)
			f = (FULL) A[1] << j | A[0] >> k;
		else
			f = A[1];
		g = (FULL) A[0] << (j + BASEB);
		d = e = topbit = (FULL)1 << (k - 1);
	} else {
		if (j)
			f = (FULL) A[m-1] << (j + BASEB) | (FULL) A[m-2] << j |
				A[m-3] >> k;
		else
			f = (FULL) A[m-1] << BASEB | A[m-2];
		g = (FULL) A[m-3] << (j + BASEB) | (FULL) A[m-4] << j;
		d = e = topbit = (FULL)1 << (BASEB + k - 1);
	}

	s = (f & topbit);
	f <<= 1;
	if (g & TOPFULL)
		f++;
	g <<= 1;
	if (s) {
		f -= 4 * d;
		e = 2 * d - 1;
	}
	else
		f -= d;
	while (d >>= 1) {
		if (!(s | f | g))
			break;
		while (d && (f & topbit) == s) {
			d >>= 1;
			f <<= 1;
			if (g & TOPFULL)
				f++;
			g <<= 1;
		}
		if (d == 0)
			break;
		if (s)
			f += e + 1;
		else
			f -= e;
		t = f & topbit;
		f <<= 1;
		if (g & TOPFULL)
			f++;
		g <<= 1;
		if (t == 0 && f < d)
			t = topbit;
		f -= d;
		if (s)
			e -= d - !t;
		else
			e += d - (t > 0);
		s = t;
	}
	if (n0 == 1) {
		A[1] = (HALF)e;
		A[0] = (HALF)f;
		m = 1;
		goto done;
	}
	if (n0 == 2) {
		A[3] = (HALF)(e >> BASEB);
		A[2] = (HALF)e;
		A[1] = (HALF)(f >> BASEB);
		A[0] = (HALF)f;
		m = 2;
		goto done;
	}
	u = (HALF)(s ? BASE1 : 0);
	if (k < BASEB) {
		A[m1 - 1] = (HALF)(e >> (BASEB - 1));
		A[m1 - 2] = ((HALF)(e << 1) | (HALF)(s > 0));
		A[m1 - 3] = (HALF)(f >> BASEB);
		A[m1 - 4] = (HALF)f;
		m = m1 - 2;
		k1 = k + 1;
	} else {
		A[m1 - 1] = 1;
		A[m1 - 2] = (HALF)(e >> (BASEB - 1));
		A[m1 - 3] = ((HALF)(e << 1) | (HALF)(s > 0));
		A[m1 - 4] = u;
		A[m1 - 5] = (HALF)(f >> BASEB);
		A[m1 - 6] = (HALF)f;
		m = m1 - 3;
		k1 = 1;
	}
	h = e >> k;
	onebit = ((e & ((FULL)1 << (k - 1))) ? TRUE : FALSE);
	j2 = BASEB - k1;
	j1 = BASEB + j2;
	while (m > n0) {
		a = A + m - 1;
		if (j2)
			f = (FULL) *a << j1 | (FULL) a[-1] << j2 | a[-2] >> k1;
		else
			f = (FULL) *a << BASEB | a[-1];
		if (u)
			f = ~f;
		x = f / h;
		if (x) {
			if (onebit && x > 2 * (f % h) + 2)
				x--;
			b = a + 1;
			i = m1 - m;
			a -= i + 1;
			if (u) {
				f = *a + x * (BASE - x);
				*a++ = (HALF)f;
				u = (HALF)(f >> BASEB);
				while (i--) {
					f = *a + x * *b++ + u;
					*a++ = (HALF)f;
					u = (HALF)(f >> BASEB);
				}
				u += *a;
				x = ~x + !u;
				if (!(x & TOPHALF))
					a[1] -= 1;
			} else {
				f = *a - x * x;
				*a++ = (HALF)f;
				u = -(HALF)(f >> BASEB);
				while (i--) {
					f = *a - x * *b++ - u;
					*a++ = (HALF)f;
					u = -(HALF)(f >> BASEB);
				}
				u = *a - u;
				x = x + u;
				if (x & TOPHALF)
					a[1] |= 1;
			}
			*a = ((HALF)(x << 1) | (HALF)(u > 0));
		} else {
			*a = u;
		}
		m--;
		if (*--a == u) {
			while (m > 1 && *--a == u)
				m--;
		}
	}
	i = n;
	a = a0;
	while (i--) {
		*a >>= 1;
		if (a[1] & 1) *a |= TOPHALF;
		a++;
	}
	s = u;
done:	if (s == 0) {
		while (m > 0 && A[m - 1] == 0)
			m--;
		if (m == 0) {
			remsign = 0;
			sqrt.v = alloc(n);
			sqrt.len = n;
			sqrt.sign = 0;
			memcpy(sqrt.v, a0, n * sizeof(HALF));
			freeh(A);
			*dest = sqrt;
			return remsign;
		}
	}
	if (rnd & 16) {
		if (s == 0) {
			if (m != n) {
				up = (m > n);
			} else {
				i = n;
				b = a0 + n;
				a = A + n;
				while (i > 0 && *--a == *--b)
					i--;
				up = (i > 0 && *a > *b);
			}
		} else {
			while (m > 1 && A[m - 1] == BASE1)
				m--;
			if (m != n) {
				up = (m < n);
			} else {
				i = n;
				b = a0 + n;
				a = A + n;
				while (i > 0 && *--a + *--b == BASE1)
					i--;
				up = ((FULL) *a + *b >= BASE);
			}
		}
	}
	else
	if (rnd & 8)
		up = (((rnd ^ *a0) & 1) ? TRUE : FALSE);
	else
		up = ((rnd & 1) ? TRUE : FALSE);
	if (up) {
		remsign = -1;
		i = n;
		a = a0;
		while (i-- && *a == BASE1)
			*a++ = 0;
		if (i >= 0) {
			(*a)++;
		} else {
			n++;
			*a = 1;
		}
	} else {
		remsign = 1;
	}
	sqrt.v = alloc(n);
	sqrt.len = n;
	sqrt.sign = 0;
	memcpy(sqrt.v, a0, n * sizeof(HALF));
	freeh(A);
	*dest = sqrt;
	return remsign;

}

/*
 * 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(ZVALUE z1, ZVALUE z2, ZVALUE *dest)
{
	ZVALUE ztry, quo, old, temp, temp2;
	ZVALUE k1;			/* holds k - 1 */
	int sign;
	long i;
	LEN highbit, k;
	SIUNION sival;

	sign = z1.sign;
	if (sign && ziseven(z2)) {
		math_error("Even root of negative number");
		/*NOTREACHED*/
	}
	if (ziszero(z2) || zisneg(z2)) {
		math_error("Non-positive root");
		/*NOTREACHED*/
	}
	if (ziszero(z1)) {	/* root of zero */
		*dest = _zero_;
		return;
	}
	if (zisunit(z2)) {	/* first root */
		zcopy(z1, dest);
		return;
	}
	if (zge31b(z2)) {	/* humongous root */
		*dest = _one_;
		dest->sign = (BOOL)((HALF)sign);
		return;
	}
	k = (LEN)ztolong(z2);
	highbit = zhighbit(z1);
	if (highbit < k) {	/* too high a root */
		*dest = _one_;
		dest->sign = (BOOL)((HALF)sign);
		return;
	}
	sival.ivalue = k - 1;
	k1.v = &sival.silow;
	/* ignore Saber-C warning #112 - get ushort from uint */
	/*	  ok to ignore on name zroot`sival */
	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;
	ztry.len = (highbit / BASEB) + 1;
	ztry.v = alloc(ztry.len);
	zclearval(ztry);
	ztry.v[ztry.len-1] = ((HALF)1 << (highbit % BASEB));
	ztry.sign = 0;
	old.v = alloc(ztry.len);
	old.len = 1;
	zclearval(old);
	old.sign = 0;
	/*
	 * Main divide and average loop
	 */
	for (;;) {
		zpowi(ztry, k1, &temp);
		zquo(z1, temp, &quo, 0);
		zfree(temp);
		i = zrel(ztry, 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, ztry) == 0)) {
				zfree(quo);
				zfree(old);
				ztry.sign = (BOOL)((HALF)sign);
				zquicktrim(ztry);
				*dest = ztry;
				return;
			}
			old.len = ztry.len;
			zcopyval(ztry, old);
		}
		/* average current try and quotient for the new try */
		zmul(ztry, k1, &temp);
		zfree(ztry);
		zadd(quo, temp, &temp2);
		zfree(temp);
		zfree(quo);
		zquo(temp2, z2, &ztry, 0);
		zfree(temp2);
	}
}


/*
 * Test to see if a number is an exact square or not.
 */
BOOL
zissquare(ZVALUE z)
{
	long n;
	ZVALUE tmp;

	/* negative values are never perfect squares */
	if (zisneg(z)) {
		return FALSE;
	}

	/* ignore trailing zero words */
	while ((z.len > 1) && (*z.v == 0)) {
		z.len--;
		z.v++;
	}

	/* zero or one is a perfect square */
	if (zisabsleone(z)) {
		return TRUE;
	}

	/* check mod 4096 values */
	if (issq_mod4k[(int)(*z.v & 0xfff)] == 0) {
		return FALSE;
	}

	/* must do full square root test now */
	n = !zsqrt(z, &tmp, 0);
	zfree(tmp);
	return (n ? TRUE : FALSE);
}