zfunc.c

Calc Software Package for Number Calc

			if (m == n) {		/* a and b same length */
				if (i) {	 /* a not equal to b */
					while (m && a0[m-1] == b0[m-1]) m--;
					if (a0[m-1] < b0[m-1]) {
						/* Swapping since a < b */
						a = a0;
						a0 = b0;
						b0 = a;
						k = j;
					}
					a = a0 + q;
					b = b0 + q;
					i = m - q;
					f = 0;
					while (i--) {
						f = (FULL) *a - *b++ - f;
						*a++ = (HALF) f;
						f >>= BASEB;
						f = -f & BASE1;
					}
				}
			} else {		/* a has more digits than b */
				a = a0 + q;
				b = b0 + q;
				i = n - q;
				f = 0;
				while (i--) {
					f = (FULL) *a - *b++ - f;
					*a++ = (HALF) f;
					f >>= BASEB;
					f = -f & BASE1;
				}
				if (f) { while (!*a) *a++ = BASE1;
					(*a)--;
				}
			}
			a0 += q;
			m -= q;
			while (m && !*a0) { /* Removing trailing zeros */
				m--;
				a0++;
			}
		}
		while (m && !a0[m-1]) m--;	/* Removing leading zeros */
	}
	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);
	if (o) memset(gcd.v, 0, o * sizeof(HALF));	/* Common zero digits */
	if (p) {			/* Left shift for common zero bits */
		i = n;
		f = 0;
		b = b0;
		a = gcd.v + o;
		while (i--) {
			f = (FULL) *b++ << p | f;
			*a++ = (HALF) f;
			f >>= BASEB;
		}
		if (f) {
			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;
}

/*
 * 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(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	ZVALUE temp1, temp2;

	zgcd(z1, z2, &temp1);
	zequo(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(ZVALUE z1, ZVALUE z2)
{
	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, (FULL)3 * 5 * 7 * 11 * 13);
	rem2 = zmodi(z2, (FULL)3 * 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, (FULL)17 * 19 * 23);
	rem2 = zmodi(z2, (FULL)17 * 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 integer floor of the log of an integer to a specified base.
 * The signs of the integers and base are ignored.
 * Example:  zlog(123456, 10) = 5.
 */
long
zlog(ZVALUE z, ZVALUE base)
{
	ZVALUE *zp;			/* current square */
	long power;			/* current power */
	ZVALUE temp;			/* temporary */
	ZVALUE squares[32];		/* table of squares of base */

	/* ignore signs */

	z.sign = 0;
	base.sign = 0;

	/*
	 * Make sure that the numbers are nonzero and the base is > 1
	 */
	if (ziszero(z) || ziszero(base) || zisone(base)) {
		math_error("Zero or too small argument argument for zlog!!!");
		/*NOTREACHED*/
	}

	/*
	 * Some trivial cases.
	 */
	power = zrel(z, base);
	if (power <= 0)
		return (power + 1);

	/* base - power of two */
	if (zisonebit(base))
		return (zhighbit(z) / zlowbit(base));

	/* base = 10 */
	if (base.len == 1 && base.v[0] == 10)
		return zlog10(z, NULL);
	/*
	 * Now loop by squaring the base each time, and see whether or
	 * not each successive square is still smaller than the number.
	 */
	zp = &squares[0];
	*zp = base;
	while (zp->len * 2 - 1 <= z.len  && zrel(z, *zp) > 0) {
		/* while square not too large */
		zsquare(*zp, zp + 1);
		zp++;
	}

	/*
	 * Now back down the squares,
	 */
	power = 0;
	for (; zp > squares; zp--) {
		if (zrel(z, *zp) >= 0) {
			zquo(z, *zp, &temp, 0);
			if (power)
				zfree(z);
			z = temp;
			power++;
		}
		zfree(*zp);
		power <<= 1;
	}
	if (zrel(z, *zp) >= 0)
		power++;
	if (power > 1)
		zfree(z);
	return power;
}


/*
 * Return the integral log base 10 of a number.
 *
 * If was_10_power != NULL, then this flag is set to TRUE if the
 * value was a power of 10, FALSE otherwise.
 */
long
zlog10(ZVALUE z, BOOL *was_10_power)
{
	ZVALUE *zp;			/* current square */
	long power;			/* current power */
	ZVALUE temp;			/* temporary */
	ZVALUE pow10;			/* power of 10 */
	FLAG rel;			/* relationship */
	int i;

	if (ziszero(z)) {
		math_error("Zero argument argument for zlog10");
		/*NOTREACHED*/
	}

	/* Ignore sign of z */
	z.sign = 0;

	/* preload power10 table if missing */
	if (power10 == NULL) {
		long v;

		/* determine power10 table size */
		for (v=1, max_power10_exp=0;
		     v <= (long)(MAXLONG/10L);
		     v *= 10L, ++max_power10_exp) {
		}

		/* create power10 table */
		power10 = malloc(sizeof(long) * (max_power10_exp+1));
		if (power10 == NULL) {
			math_error("cannot malloc power10 table");
			/*NOTREACHED*/
		}

		/* load power10 table */
		for (i=0, v = 1L; i <= max_power10_exp; ++i, v *= 10L) {
			power10[i] = v;
		}
	}

	/* assume not a power of ten unless we find out otherwise */
	if (was_10_power != NULL) {
		*was_10_power = FALSE;
	}

	/* quick exit for small values */
	if (! zgtmaxlong(z)) {
		long value = ztolong(z);

		for (i=0; i <= max_power10_exp; ++i) {
			if (value == power10[i]) {
				if (was_10_power != NULL) {
					*was_10_power = TRUE;
				}
				return i;
			} else if (value < power10[i]) {
				return i-1;
			}
		}
	}

	/*
	 * Loop by squaring the base each time, and see whether or
	 * not each successive square is still smaller than the number.
	 */
	zp = &_tenpowers_[0];
	*zp = _ten_;
	while (((zp->len * 2) - 1) <= z.len) {	/* while square not too large */
		if (zp >= &_tenpowers_[TEN_MAX]) {
			math_error("Maximum storable power of 10 reached!");
			/*NOTREACHED*/
		}
		if (zp[1].len == 0)
			zsquare(*zp, zp + 1);
		zp++;
	}

	/*
	 * Now back down the squares, and multiply them together to see
	 * exactly how many times the base can be raised by.
	 */
	/* find the tenpower table entry < z */
	do {
		rel = zrel(*zp, z);
		if (rel == 0) {
			/* quick exit - we match a tenpower entry */
			if (was_10_power != NULL) {
				*was_10_power = TRUE;
			}
			return (1L << (zp - _tenpowers_));
		}
	} while (rel > 0 && --zp >= _tenpowers_);
	if (zp < _tenpowers_) {
		math_error("fell off bottom of tenpower table!");
		/*NOTREACHED*/
	}

	/* the tenpower value is now our starting comparison value */
	zcopy(*zp, &pow10);
	power = (1L << (zp - _tenpowers_));

	/* try to build up a power of 10 from tenpower table entries */
	while (--zp >= _tenpowers_) {

		/* try the next lower tenpower value */
		zmul(pow10, *zp, &temp);
		rel = zrel(temp, z);
		if (rel == 0) {
			/* exact power of 10 match */
			power += (1L << (zp - _tenpowers_));
			if (was_10_power != NULL) {
				*was_10_power = TRUE;
			}
			return power;

		/* ignore this entry if we went too large */
		} else if (rel > 0) {
			zfree(temp);

		/* otherwise increase power and keep going */
		} else {
			power += (1L << (zp - _tenpowers_));
			zfree(pow10);
			pow10 = 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(ZVALUE z1, ZVALUE z2)
{
	long count;		/* number of factors removed */
	ZVALUE tmp;		/* ignored return value */

	if (ziszero(z1) || ziszero(z2) || zisunit(z2))
		return 0;
	count = zfacrem(z1, z2, &tmp);
	zfree(tmp);
	return count;
}


/*
 * Remove all occurrences 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(ZVALUE z1, ZVALUE z2, ZVALUE *rem)
{
	register ZVALUE *zp;		/* current square */
	long count;			/* total count of divisions */
	long worth;			/* worth of current square */
	long lowbit;			/* for zlowbit(z2) */
	ZVALUE temp1, temp2, temp3;	/* temporaries */
	ZVALUE squares[32];		/* table of squares of factor */

	z1.sign = 0;
	z2.sign = 0;
	/*
	 * Reject trivial cases.
	 */
	if ((z1.len < z2.len) || (zisodd(z1) && ziseven(z2)) ||
		ziszero(z2) || zisone(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)) {
		lowbit = zlowbit(z2);
		count = zlowbit(z1) / lowbit;
		rem->v = alloc(z1.len);
		rem->len = z1.len;
		rem->sign = 0;
		zcopyval(z1, *rem);
		zshiftr(*rem, count * lowbit);
		ztrim(rem);
		return count;
	}
	/*
	 * See if the factor goes in even once.
	 */
	zdiv(z1, z2, &temp1, &temp2, 0);
	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, 0);
		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.
	 */
	/*
	 * We prevent the zp pointer from walking behind squares
	 * by stopping one short of the end and running the loop one
	 * more time.
	 *
	 * We could stop the loop with just zp >= squares, but stopping
	 * short and running the loop one last time manually helps make
	 * code checkers such as insure happy.
	 */
	for (; zp > squares; zp--, worth /= 2) {
		if (zp->len <= z1.len) {
			zdiv(z1, *zp, &temp1, &temp2, 0);
			if (ziszero(temp2)) {
				temp3 = z1;
				z1 = temp1;
				temp1 = temp3;
				count += worth;
			}
			zfree(temp1);
			zfree(temp2);
		}
		if (zp != squares)
			zfree(*zp);
	}
	/* run the loop manually one last time */
	if (zp == squares) {
		if (zp->len <= z1.len) {
			zdiv(z1, *zp, &temp1, &temp2, 0);
			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.  Returns the number
 * of divisions made, and if this is positive, stores result at res.
 */