zmul.c

Calc Software Package for Number Calc

		*h1++ = sival.silow;
		*h2++ = sival.silow;
		carryACBD = sival.sihigh;
	}

	/*
	 * Now add in the bottom half of B*D and the top half of A*C.
	 * These additions are straightforward, except that A*C should
	 * be done first because of possible carries from B*D, and the
	 * top half of A*C might not exist.  Add in one of the carries
	 * from the previous addition while we are at it.
	 */
	h1 = baseAC + shift;
	hd = baseAC;
	carry = carryACBD;
	len = sizetotal - 3 * shift;
	while (len--) {
		sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
		*hd++ = sival.silow;
		carry = sival.sihigh;
	}
	while (carry) {
		sival.ivalue = ((FULL) *hd) + carry;
		*hd++ = sival.silow;
		carry = sival.sihigh;
	}

	h1 = baseBD;
	hd = baseBD + shift;
	carry = 0;
	len = shift;
	while (len--) {
		sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
		*hd++ = sival.silow;
		carry = sival.sihigh;
	}
	while (carry) {
		sival.ivalue = ((FULL) *hd) + carry;
		*hd++ = sival.silow;
		carry = sival.sihigh;
	}

	/*
	 * Now finally add in the other delayed carry from the
	 * above addition.
	 */
	hd = baseAC + shift;
	while (carryACBD) {
		sival.ivalue = ((FULL) *hd) + carryACBD;
		*hd++ = sival.silow;
		carryACBD = sival.sihigh;
	}

	/*
	 * Now finally add or subtract (A-B)*(D-C) into the final result at
	 * the correct position (S), according to whether it is positive or
	 * negative.  When subtracting, the answer cannot go negative.
	 */
	h1 = baseABDC;
	hd = ans + shift;
	carry = 0;
	len = sizeABDC;
	if (neg) {
		while (len--) {
			sival.ivalue = BASE1 - ((FULL) *hd) +
				((FULL) *h1++) + carry;
			*hd++ = (HALF)(BASE1 - sival.silow);
			carry = sival.sihigh;
		}
		while (carry) {
			sival.ivalue = BASE1 - ((FULL) *hd) + carry;
			*hd++ = (HALF)(BASE1 - sival.silow);
			carry = sival.sihigh;
		}
	} else {
		while (len--) {
			sival.ivalue = ((FULL) *h1++) + ((FULL) *hd) + carry;
			*hd++ = sival.silow;
			carry = sival.sihigh;
		}
		while (carry) {
			sival.ivalue = ((FULL) *hd) + carry;
			*hd++ = sival.silow;
			carry = sival.sihigh;
		}
	}

	/*
	 * Finally determine the size of the final result and return that.
	 */
	len = sizetotal;
	hd = &ans[len - 1];
	while ((*hd == 0) && (len > 1)) {
		hd--;
		len--;
	}
	tempbuf = temp;
	return len;
}


/*
 * Square a number by using the following formula recursively:
 *	(A*S+B)^2 = (S^2+S)*A^2 + (S+1)*B^2 - S*(A-B)^2
 * where S is a power of 2^16, and so multiplies by it are shifts,
 * and A and B are the left and right halfs of the number to square.
 */
void
zsquare(ZVALUE z, ZVALUE *res)
{
	LEN len;

	if (ziszero(z)) {
		*res = _zero_;
		return;
	}
	if (zisunit(z)) {
		*res = _one_;
		return;
	}

	/*
	 * Allocate a temporary array if necessary for the recursion to use.
	 * The array needs to be allocated large enough for all of the
	 * temporary results to fit in.	 This size is about 3 times the
	 * size of the original number, since each recursion level uses 3/2
	 * of the size of its given number, and whose size is 1/2 the size
	 * of the previous level.  The sum of the infinite series is 3.
	 * Allocate some extra words for rounding up the sizes.
	 */
	len = 3 * z.len + 32;
	tempbuf = zalloctemp(len);

	res->sign = 0;
	res->v = alloc((z.len+2) * 2);
	/*
	 * Without the memset below, Purify reports that dosquare()
	 *	 will read uninitialized memory at the dosquare() line below
	 *	 the comment:
	 *
	 *		uninitialized memory read (see zsquare)
	 *
	 * This problem occurs during regression test #622 and may
	 * be duplicated by executing:
	 *
	 *	config("sq2", 2);
	 *	0xffff0000ffffffff00000000ffff0000000000000000ffff^2;
	 */
	memset((char *)res->v, 0, ((z.len+2) * 2)*sizeof(HALF));
	res->len = dosquare(z.v, z.len, res->v);
}


/*
 * Recursive routine to square a number by splitting it up into two numbers
 * of half the size, and using the results of squaring the subpieces.
 * The result is placed in the indicated array, which must be large
 * enough for the result (size * 2).  Returns the size of the result.
 * This uses a temporary array which must be 3 times as large as the
 * size of the number at the top level recursive call.
 *
 * given:
 *	vp		value to be squared
 *	size		length of value to square
 *	ans		location for result
 */
static LEN
dosquare(HALF *vp, LEN size, HALF *ans)
{
	LEN shift;		/* amount numbers are shifted by */
	LEN sizeA;		/* size of left half of number to square */
	LEN sizeB;		/* size of right half of number to square */
	LEN sizeAA;		/* size of square of left half */
	LEN sizeBB;		/* size of square of right half */
	LEN sizeAABB;		/* size of sum of squares of A and B */
	LEN sizeAB;		/* size of difference of A and B */
	LEN sizeABAB;		/* size of square of difference of A and B */
	LEN subsize;		/* size of difference of halfs */
	LEN copysize;		/* size of number left to copy */
	LEN sumsize;		/* size of sum */
	LEN sizetotal;		/* total size of square */
	LEN len;		/* temporary length */
	LEN len1;		/* another temporary length */
	FULL carry;		/* carry digit for small multiplications */
	FULL digit;		/* single digit multiplying by */
	HALF *temp;		/* base for temporary calculations */
	HALF *baseA;		/* base of left half of number */
	HALF *baseB;		/* base of right half of number */
	HALF *baseAA;		/* base of square of left half of number */
	HALF *baseBB;		/* base of square of right half of number */
	HALF *baseAABB;		/* base of sum of squares of A and B */
	HALF *baseAB;		/* base of difference of A and B */
	HALF *baseABAB;		/* base of square of difference of A and B */
	register HALF *hd, *h1, *h2, *h3;	/* for inner loops */
	SIUNION sival;		/* for addition of digits */

	/*
	 * First trim the number of leading zeroes.
	 */
	hd = &vp[size - 1];
	while ((*hd == 0) && (size > 1)) {
		size--;
		hd--;
	}
	sizetotal = size + size;

	/*
	 * If the number has only a small number of digits, then use the
	 * (almost) normal multiplication method.  Multiply each halfword
	 * only by those halfwards further on in the number.  Missed terms
	 * will then be the same pairs of products repeated, and the squares
	 * of each halfword.  The first case is handled by doubling the
	 * result.  The second case is handled explicitly.  The number of
	 * digits where the algorithm changes is settable from 2 to maxint.
	 */
	if (size < conf->sq2) {
		hd = ans;
		len = sizetotal;
		while (len--)
			*hd++ = 0;

		h2 = vp;
		hd = ans + 1;
		for (len = size; len--; hd += 2) {
			digit = (FULL) *h2++;
			if (digit == 0)
				continue;
			h3 = h2;
			h1 = hd;
			carry = 0;
			len1 = len;
			while (len1 >= 4) {	/* expand the loop some */
				len1 -= 4;
				sival.ivalue = (digit * ((FULL) *h3++))
					+ ((FULL) *h1) + carry;
				*h1++ = sival.silow;
				sival.ivalue = (digit * ((FULL) *h3++))
					+ ((FULL) *h1) + ((FULL) sival.sihigh);
				*h1++ = sival.silow;
				sival.ivalue = (digit * ((FULL) *h3++))
					+ ((FULL) *h1) + ((FULL) sival.sihigh);
				*h1++ = sival.silow;
				sival.ivalue = (digit * ((FULL) *h3++))
					+ ((FULL) *h1) + ((FULL) sival.sihigh);
				*h1++ = sival.silow;
				carry = sival.sihigh;
			}
			while (len1--) {
				sival.ivalue = (digit * ((FULL) *h3++))
					+ ((FULL) *h1) + carry;
				*h1++ = sival.silow;
				carry = sival.sihigh;
			}
			while (carry) {
				sival.ivalue = ((FULL) *h1) + carry;
				*h1++ = sival.silow;
				carry = sival.sihigh;
			}
		}

		/*
		 * Now double the result.
		 * There is no final carry to worry about because we
		 * handle all digits of the result which must fit.
		 */
		carry = 0;
		hd = ans;
		len = sizetotal;
		while (len--) {
			digit = ((FULL) *hd);
			sival.ivalue = digit + digit + carry;
			/* ignore Saber-C warning #112 - get ushort from uint */
			/*	  ok to ignore on name dosquare`sival */
			*hd++ = sival.silow;
			carry = sival.sihigh;
		}

		/*
		 * Now add in the squares of each halfword.
		 */
		carry = 0;
		hd = ans;
		h3 = vp;
		len = size;
		while (len--) {
			digit = ((FULL) *h3++);
			sival.ivalue = digit * digit + ((FULL) *hd) + carry;
			*hd++ = sival.silow;
			carry = sival.sihigh;
			sival.ivalue = ((FULL) *hd) + carry;
			*hd++ = sival.silow;
			carry = sival.sihigh;
		}
		while (carry) {
			sival.ivalue = ((FULL) *hd) + carry;
			*hd++ = sival.silow;
			carry = sival.sihigh;
		}

		/*
		 * Finally return the size of the result.
		 */
		len = sizetotal;
		hd = &ans[len - 1];
		while ((*hd == 0) && (len > 1)) {
			len--;
			hd--;
		}
		return len;
	}

	/*
	 * The number to be squared is large.
	 * Allocate temporary space and determine the sizes and
	 * positions of the values to be calculated.
	 */
	temp = tempbuf;
	tempbuf += (3 * (size + 1) / 2);

	sizeA = size / 2;
	sizeB = size - sizeA;
	shift = sizeB;
	baseA = vp + sizeB;
	baseB = vp;
	baseAA = &ans[shift * 2];
	baseBB = ans;
	baseAABB = temp;
	baseAB = temp;
	baseABAB = &temp[shift];

	/*
	 * Trim the second number of leading zeroes.
	 */
	hd = &baseB[sizeB - 1];
	while ((*hd == 0) && (sizeB > 1)) {
		sizeB--;
		hd--;
	}

	/*
	 * Now to proceed to calculate the result using the formula.
	 *	(A*S+B)^2 = (S^2+S)*A^2 + (S+1)*B^2 - S*(A-B)^2.
	 * The steps are done in the following order:
	 *	S^2*A^2 + B^2
	 *	A^2 + B^2
	 *	(S^2+S)*A^2 + (S+1)*B^2
	 *	(A-B)^2
	 *	(S^2+S)*A^2 + (S+1)*B^2 - S*(A-B)^2.
	 *
	 * Begin by forming the squares of two the halfs concatenated
	 * together in the final result location.  Make sure that the
	 * highest words of the results are zero.
	 */
	sizeBB = dosquare(baseB, sizeB, baseBB);
	hd = &baseBB[sizeBB];
	len = shift * 2 - sizeBB;
	while (len--)
		*hd++ = 0;

	sizeAA = dosquare(baseA, sizeA, baseAA);
	hd = &baseAA[sizeAA];
	len = sizetotal - shift * 2 - sizeAA;
	while (len--)
		*hd++ = 0;

	/*
	 * Sum the two squares into a temporary location.
	 */
	if (sizeAA >= sizeBB) {
		h1 = baseAA;
		h2 = baseBB;
		sizeAABB = sizeAA;
		sumsize = sizeBB;
	} else {
		h1 = baseBB;
		h2 = baseAA;
		sizeAABB = sizeBB;
		sumsize = sizeAA;
	}
	copysize = sizeAABB - sumsize;

	hd = baseAABB;
	carry = 0;
	while (sumsize--) {
		sival.ivalue = ((FULL) *h1++) + ((FULL) *h2++) + carry;
		*hd++ = sival.silow;
		carry = sival.sihigh;
	}
	while (copysize--) {
		sival.ivalue = ((FULL) *h1++) + carry;
		*hd++ = sival.silow;
		carry = sival.sihigh;
	}
	if (carry) {
		*hd = (HALF)carry;
		sizeAABB++;
	}

	/*
	 * Add the sum back into the previously calculated squares
	 * shifted over to the proper location.
	 */
	h1 = baseAABB;
	hd = ans + shift;
	carry = 0;
	len = sizeAABB;
	while (len--) {
		sival.ivalue = ((FULL) *hd) + ((FULL) *h1++) + carry;
		*hd++ = sival.silow;
		carry = sival.sihigh;
	}
	while (carry) {
		/* uninitialized memory read (see zsquare) */
		sival.ivalue = ((FULL) *hd) + carry;
		*hd++ = sival.silow;
		carry = sival.sihigh;
	}

	/*
	 * Calculate the absolute value of the difference of the two halfs
	 * into a temporary location.
	 */
	if (sizeA == sizeB) {
		len = sizeA;
		h1 = &baseA[len - 1];
		h2 = &baseB[len - 1];
		while ((len > 1) && (*h1 == *h2)) {
			len--;
			h1--;
			h2--;
		}
	}
	if ((sizeA > sizeB) || ((sizeA == sizeB) && (*h1 > *h2))) {
		h1 = baseA;
		h2 = baseB;
		sizeAB = sizeA;
		subsize = sizeB;
	} else {
		h1 = baseB;
		h2 = baseA;
		sizeAB = sizeB;
		subsize = sizeA;
	}
	copysize = sizeAB - subsize;

	hd = baseAB;
	carry = 0;
	while (subsize--) {
		sival.ivalue = BASE1 - ((FULL) *h1++) + ((FULL) *h2++) + carry;
		*hd++ = (HALF)(BASE1 - sival.silow);
		carry = sival.sihigh;
	}
	while (copysize--) {
		sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
		*hd++ = (HALF)(BASE1 - sival.silow);
		carry = sival.sihigh;
	}

	hd = &baseAB[sizeAB - 1];
	while ((*hd == 0) && (sizeAB > 1)) {
		sizeAB--;
		hd--;
	}

	/*
	 * Now square the number into another temporary location,
	 * and subtract that from the final result.
	 */
	sizeABAB = dosquare(baseAB, sizeAB, baseABAB);

	h1 = baseABAB;
	hd = ans + shift;
	carry = 0;
	while (sizeABAB--) {
		sival.ivalue = BASE1 - ((FULL) *hd) + ((FULL) *h1++) + carry;
		*hd++ = (HALF)(BASE1 - sival.silow);
		carry = sival.sihigh;
	}
	while (carry) {
		sival.ivalue = BASE1 - ((FULL) *hd) + carry;
		*hd++ = (HALF)(BASE1 - sival.silow);
		carry = sival.sihigh;
	}

	/*
	 * Return the size of the result.
	 */
	len = sizetotal;
	hd = &ans[len - 1];
	while ((*hd == 0) && (len > 1)) {
		len--;
		hd--;
	}
	tempbuf = temp;
	return len;
}


/*
 * Return a pointer to a buffer to be used for holding a temporary number.
 * The buffer will be at least as large as the specified number of HALFs,
 * and remains valid until the next call to this routine.  The buffer cannot
 * be freed by the caller.  There is only one temporary buffer, and so as to
 * avoid possible conflicts this is only used by the lowest level routines
 * such as divide, multiply, and square.
 *
 * given:
 *	len		required number of HALFs in buffer
 */
HALF *
zalloctemp(LEN len)
{
	HALF *hp;
	static LEN buflen;	/* current length of temp buffer */
	static HALF *bufptr;	/* pointer to current temp buffer */

	if (len <= buflen)
		return bufptr;

	/*
	 * We need to grow the temporary buffer.
	 * First free any existing buffer, and then allocate the new one.
	 * While we are at it, make the new buffer bigger than necessary
	 * in order to reduce the number of reallocations.
	 */
	len += 100;
	if (buflen) {
		buflen = 0;
		free(bufptr);
	}
	/* don't call alloc() because _math_abort_ may not be set right */
	hp = (HALF *) malloc((len+1) * sizeof(HALF));
	if (hp == NULL) {
		math_error("No memory for temp buffer");
		/*NOTREACHED*/
	}
	bufptr = hp;
	buflen = len;
	return hp;
}