Initial commit.

1 files changed, 858 insertions, 0 deletions

diff --git a/ic-reals-6.3/base/util.c b/ic-reals-6.3/base/util.c
new file mode 100644
index 0000000..f523cd8
--- /dev/null
+++ b/ic-reals-6.3/base/util.c
@@ -0,0 +1,858 @@
+/*
+ * Copyright (C) 2000, Imperial College
+ *
+ * This file is part of the Imperial College Exact Real Arithmetic Library.
+ * See the copyright notice included in the distribution for conditions
+ * of use.
+ */
+
+#include <stdio.h>
+#include "real.h"
+#include <gmp-impl.h>
+#include <longlong.h>
+#include "real-impl.h"
+
+/*
+ * Utilities used internally by the library
+ */
+/*
+ * This is used only when compiled with -DTRACE=traceOn.
+ */
+int traceOn = 0;
+
+/*
+ * Some temporary big numbers which are shared
+ */
+mpz_t tmpa_z, tmpb_z, tmpc_z, tmpd_z, tmpe_z, tmpf_z;
+mpz_t zero_z;
+
+Matrix bigTmpMat;
+Tensor bigTmpTen;
+
+void
+debugTrace(int b)
+{
+	traceOn = b;
+}
+
+/*
+ * The library shares bignum temporart storage. This needs to change
+ * when garbage collection goes back in.
+ */
+void
+initTmp()
+{
+	mpz_init(tmpa_z);
+	mpz_init(tmpb_z);
+	mpz_init(tmpc_z);
+	mpz_init(tmpd_z);
+	mpz_init(tmpe_z);
+	mpz_init(tmpf_z);
+
+	mpz_init_set_ui(zero_z, 0);
+
+	mpz_init(bigTmpMat[0][0]);
+	mpz_init(bigTmpMat[0][1]);
+	mpz_init(bigTmpMat[1][0]);
+	mpz_init(bigTmpMat[1][1]);
+
+	mpz_init(bigTmpTen[0][0]);
+	mpz_init(bigTmpTen[0][1]);
+	mpz_init(bigTmpTen[1][0]);
+	mpz_init(bigTmpTen[1][1]);
+	mpz_init(bigTmpTen[2][0]);
+	mpz_init(bigTmpTen[2][1]);
+	mpz_init(bigTmpTen[3][0]);
+	mpz_init(bigTmpTen[3][1]);
+}
+
+void
+multVectorPairTimesVector(Vector vec0, Vector vec1, Vector vec)
+{
+	/* ae + cf */
+	mpz_mul(tmpa_z, vec0[0], vec[0]);
+	mpz_mul(tmpe_z, vec1[0], vec[1]);
+	mpz_add(tmpa_z, tmpa_z, tmpe_z);
+
+	/* be + df */
+	mpz_mul(tmpb_z, vec0[1], vec[0]);
+	mpz_mul(tmpe_z, vec1[1], vec[1]);
+	mpz_add(tmpb_z, tmpb_z, tmpe_z);
+
+	/*
+	 * Computed the product, now replace with the original matrix
+	 * with a vector.
+	 */
+	MPZ_SWAP(tmpa_z, vec0[0]);
+	MPZ_SWAP(tmpb_z, vec0[1]);
+
+	mpz_set_ui(vec1[0], 0);
+	mpz_set_ui(vec1[1], 0);
+}
+
+/*
+ * Given vectors (a,b) and (c,d) and the matrix is ((e,f), (g,h)). This
+ * computes the product as if the first two vectors are the columns of a
+ * matrix. The result is a pair of vectors:
+ *
+ *    (ae + cf, be + df) and (ag + ch, bg + dh)
+ *
+ * It is important to note that the result overwrites the original vectors
+ * using the SWAP macro.
+ */
+void
+multVectorPairTimesMatrix(Vector vec0, Vector vec1, Matrix mat)
+{
+	/* ae + cf */
+	mpz_mul(tmpa_z, vec0[0], mat[0][0]);
+	mpz_mul(tmpe_z, vec1[0], mat[0][1]);
+	mpz_add(tmpa_z, tmpa_z, tmpe_z);
+
+	/* be + df */
+	mpz_mul(tmpb_z, vec0[1], mat[0][0]);
+	mpz_mul(tmpe_z, vec1[1], mat[0][1]);
+	mpz_add(tmpb_z, tmpb_z, tmpe_z);
+
+	/* ag + ch */
+	mpz_mul(tmpc_z, vec0[0], mat[1][0]);
+	mpz_mul(tmpe_z, vec1[0], mat[1][1]);
+	mpz_add(tmpc_z, tmpc_z, tmpe_z);
+
+	/* bg + dh */
+	mpz_mul(tmpd_z, vec0[1], mat[1][0]);
+	mpz_mul(tmpe_z, vec1[1], mat[1][1]);
+	mpz_add(tmpd_z, tmpd_z, tmpe_z);
+
+	/*
+	 * Computed the product, now replace with the original vectors
+	 */
+	MPZ_SWAP(tmpa_z, vec0[0]);
+	MPZ_SWAP(tmpb_z, vec0[1]);
+	MPZ_SWAP(tmpc_z, vec1[0]);
+	MPZ_SWAP(tmpd_z, vec1[1]);
+}
+
+/*
+ * Same as the above except the columns of the matrix are (e,f) and (g,h)
+ */
+void
+multVectorPairTimesMatrix_Z(Vector vec0, Vector vec1,
+    mpz_t e, mpz_t f, mpz_t g, mpz_t h)
+{
+	/* ae + cf */
+	mpz_mul(tmpa_z, vec0[0], e);
+	mpz_mul(tmpe_z, vec1[0], f);
+	mpz_add(tmpa_z, tmpa_z, tmpe_z);
+
+	/* be + df */
+	mpz_mul(tmpb_z, vec0[1], e);
+	mpz_mul(tmpe_z, vec1[1], f);
+	mpz_add(tmpb_z, tmpb_z, tmpe_z);
+
+	/* ag + ch */
+	mpz_mul(tmpc_z, vec0[0], g);
+	mpz_mul(tmpe_z, vec1[0], h);
+	mpz_add(tmpc_z, tmpc_z, tmpe_z);
+
+	/* bg + dh */
+	mpz_mul(tmpd_z, vec0[1], g);
+	mpz_mul(tmpe_z, vec1[1], h);
+	mpz_add(tmpd_z, tmpd_z, tmpe_z);
+
+	/*
+	 * Computed the product, now replace with the original vectors
+	 */
+	MPZ_SWAP(tmpa_z, vec0[0]);
+	MPZ_SWAP(tmpb_z, vec0[1]);
+	MPZ_SWAP(tmpc_z, vec1[0]);
+	MPZ_SWAP(tmpd_z, vec1[1]);
+}
+
+/*
+ * Not needed for GMP 3.1
+inline void
+mpz_mul_si(mpz_t x, mpz_t y, int z)
+{
+	if (z >= 0)
+		mpz_mul_ui(x, y, (unsigned) z);
+	else {
+		mpz_mul_ui(x, y, (unsigned) (-z));
+		mpz_neg(x, x);
+	}
+}
+*/
+
+/*
+ * This is the same as above except the entries in the matrix all
+ * fit in a machine word.
+ */
+void
+multVectorPairTimesSmallMatrix(Vector vec0, Vector vec1, SmallMatrix mat)
+{
+	/* ae + cf */
+	mpz_mul_si(tmpa_z, vec0[0], mat[0][0]);
+	mpz_mul_si(tmpe_z, vec1[0], mat[0][1]);
+	mpz_add(tmpa_z, tmpa_z, tmpe_z);
+
+	/* be + df */
+	mpz_mul_si(tmpb_z, vec0[1], mat[0][0]);
+	mpz_mul_si(tmpe_z, vec1[1], mat[0][1]);
+	mpz_add(tmpb_z, tmpb_z, tmpe_z);
+
+	/* ag + ch */
+	mpz_mul_si(tmpc_z, vec0[0], mat[1][0]);
+	mpz_mul_si(tmpe_z, vec1[0], mat[1][1]);
+	mpz_add(tmpc_z, tmpc_z, tmpe_z);
+
+	/* bg + dh */
+	mpz_mul_si(tmpd_z, vec0[1], mat[1][0]);
+	mpz_mul_si(tmpe_z, vec1[1], mat[1][1]);
+	mpz_add(tmpd_z, tmpd_z, tmpe_z);
+
+	/*
+	 * Computed the product, now replace with the original vectors
+	 */
+	MPZ_SWAP(tmpa_z, vec0[0]);
+	MPZ_SWAP(tmpb_z, vec0[1]);
+	MPZ_SWAP(tmpc_z, vec1[0]);
+	MPZ_SWAP(tmpd_z, vec1[1]);
+}
+
+/*
+ * This does the work of taking some number of digits from something bearing
+ * information such as an LFT. There are however, possibilities other than
+ * LFTs. See realSqrt for another example of the use of this function.
+ * The digits are deposited in a digsX structure augmenting any digits
+ * already in the structure.
+ *
+ * The arguments are as follows:
+ *
+ *   digsX - the place where the digits are to be deposited
+ *
+ *   emitDigit - this is a pointer to a function which is used to
+ *     obtain a single digit from the object and compute the residual of
+ *     that object once than digits is emitted.
+ *
+ *   info - this is the object which holds informations such as a vector,
+ *     matrix or tensor and from which we are taking digits. This function 
+ *     does not inspect ``info'' directly. It doesn't care what ``info'' is
+ *     and only passes it to emitDigit
+ *
+ *   digitsNeeded - this is what is says - the number of digits we should try
+ *     to emit. Often this is chosen to be MAXINT, or in other words take
+ *     all the digits possible.
+ *
+ * The function returns the number of digits it has managed to emit (which may
+ * be and is often less than we asked for).
+ */
+int
+emitDigits(DigsX *digsX, edf emitDigit, void *info, int digitsNeeded)
+{
+	Digit d;
+	int word;
+	unsigned char count;
+	int total;
+	bool ok;
+
+	total = 0;
+	ok = TRUE;
+	while (digitsNeeded > 0 && ok) {
+		/*
+		 * We now extract digits from info until we have filled
+		 * a machine word or until no digit can be extracted.
+		 */
+		word = 0;
+		count = 0;
+		while (count < DIGITS_PER_WORD && digitsNeeded > 0 && ok) {
+			ok = (*emitDigit)(info, &d);
+			if (ok) {
+				count++;
+				digitsNeeded--;
+				word = (word << 1) + d;
+			}
+		}
+
+		/*
+		 * If we get here and count > 0, then we need to augment the 
+		 * word stored in the DigsX.
+		 */
+		if (count > 0) {
+			total += count;
+
+			/*
+			 * First check to see if we are already dealing with a large word.
+			 */
+#ifdef PACK_DIGITS
+			if (digsX->count > DIGITS_PER_WORD) {
+#endif
+				mpz_mul_2exp(digsX->word.big, digsX->word.big, count);
+				if (word >= 0) {	/* there is no mpz_add_si function */
+					mpz_add_ui(digsX->word.big, digsX->word.big, word);
+				}
+				else {
+					mpz_sub_ui(digsX->word.big, digsX->word.big, -word);
+				}
+				digsX->count += count;
+#ifdef PACK_DIGITS
+			}
+			else {
+				digsX->count += count;
+	
+				/*
+				 * Now see if we are about to overflow the machine word.
+				 */
+				if (digsX->count > DIGITS_PER_WORD) {
+					mpz_init_set_si(digsX->word.big, digsX->word.small);
+					mpz_mul_2exp(digsX->word.big, digsX->word.big, count);
+					if (word >= 0) {	/* there is no mpz_add_si function */
+						mpz_add_ui(digsX->word.big, digsX->word.big, word);
+					}
+					else {
+						mpz_sub_ui(digsX->word.big, digsX->word.big, -word);
+					}
+				}
+				else
+					digsX->word.small = (digsX->word.small << count) + word;
+			}
+#endif
+		}
+	}
+	return total;
+}
+
+/*
+ * The assumption here is that the number of digits in digsX <= DIGITS_PER_WORD. */
+void
+makeSmallMatrixFromDigits(SmallMatrix mat, DigsX *digsX)
+{
+	int twoPowN;
+
+	if (digsX->count <= 0) {
+		mat[0][0] = 1;
+		mat[0][1] = 0;
+		mat[1][0] = 0;
+		mat[1][1] = 1;
+	}
+	else {
+		twoPowN = 1 << digsX->count;
+		mat[0][0] = twoPowN + digsX->word.small + 1;
+		mat[0][1] = twoPowN - digsX->word.small - 1;
+		mat[1][0] = twoPowN + digsX->word.small - 1;
+		mat[1][1] = twoPowN - digsX->word.small + 1;
+	}
+}
+
+/*
+ * Same as the function above except for big integers.
+ * In this case we assume that the number of digits with which to form
+ * the matrix is > DIGITS_PER_WORD so be warned.
+ */
+void
+makeMatrixFromDigits(Matrix mat, DigsX *digsX)
+{
+	mpz_set_ui(tmpe_z, (unsigned long) 1);
+	mpz_mul_2exp(tmpe_z, tmpe_z, (unsigned long) digsX->count);  /* 2^n */
+	switch (mpz_sgn(digsX->word.big)) {
+		case 0 :   /* not sure I should bother making this a special case */
+			mpz_add_ui(mat[0][0], tmpe_z, (unsigned long) 1);
+			mpz_sub_ui(mat[0][1], tmpe_z, (unsigned long) 1);
+			mpz_sub_ui(mat[1][0], tmpe_z, (unsigned long) 1);
+			mpz_add_ui(mat[1][1], tmpe_z, (unsigned long) 1);
+			break;
+		case 1 :
+		case -1 :
+			mpz_add(mat[0][0], tmpe_z, digsX->word.big);
+			mpz_set(mat[1][0], mat[0][0]);
+			mpz_sub(mat[0][1], tmpe_z, digsX->word.big);
+			mpz_set(mat[1][1], mat[0][1]);
+			mpz_add_ui(mat[0][0], mat[0][0], (unsigned long) 1);
+			mpz_sub_ui(mat[0][1], mat[0][1], (unsigned long) 1);
+			mpz_sub_ui(mat[1][0], mat[1][0], (unsigned long) 1);
+			mpz_add_ui(mat[1][1], mat[1][1], (unsigned long) 1);
+			break;
+		default :
+			Error(FATAL, E_INT, "makeMatrixFromDigits",
+						"bad value from mpz_sign");
+			break;
+	}
+}
+
+/*
+ * This makes a vector ``canonical'' by removing common factors from the
+ * numerator and denominator.
+ */
+void
+canonVector(Vector vec)
+{
+	mpz_gcd(tmpa_z, vec[0], vec[1]);
+	if (mpz_sgn(tmpa_z) != 0) {
+		mpz_divexact(vec[0], vec[0], tmpa_z);
+		mpz_divexact(vec[1], vec[1], tmpa_z);
+	}
+}
+#define MIN_AB(a,b) ((a)>0 ? ((a)>(b) ? (a) : (b)) : ((b)>0 ? (b) : MAXINT))
+#define SIZE_GT_ZERO(size, w)	((size) > 0 ? (w) : 0)
+
+/*
+ * A vector, matrix or tensor is normalized when there are no negative
+ * entries and there are no common factors. In practice, ensuring that
+ * there are no negative entries is dealt with by the functions which 
+ * emit signs and digits. Also, rather than looking for gcd's for all the
+ * entries, we only look for exponents of 2. See Reinhold Heckman's
+ * notes for a justification of this.
+ */
+int
+normalizeVector(Vector vec)
+{
+	mp_limb_t word;
+	mp_size_t size_a, size_b, min_size;
+	int count;
+	int i;
+
+	size_a = ABS(vec[0][0]._mp_size);
+	size_b = ABS(vec[1][0]._mp_size);
+	min_size = MIN_AB(size_a, size_b);
+
+	if (min_size == MAXINT)
+		return 0;
+
+	/*
+	 * The trick with normalization is to find the largest
+	 * exponent of 2 which divides all four entries in the matrix.
+	 * When looking for the least significant bit which is set,
+	 * we ignore those vector entries which are 0.
+	 */
+	for (i = 0; i < min_size; i++) {
+		word = SIZE_GT_ZERO(size_a, vec[0][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_b, vec[1][0]._mp_d[i]);
+		if (word != 0) {
+			count_trailing_zeros(count, word);
+			count = count + (i * mp_bits_per_limb);
+			if (count > 0) {
+				if (size_a > 0)
+					mpz_tdiv_q_2exp(vec[0], vec[0], count);
+				if (size_b > 0)
+					mpz_tdiv_q_2exp(vec[1], vec[1], count);
+			}
+			return count;
+		}
+	}
+	Error(FATAL, E_INT, "normalizeVector", "zero entry with non-zero size");
+	return 0;
+}
+
+int
+normalizeMatrix(Matrix mat)
+{
+	mp_limb_t word;
+	mp_size_t size_a, size_b, size_c, size_d, min_ab, min_cd, min_size;
+	int count;
+	int i;
+	
+	size_a = ABS(mat[0][0][0]._mp_size);
+	size_b = ABS(mat[0][1][0]._mp_size);
+	size_c = ABS(mat[1][0][0]._mp_size);
+	size_d = ABS(mat[1][1][0]._mp_size);
+
+	/*
+	 * GMP respresents zero as _mp_size = 0. Such matrix entries
+	 * can be ignored for the purposes of normalization.
+	 */ 
+	min_ab = MIN_AB(size_a, size_b);
+	min_cd = MIN_AB(size_c, size_d);
+	min_size = MIN(min_ab, min_cd);
+
+	if (min_size == MAXINT)
+		return 0;
+
+	/*
+	 * The trick with normalization is to find the largest
+	 * exponent of 2 which divides all four entries in the matrix.
+	 * When looking for the least significant bit which is set,
+	 * we ignore those matrix entries which are 0.
+	 */
+	for (i = 0; i < min_size; i++) {
+		word = SIZE_GT_ZERO(size_a, mat[0][0][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_b, mat[0][1][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_c, mat[1][0][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_d, mat[1][1][0]._mp_d[i]);
+		if (word != 0) {
+			count_trailing_zeros(count, word);
+			count = count + (i * mp_bits_per_limb);
+			if (count > 0) {
+				if (size_a > 0)
+					mpz_tdiv_q_2exp(mat[0][0], mat[0][0], count);
+				if (size_b > 0)
+					mpz_tdiv_q_2exp(mat[0][1], mat[0][1], count);
+				if (size_c > 0)
+					mpz_tdiv_q_2exp(mat[1][0], mat[1][0], count);
+				if (size_d > 0)
+					mpz_tdiv_q_2exp(mat[1][1], mat[1][1], count);
+			}
+			return count;
+		}
+	}
+	Error(FATAL, E_INT, "normalizeMatrix", "zero entry with non-zero size");
+	return 0;
+}
+
+int
+normalizeTensor(Tensor ten)
+{
+	mp_limb_t word;
+	mp_size_t size_a, size_b, size_c, size_d, min_ab, min_cd, min_size;
+	mp_size_t size_e, size_f, size_g, size_h, min_ef, min_gh;
+	mp_size_t min_abcd, min_efgh;
+	int count;
+	int i;
+	
+	size_a = ABS(ten[0][0][0]._mp_size);
+	size_b = ABS(ten[0][1][0]._mp_size);
+	size_c = ABS(ten[1][0][0]._mp_size);
+	size_d = ABS(ten[1][1][0]._mp_size);
+	size_e = ABS(ten[2][0][0]._mp_size);
+	size_f = ABS(ten[2][1][0]._mp_size);
+	size_g = ABS(ten[3][0][0]._mp_size);
+	size_h = ABS(ten[3][1][0]._mp_size);
+
+	/*
+	 * GMP respresents zero as _mp_size = 0. Such matrix entries
+	 * can be ignored for the purposes of normalization. 
+	 *
+	 * The following are macros and hence we prefer not to nest them.
+	 */ 
+	min_ab = MIN_AB(size_a, size_b);
+	min_cd = MIN_AB(size_c, size_d);
+	min_ef = MIN_AB(size_e, size_f);
+	min_gh = MIN_AB(size_g, size_h);
+	min_abcd = MIN(min_ab, min_cd);
+	min_efgh = MIN(min_ef, min_gh);
+	min_size = MIN(min_abcd, min_efgh);
+
+	if (min_size == MAXINT)
+		return 0;
+
+	/*
+	 * The trick with normalization is to find the largest
+	 * exponent of 2 which divides all four entries in the tensor.
+	 * When looking for the least significant bit which is set,
+	 * we ignore those tensor entries which are 0.
+	 */
+	for (i = 0; i < min_size; i++) {
+		word = SIZE_GT_ZERO(size_a, ten[0][0][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_b, ten[0][1][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_c, ten[1][0][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_d, ten[1][1][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_e, ten[2][0][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_f, ten[2][1][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_g, ten[3][0][0]._mp_d[i])
+			 | SIZE_GT_ZERO(size_h, ten[3][1][0]._mp_d[i]);
+		if (word != 0) {
+			count_trailing_zeros(count, word);
+			count = count + (i * mp_bits_per_limb);
+			if (count > 0) {
+				if (size_a > 0)
+					mpz_tdiv_q_2exp(ten[0][0], ten[0][0], count);
+				if (size_b > 0)
+					mpz_tdiv_q_2exp(ten[0][1], ten[0][1], count);
+				if (size_c > 0)
+					mpz_tdiv_q_2exp(ten[1][0], ten[1][0], count);
+				if (size_d > 0)
+					mpz_tdiv_q_2exp(ten[1][1], ten[1][1], count);
+				if (size_e > 0)
+					mpz_tdiv_q_2exp(ten[2][0], ten[2][0], count);
+				if (size_f > 0)
+					mpz_tdiv_q_2exp(ten[2][1], ten[2][1], count);
+				if (size_g > 0)
+					mpz_tdiv_q_2exp(ten[3][0], ten[3][0], count);
+				if (size_h > 0)
+					mpz_tdiv_q_2exp(ten[3][1], ten[3][1], count);
+			}
+			return count;
+		}
+	}
+	Error(FATAL, E_INT, "normalizeTensor", "zero entry with non-zero size");
+	return 0;
+}
+
+int
+vectorSign(Vector v)
+{
+	int sum;
+
+	sum = mpz_sgn(v[0]) + mpz_sgn(v[1]);
+	if (sum > 0)
+		return 1;
+	else
+		if (sum < 0)
+			return -1;
+		else
+			return 0;
+}
+
+/*
+ * This function returns 1 if there are no 0 columns and all entries
+ * are >= 0, -1 if there are no 0 columns and all entries are <= 0 and
+ * 0 otherwise (ie when the signs are mixed or there are 0 columns).
+ */
+int
+matrixSign(Matrix m)
+{
+	int vecSign;
+
+	if (((vecSign = vectorSign(m[0])) != 0) && (vectorSign(m[1]) == vecSign))
+		return vecSign;
+	else
+		return 0;
+}
+
+/*
+ * This function returns 1 if there are no 0 columns and all entries
+ * are >= 0, -1 if there are no 0 columns and all entries are <= 0 and
+ * 0 otherwise (ie when the signs are mixed and/or there are 0 columns).
+ */
+int
+tensorSign(Tensor t)
+{
+	int vecSign;
+
+	if (((vecSign = vectorSign(t[0])) != 0)
+			&& (vectorSign(t[1]) == vecSign)
+			&& (vectorSign(t[2]) == vecSign)
+			&& (vectorSign(t[3]) == vecSign))
+		return vecSign;
+	else
+		return 0;
+}
+
+/*
+ * A vector is positive when at least one value is not 0 and both
+ * are greater than or equal to zero
+ */
+bool
+vectorIsPositive(Vector v)
+{
+	return vectorSign(v) == 1;
+}
+
+/*
+ * A matrix is ``positive'' when it contains no zero columns and
+ * all the entries are >= 0
+ */
+bool
+matrixIsPositive(Matrix m)
+{
+	return matrixSign(m) == 1;
+}
+
+/*
+ * A matrix is ``positive'' when it contains no zero columns and
+ * all the entries are >= 0
+ */
+bool
+tensorIsPositive(Tensor t)
+{
+	return tensorSign(t) == 1;
+}
+
+#define NEG_MPZ(x) ((x)->_mp_size = -((x)->_mp_size))
+
+void
+negateMatrix(Matrix m)
+{
+	NEG_MPZ(m[0][0]);
+	NEG_MPZ(m[0][1]);
+	NEG_MPZ(m[1][0]);
+	NEG_MPZ(m[1][1]);
+}
+
+void
+negateTensor(Tensor t)
+{
+	NEG_MPZ(t[0][0]);
+	NEG_MPZ(t[0][1]);
+	NEG_MPZ(t[1][0]);
+	NEG_MPZ(t[1][1]);
+	NEG_MPZ(t[2][0]);
+	NEG_MPZ(t[2][1]);
+	NEG_MPZ(t[3][0]);
+	NEG_MPZ(t[3][1]);
+}
+
+/*
+ * A tensor is refining when it has no zero columns
+ * (where both entries 0) and all non-zero entries have the same sign.
+ */
+bool
+tensorIsRefining(Tensor t)
+{
+	return (tensorSign(t) != 0);
+}
+
+void
+absorbSignIntoVectorPair(Vector vec0, Vector vec1, Sign sign)
+{
+	switch (sign) {
+	case SPOS :
+		break;
+	case SNEG :			/* ((c, d), (-a, -b)) */
+		mpz_neg(vec0[0], vec0[0]);
+		mpz_neg(vec0[1], vec0[1]);
+		MPZ_SWAP(vec0[0], vec1[0]);
+		MPZ_SWAP(vec0[1], vec1[1]);
+		break;
+	case SINF :			/* ((a-c,b-d), (a+c,b+d)) */
+		mpz_set(tmpa_z, vec0[0]);				/* tmp = a */
+		mpz_sub(vec0[0], tmpa_z, vec1[0]);	/* a = tmp - c */
+		mpz_add(vec1[0], tmpa_z, vec1[0]);	/* c = tmp + c */
+		mpz_set(tmpa_z, vec0[1]);				/* tmp = b */
+		mpz_sub(vec0[1], tmpa_z, vec1[1]);	/* b = tmp - d */
+		mpz_add(vec1[1], tmpa_z, vec1[1]);	/* d = tmp + d */
+		break;
+	case SZERO :		/* ((a+c,b+d), (c-a,d-b)) */
+		mpz_set(tmpa_z, vec0[0]);				/* tmp = a */
+		mpz_add(vec0[0], tmpa_z, vec1[0]);	/* a = tmp + c */
+		mpz_sub(vec1[0], vec1[0], tmpa_z);	/* c = c - tmp */
+		mpz_set(tmpa_z, vec0[1]);				/* tmp = b */
+		mpz_add(vec0[1], tmpa_z, vec1[1]);	/* b = tmp + d */
+		mpz_sub(vec1[1], vec1[1], tmpa_z);	/* d = d - tmp */
+		break;
+	default :
+		Error(FATAL, E_INT, "absorbSignIntoVectorPair", "bad sign");
+	}
+}
+
+Real
+derefToStrm(Real x)
+{
+	if (x != NULL) {
+		switch (x->gen.tag.type) {
+		case DIGSX :
+		case SIGNX :
+			break;
+		case VECTOR :
+			return derefToStrm(x->vec.strm);
+		case MATX :
+			return derefToStrm(x->matX.strm);
+		case TENXY :
+			return derefToStrm(x->tenXY.strm);
+		case ALT :
+			return derefToStrm(x->alt.redirect);
+		case CLOSURE :
+			return derefToStrm(x->cls.redirect);
+		default :
+			Error(FATAL, E_INT, "derefToStrm", "invalid real");
+		}
+	}
+	return x;
+}
+
+char *
+comparisonToString(Comparison d)
+{
+	switch (d) {
+	case LT :
+		return "lt";
+		break;
+	case GT :
+		return "gt";
+		break;
+	case EQ :
+		return "eq";
+		break;
+	default :
+	  return NULL;
+	  break;
+	}
+}
+
+char *
+signToString(Sign s)
+{
+	switch (s) {
+	case SIGN_UNKN :
+		return "unkn";
+	case SPOS :
+		return "spos";
+		break;
+	case SNEG :
+		return "sneg";
+		break;
+	case SZERO :
+		return "szer";
+		break;
+	case SINF :
+		return "sinf";
+		break;
+	default :
+	  return NULL;
+		break;
+	}
+}
+
+char *
+digitToString(Digit d)
+{
+	switch (d) {
+	case DPOS :
+		return "dpos";
+		break;
+	case DNEG :
+		return "dneg";
+		break;
+	case DZERO :
+		return "dzer";
+		break;
+	default :
+	  return NULL;
+		break;
+	}
+}
+
+char *
+typeToString(unsigned type)
+{
+	switch (type) {
+		case ALT :
+			return "alt   ";
+		case VECTOR :
+			return "vector";
+		case MATX :
+			return "matrix";
+		case TENXY :
+			return "tensor";
+		case SIGNX :
+			return "sign  ";
+		case DIGSX :
+			return "digits";
+		case CLOSURE :
+			return "closure";
+		case BOOLX :
+			return "boolx";
+		case BOOLXY :
+			return "boolxy";
+		case PREDX :
+			return "predx";
+		default :
+			Error(FATAL, E_INT, "typeToString", "bad type: %d", type);
+			return NULL;
+			break;
+	}
+}
+
+char *
+boolValToString(unsigned boolVal)
+{
+	switch (boolVal) {
+	case LAZY_TRUE :
+		return "true ";
+	case LAZY_FALSE :
+		return "false";
+	case LAZY_UNKNOWN :
+		return "unkn ";
+	default :
+	  return NULL;
+		Error(FATAL, E_INT, "boolValToString", "bad boolean value");
+	}
+}
+