zmath.c

来自「Calc Software Package for Number Calc」· C语言 代码 · 共 2,000 行 · 第 1/3 页

 */
void
zreduce(ZVALUE z1, ZVALUE z2, ZVALUE *z1res, ZVALUE *z2res)
{
	ZVALUE tmp;

	if (zisabsleone(z1) || zisabsleone(z2))
		tmp = _one_;
	else
		zgcd(z1, z2, &tmp);
	if (zisunit(tmp)) {
		zcopy(z1, z1res);
		zcopy(z2, z2res);
	} else {
		zequo(z1, tmp, z1res);
		zequo(z2, tmp, z2res);
	}
	zfree(tmp);
}



/*
 * Compute the quotient and remainder for division of an integer by an
 * integer, rounding criteria determined by rnd.  Returns the sign of
 * the remainder.
 */
long
zdiv(ZVALUE z1, ZVALUE z2, ZVALUE *quo, ZVALUE *rem, long rnd)
{
	register HALF *a, *b;
	HALF s, u;
	HALF *A, *B, *a1, *b0;
	FULL f, g, h, x;
	BOOL adjust, onebit;
	LEN m, n, len, i, p, j1, j2, k;
	long t, val;

	if (ziszero(z2)) {
		math_error("Division by zero in zdiv");
		/*NOTREACHED*/
	}
	m = z1.len;
	n = z2.len;
	B = z2.v;
	s = 0;
	if (m < n) {
		A = alloc(n + 1);
		memcpy(A, z1.v, m * sizeof(HALF));
		for (i = m; i <= n; i++)
			A[i] = 0;
		a1 = A + n;
		len = 1;
		goto done;
	}
	A = alloc(m + 2);
	memcpy(A, z1.v, m * sizeof(HALF));
	A[m] = 0;
	A[m + 1] = 0;
	len = m - n + 1;	/* quotient length will be len or len +/- 1 */
	a1 = A + n;		/* start of digits for quotient */
	b0 = B;
	p = n;
	while (!*b0) {		/* b0: working start for divisor */
		++b0;
		--p;
	}
	if (p == 1) {
		u = *b0;
		if (u == 1) {
			for (; m >= n; m--)
				A[m] = A[m - 1];
			A[m] = 0;
			goto done;
		}
		f = 0;
		a = A + m;
		i = len;
		while (i--) {
			f = f << BASEB | *--a;
			a[1] = (HALF)(f / u);
			f = f % u;
		}
		*a = (HALF)f;
		m = n;
		goto done;
	}
	f = B[n - 1];
	k = 1;
	while (f >>= 1)
		k++;		/* k: number of bits in top divisor digit */
	j1 = BASEB - k;
	j2 = BASEB + j1;
	h = j1 ? ((FULL) B[n - 1] << j1 | B[n - 2] >> k) : B[n-1];
	onebit = (BOOL)((B[n - 2] >> (k - 1)) & 1);
	m++;
	while (m > n) {
		m--;
		f = (FULL) A[m] << j2 | (FULL) A[m - 1] << j1;
		if (j1) f |= A[m - 2] >> k;
		if (s) f = ~f;
		x = f / h;
		if (x) {
			if (onebit && x > TOPHALF + f % h)
				x--;
			a = A + m - p;
			b = b0;
			u = 0;
			i = p;
			if (s) {
				while (i--) {
					f = (FULL) *a + u + x * *b++;
					*a++ = (HALF) f;
					u = (HALF) (f >> BASEB);
				}
				s = *a + u;
				A[m] = (HALF) (~x + !s);
			} else {
				while (i--) {
					f = (FULL) *a - u - x * *b++;
					*a++ = (HALF) f;
					u = -(HALF)(f >> BASEB);
				}
				s = *a - u;
				A[m] = (HALF)(x + s);
			}
		}
		else
			A[m] = s;
	}
done:	while (m > 0 && A[m - 1] == 0)
		m--;
	if (m == 0 && s == 0) {
		*rem = _zero_;
		val = 0;
		if (a1[len - 1] == 0)
			len--;
		if (len == 0) {
			*quo = _zero_;
		} else {
			quo->len = len;
			quo->v = alloc(len);
			memcpy(quo->v, a1, len * sizeof(HALF));
			quo->sign = z1.sign ^ z2.sign;
		}
		freeh(A);
		return val;
	}
	if (rnd & 8)
		adjust = (((*a1 ^ rnd) & 1) ? TRUE : FALSE);
	else
		adjust = (((rnd & 1) ^ z1.sign ^ z2.sign) ? TRUE : FALSE);
	if (rnd & 2)
		adjust ^= z1.sign ^ z2.sign;
	if (rnd & 4)
		adjust ^= z2.sign;
	if (rnd & 16) {			/* round-off case */
		a = A + n;
		b = B + n;
		i = n + 1;
		f = g = 0;
		t = -1;
		if (s) {
			while (--i > 0) {
				g = (FULL) *--a + (*--b >> 1 | f);
				f = *b & 1 ? TOPHALF : 0;
				if (g != BASE1)
					break;
			}
			if (g == BASE && f == 0) {
				while ((--i > 0) && ((*--a | *--b) == 0));
				t = (i > 0);
			} else if (g >= BASE) {
				t = 1;
			}
		} else {
			while (--i > 0) {
				g = (FULL) *--a - (*--b >> 1 | f);
				f = *b & 1 ? TOPHALF : 0;
				if (g != 0)
					break;
			}
			if (g > 0 && g < BASE)
				t = 1;
			else if (g == 0 && f == 0)
				t = 0;
		}
		if (t)
			adjust = (t > 0);
	}
	if (adjust) {
		i = len;
		a = a1;
		while (i > 0 && *a == BASE1) {
			i--;
			*a++ = 0;
		}
		(*a)++;
		if (i == 0)
			len++;
	}
	if (s && adjust) {
		i = 0;
		while (A[i] == 0)
			i++;
		A[i] = -A[i];
		while (++i < n)
			A[i] = ~A[i];
		m = n;
		while (A[m - 1] == 0)
			m--;
	}
	if (!s && adjust) {
		a = A;
		b = B;
		i = n;
		u = 0;
		while (i--) {
			f = (FULL) *b++ - *a - (FULL) u;
			*a++ = (HALF) f;
			u = -(HALF)(f >> BASEB);
		}
		m = n;
		while (A[m - 1] == 0)
			m--;
	}
	if (s && !adjust) {
		a = A;
		b = B;
		i = n;
		f = 0;
		while (i--) {
			f = (FULL) *b++ + *a + (f >> BASEB);
			*a++ = (HALF) f;
		}
		m = n;
		while (A[m-1] == 0)
			m--;
	}
	rem->len = m;
	rem->v = alloc(m);
	memcpy(rem->v, A, m * sizeof(HALF));
	rem->sign = z1.sign ^ adjust;
	val = rem->sign ? -1 : 1;
	if (a1[len - 1] == 0)
		len--;
	if (len == 0) {
		*quo = _zero_;
	} else {
		quo->len = len;
		quo->v = alloc(len);
		memcpy(quo->v, a1, len * sizeof(HALF));
		quo->sign = z1.sign ^ z2.sign;
	}
	freeh(A);
	return val;
}


/*
 * Compute and store at a specified location the integer quotient
 * of two integers, the type of rounding being determined by rnd.
 * Returns the sign of z1/z2 - calculated quotient.
 */
long
zquo(ZVALUE z1, ZVALUE z2, ZVALUE *res, long rnd)
{
	ZVALUE tmp;
	long val;

	val = zdiv(z1, z2, res, &tmp, rnd);
	if (z2.sign)
		val = -val;
	zfree(tmp);
	return val;
}


/*
 * Compute and store at a specified location the remainder for
 * division of an integer by an integer, the type of rounding
 * used being determined by rnd.  Returns the sign of the remainder.
 */
long
zmod(ZVALUE z1, ZVALUE z2, ZVALUE *res, long rnd)
{
	ZVALUE tmp;
	long val;

	val = zdiv(z1, z2, &tmp, res, rnd);
	zfree(tmp);
	return val;
}


/*
 * Computes the quotient z1/z2 on the assumption that this is exact.
 * There is no thorough check on the exactness of the division
 * so this should not be called unless z1/z2 is an integer.
 */
void
zequo(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	LEN i, j, k, len, m, n, o, p;
	HALF *a, *a0, *A, *b, *B, u, v, w, x;
	FULL f, g;

	if (ziszero(z1)) {
		*res = _zero_;
		return;
	}
	if (ziszero(z2)) {
		math_error("Division by zero");
		/*NOTREACHED*/
	}
	if (zisone(z2)) {
		zcopy(z1, res);
		return;
	}
	if (zhighbit(z1) < zhighbit(z2)) {
		math_error("Bad call to zequo");
		/*NOTREACHED*/
	}
	B = z2.v;
	o = 0;
	while (!*B) {
		++B;
		++o;
	}
	m = z1.len - o;
	n = z2.len - o;
	len = m - n + 1;		/* Maximum length of quotient */
	v = *B;
	A = alloc(len+1);
	memcpy(A, z1.v + o, len * sizeof(HALF));
	A[len] = 0;
	if (n == 1) {
		if (v > 1) {
			a = A + len;
			i = len;
			f = 0;
			while (i--) {
				f = f << BASEB | *--a;
				*a = (HALF)(f / v);
				f %= v;
			}
		}
	} else {
		k = 0;
		while (!(v & 1)) {
			k++;
			v >>= 1;
		}
		j = BASEB - k;
		if (n > 1 && k > 0)
			v |= B[1] << j;
		u = v - 1;
		w = x = 1;
		while (u) {	/* To find w = inverse of v modulo BASE */
			do {
				v <<= 1;
				x <<= 1;
			}
			while (!(u & x));
			u += v;
			w |= x;
		}
		a0 = A;
		p = len;
		while (p > 1) {
			if (!*a0) {
				while (!*++a0 && p > 1)
					p--;
				--a0;
			}
			if (p == 1)
				break;
			x = k ? w * (*a0 >> k | a0[1] << j) : w * *a0;
			g = x;
			if (x) {
				a = a0;
				b = B;
				u = 0;
				i = n;
				if (i > p)
					i = p;
				while (i--) {
					f = (FULL) *a - g * *b++ - (FULL) u;
					*a++ = (HALF)f;
					u = -(HALF)(f >> BASEB);
				}
				if (u && p > n) {
					i = p - n;
					while (u && i--) {
						f = (FULL) *a - u;
						*a++ = (HALF) f;
						u = -(HALF)(f >> BASEB);
					}
				}
			}
			*a0++ = x;
			p--;
		}
		if (k == 0) {
			*a0 = w * *a0;
		} else {
			u = (HALF)(w * *a0) >> k;
			x = (HALF)(((FULL) z1.v[z1.len - 1] << BASEB
				| z1.v[z1.len - 2])
				/((FULL) B[n-1] << BASEB | B[n-2]));
			if ((x ^ u) & 1) x--;
			*a0 = x;
		}
	}
	if (!A[len - 1]) len--;
	res->v = A;
	res->len = len;
	res->sign = z1.sign != z2.sign;
}



/*
 * Return the quotient and remainder of an integer divided by a small
 * number.  A nonzero remainder is only meaningful when both numbers
 * are positive.
 */
long
zdivi(ZVALUE z, long n, ZVALUE *res)
{
	HALF *h1, *sd;
	FULL val;
	HALF divval[2];
	ZVALUE div;
	ZVALUE dest;
	LEN len;

	if (n == 0) {
		math_error("Division by zero");
		/*NOTREACHED*/
	}
	if (ziszero(z)) {
		*res = _zero_;
		return 0;
	}
	if (n < 0) {
		n = -n;
		z.sign = !z.sign;
	}
	if (n == 1) {
		zcopy(z, res);
		return 0;
	}
	/*
	 * If the division is by a large number, then call the normal
	 * divide routine.
	 */
	if (n & ~BASE1) {
		div.sign = 0;
		div.v = divval;
		divval[0] = (HALF) n;
#if LONG_BITS > BASEB
		divval[1] = (HALF)(((FULL) n) >> BASEB);
		div.len = 2;
#else
		div.len = 1;
#endif
		zdiv(z, div, res, &dest, 0);
		n = ztolong(dest);
		zfree(dest);
		return n;
	}
	/*
	 * Division is by a small number, so we can be quick about it.
	 */
	len = z.len;
	dest.sign = z.sign;
	dest.len = len;
	dest.v = alloc(len);
	h1 = z.v + len;
	sd = dest.v + len;
	val = 0;
	while (len--) {
		val = ((val << BASEB) + ((FULL) *--h1));
		*--sd = (HALF)(val / n);
		val %= n;
	}
	zquicktrim(dest);
	*res = dest;
	return (long)val;
}



/*
 * Calculate the mod of an integer by a small number.
 * This is only defined for positive moduli.
 */
long
zmodi(ZVALUE z, long n)
{
	register HALF *h1;
	FULL val;
	HALF divval[2];
	ZVALUE div;
	ZVALUE temp;
	long len;

	if (n == 0) {
		math_error("Division by zero");
		/*NOTREACHED*/
	}
	if (n < 0) {
		math_error("Non-positive modulus");
		/*NOTREACHED*/
	}
	if (ziszero(z) || (n == 1))
		return 0;
	if (zisone(z))
		return 1;
	/*
	 * If the modulus is by a large number, then call the normal
	 * modulo routine.
	 */
	if (n & ~BASE1) {
		div.sign = 0;
		div.v = divval;
		divval[0] = (HALF) n;
#if LONG_BITS > BASEB
		divval[1] = (HALF)(((FULL) n) >> BASEB);
		div.len = 2;
#else
		div.len = 1;
#endif
		zmod(z, div, &temp, 0);
		n = ztolong(temp);
		zfree(temp);
		return n;
	}
	/*
	 * The modulus is by a small number, so we can do this quickly.
	 */
	len = z.len;
	h1 = z.v + len;
	val = 0;
	while (len-- > 0)
		val = ((val << BASEB) + ((FULL) *--h1)) % n;
	if (val && z.sign)
		val = n - val;
	return (long)val;
}


/*
 * Return whether or not one number exactly divides another one.
 * Returns TRUE if division occurs with no remainder.
 * z1 is the number to be divided by z2.
 */

BOOL
zdivides(ZVALUE z1, ZVALUE z2)
{
	LEN i, j, k, m, n;
	HALF u, v, w, x;
	HALF *a, *a0, *A, *b, *B, *c, *d;
	FULL f;
	BOOL ans;

	if (zisunit(z2) || ziszero(z1)) return TRUE;
	if (ziszero(z2)) return FALSE;

	m = z1.len;
	n = z2.len;
	if (m < n) return FALSE;

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

	while (!*d) {
		if (*c++) return FALSE;
		d++;
		m--;
		n--;
	}

	j = 0;
	u = *c;
	v = *d;
	while (!(v & 1)) {	/* Counting and checking zero bits */
		if (u & 1) return FALSE;
		u >>= 1;
		v >>= 1;
		j++;
	}

	if (n == 1 && v == 1) return TRUE;
	if (!*c) {			/* Removing any further zeros */
		while(!*++c)
			m--;
		c--;
	}

	if (m < n) return FALSE;

	if (j) {
		B = alloc((LEN)n);	/* Array for shifted z2 */
		d += n;
		b = B + n;
		i = n;
		f = 0;
		while(i--) {
			f = f << BASEB | *--d;
			*--b = (HALF)(f >> j);
		}
		if (!B[n - 1]) n--;
	}
	else B = d;
	u = *B;
	v = x = 1;
	w = 0;
	while (x) {			/* Finding minv(*B, BASE) */
		if (v & x) {
			v -= x * u;
			w |= x;
		}
		x <<= 1;
	}

	A = alloc((LEN)(m + 2));	/* Main working array */
	memcpy(A, c, m * sizeof(HALF));
	A[m + 1] = 1;
	A[m] = 0;

	k = m - n + 1;			/* Length of presumed quotient */

	a0 = A;

	while (k--) {
		if (*a0) {
			x = w * *a0;
			a = a0;
			b = B;
			i = n;
			u = 0;
			while (i--) {
				f = (FULL) *a - (FULL) x * *b++ - u;
				*a++ = (HALF)f;
				u = -(HALF)(f >> BASEB);
			}
			f = (FULL) *a - u;
			*a++ = (HALF)f;
			u = -(HALF)(f >> BASEB);
			if (u) {
				while (*a == 0) *a++ = BASE1;
				(*a)--;
			}
		}
		a0++;
	}
	ans = FALSE;
	if (A[m + 1]) {
		a = A + m;
		while (--n && !*--a);
		if (!n) ans = TRUE;
	}
	freeh(A);
	if (j) freeh(B);