_times_power.c

操作系统SunOS 4.1.3版本的源码

#ifndef lint
static char     sccsid[] = "@(#)_times_power.c 1.1 92/07/30 SMI";
#endif

/*
 * Copyright (c) 1988 by Sun Microsystems, Inc.
 */

#include "base_conversion.h"
#include <malloc.h>

void
_copy_big_float_digits(p1, p2, n)
	_BIG_FLOAT_DIGIT *p1, *p2;
	short unsigned  n;

{				/* Copies p1[n] = p2[n] */
	short unsigned  i;

	for (i = 0; i < n; i++)
		*p1++ = *p2++;
}

void
_free_big_float(p)
	_big_float     *p;

{
	/* Central routine to call free for base conversion.	 */

	char           *freearg = (char *) p;

	(void) free(freearg);
#ifdef DEBUG
	printf(" free called with %X \n", freearg);
#endif
}

void
_base_conversion_abort(ern, bcastring)
	int             ern;
	char           *bcastring;

{
	char            pstring[160];

	errno = ern;
	(void) sprintf(pstring, " libc base conversion file %s line %d: %s", __FILE__, __LINE__, bcastring);
	perror(pstring);
	abort();
}

void
_big_float_times_power(
#ifdef __STDC__
/* function to multiply a big_float times a positive power of two or ten. */
		       _big_float * pbf;	/* Operand x, to be replaced
						 * by the product x * mult **
						 * n. */
	int             mult;	/* if mult is two, x is base 10**4; if mult
				 * is ten, x is base 2**16 */
	int             n;
	int             precision;	/* Number of bits of precision
					 * ultimately required (mult=10) or
					 * number of digits of precision
					 * ultimately required (mult=2).
					 * Extra bits are allowed internally
					 * to permit correct rounding. */
/*
 * if precision is < 0, then -pfb->bexponent/{4 or 16} digits are discarded
 * on the last product.
 */
	_big_float    **pnewbf;	/* Return result *pnewbf is set to: pbf if
				 * uneventful BIG_FLOAT_TIMES_TOOBIG  if n is
				 * bigger than the tables permit ;
				 * BIG_FLOAT_TIMES_NOMEM   if pbf->blength
				 * was insufficient to hold the product, and
				 * malloc failed to produce a new block ;
				 * &newbf                  if pbf->blength
				 * was insufficient to hold the product, and
				 * a new _big_float was allocated by malloc.
				 * newbf holds the product.  It's the
				 * caller's responsibility to free this space
				 * when no longer needed. */
#else
		                       pbf, mult, n, precision, pnewbf
#endif
)
	_big_float     *pbf;
	int             mult, n, precision;
	_big_float    **pnewbf;

{
	short unsigned  base, sumlz = 0;
	unsigned short  productsize, trailing_zeros_to_delete, needed_precision, *pp, *table[3], max[3], *start[3], *lz[3], tablepower[3];
	int             i, j, itlast;
	_big_float     *pbfold = pbf;
	int             discard;

	if (precision >= 0)
		discard = -1;
	else {
		precision = -precision;
		discard = 0;
	}
	switch (mult) {
	case 2:		/* *pbf is in base 10**4 so multiply by a
				 * power of two */
		base = 10;
		max[0] = _max_tiny_powers_two;
		max[1] = _max_small_powers_two;
		max[2] = _max_big_powers_two;
		table[0] = _tiny_powers_two;
		table[1] = _small_powers_two;
		table[2] = _big_powers_two;
		lz[0] = 0;
		lz[1] = 0;
		lz[2] = 0;
		start[0] = _start_tiny_powers_two;
		start[1] = _start_small_powers_two;
		start[2] = _start_big_powers_two;
		needed_precision = 2 + (precision + 1) / 4;	/* Precision is in base
								 * ten; counts round and
								 * sticky. */
		break;
	case 10:		/* *pbf is in base 2**16 so multiply by a
				 * power of ten */
		base = 2;
		max[0] = _max_tiny_powers_ten;
		max[1] = _max_small_powers_ten;
		max[2] = _max_big_powers_ten;
		table[0] = _tiny_powers_ten;
		table[1] = _small_powers_ten;
		table[2] = _big_powers_ten;
		start[0] = _start_tiny_powers_ten;
		start[1] = _start_small_powers_ten;
		start[2] = _start_big_powers_ten;
		lz[0] = _leading_zeros_tiny_powers_ten;
		lz[1] = _leading_zeros_small_powers_ten;
		lz[2] = _leading_zeros_big_powers_ten;
		needed_precision = 2 + (precision + 1) / 16;	/* Precision is in base
								 * two; counts round and
								 * sticky. */
		break;
	}
	for (i = 0; i < 3; i++) {
		tablepower[i] = n % max[i];
		n = n / max[i];
	}
	for (itlast = 2; (itlast >= 0) && (tablepower[itlast] == 0); itlast--);
	/* Determine last i; could be 0, 1, or 2.	 */
	if (n > 0) {		/* The tables aren't big enough to accomodate
				 * mult**n, but it doesn't matter since the
				 * result would undoubtedly overflow even
				 * binary quadruple precision format.  Return
				 * an error code. */
		(void) printf("\n _times_power failed due to exponent %d %d %d leftover: %d \n", tablepower[0], tablepower[1], tablepower[2], n);
		*pnewbf = BIG_FLOAT_TIMES_TOOBIG;
		goto ret;
	}
	productsize = pbf->blength;
	for (i = 0; i < 3; i++)
		productsize += (start[i])[tablepower[i] + 1] - (start[i])[tablepower[i]];

	if (productsize < needed_precision)
		needed_precision = productsize;

	if (productsize <= pbf->bsize) {
		*pnewbf = pbf;	/* Work with *pnewbf from now on. */
	} else {		/* Need more significance than *pbf can hold. */
		char           *mallocresult;
		int             mallocarg;

		mallocarg = sizeof(_big_float) + sizeof(_BIG_FLOAT_DIGIT) * (productsize - _BIG_FLOAT_SIZE);
		mallocresult = malloc(mallocarg);
#ifdef DEBUG
		printf(" malloc arg %X result %X \n", mallocarg, (int) mallocresult);
#endif
		if (mallocresult == (char *) 0) {	/* Not enough memory
							 * left, bail out. */
			*pnewbf = BIG_FLOAT_TIMES_NOMEM;
			goto ret;
		}
		*pnewbf = (_big_float *) mallocresult;
		_copy_big_float_digits((*pnewbf)->bsignificand, pbf->bsignificand, pbf->blength);
		(*pnewbf)->blength = pbf->blength;
		(*pnewbf)->bexponent = pbf->bexponent;
		pbf = *pnewbf;
		pbf->bsize = productsize;
	}

	/* pbf now points to the input and the output big_floats.	 */

	for (i = 0; i <= itlast; i++)
		if (tablepower[i] != 0) {	/* Step through each of the
						 * tables. */
			unsigned        lengthx, lengthp;

			/* Powers of 10**4 have leading zeros in base 2**16. */
			lengthp = (start[i])[tablepower[i] + 1] - (start[i])[tablepower[i]];
			lengthx = pbf->blength;

			if (discard >= 0)
				switch (base) {
				case 2:
					discard = (-pbf->bexponent) / 16;
					break;
				case 10:
					discard = (-pbf->bexponent) / 4;
					break;
				}

#ifdef DEBUG
			{
				long            basexp;
				int             id;

				printf(" step %d x operand length %d \n", i, lengthx);
				_display_big_float(pbf, base);
				printf(" step %d p operand length %d power %d \n", i, lengthp, tablepower[i]);
				basexp = (base == 2) ? (lz[i])[tablepower[i]] : 0;
				for (id = 0; id < lengthp; id++) {
					printf("+ %d * ", (table[i])[id + (start[i])[tablepower[i]]]);
					if (base == 2)
						printf("2**%d", 16 * (basexp + id));
					if (base == 10)
						printf("10**%d", 4 * (basexp + id));
					if ((id % 4) == 3)
						printf("\n");
				}
				printf("\n");
			}
			if ((i == itlast) && (discard >= 0))
				printf(" alternative discard %d digits \n", discard);
#endif

			if (base == 2) {
				sumlz += (lz[i])[tablepower[i]];
				pbf->bexponent += 16 * (lz[i])[tablepower[i]];
			}
			if (lengthp == 1) {	/* Special case - multiply by
						 * <= 10**4 or 2**13 */
				switch (base) {
				case 10:
					_multiply_base_ten_by_two(pbf, tablepower[i]);
					break;
				case 2:
					_multiply_base_two(pbf, (_BIG_FLOAT_DIGIT) ((table[i])[tablepower[i]]), (unsigned long) 0);
					break;
				}
#ifdef DEBUG
				assert(pbf->blength <= pbf->bsize);
#endif
			} else if (lengthx == 1) {	/* Special case of short
							 * multiplicand. */
				_BIG_FLOAT_DIGIT multiplier = pbf->bsignificand[0];

				_copy_big_float_digits(pbf->bsignificand, (unsigned short *) &((table[i])[(start[i])[tablepower[i]]]), lengthp);
				pbf->blength = lengthp;
				switch (base) {
				case 10:
					_multiply_base_ten(pbf, multiplier);
					break;
				case 2:
					_multiply_base_two(pbf, multiplier, (unsigned long) 0);
					break;
				}
#ifdef DEBUG
				assert(pbf->blength <= pbf->bsize);
#endif
			} else {/* General case. */
				short unsigned  canquit;
				short unsigned  excess;

				/*
				 * The result will be accumulated in *pbf
				 * from most significant to least
				 * significant.
				 */

				/* Generate criterion for early termination.	 */
				switch (base) {
				case 2:
					canquit = 65536;
					break;
				case 10:
					canquit = 10000;
					break;
				}
				canquit -= 3 + ((lengthx < lengthp) ? lengthx : lengthp);

				pbf->bsignificand[lengthx + lengthp - 1] = 0;	/* Only gets filled by
										 * carries. */
				for (j = lengthx + lengthp - 2; j >= 0; j--) {
					int             istart = j - lengthp + 1, istop = lengthx - 1;
					short unsigned  lengthprod;
					_BIG_FLOAT_DIGIT product[3];

					pp = (unsigned short *) &((table[i])[(start[i])[tablepower[i]]]);
					if (j < istop)
						istop = j;
					if (0 > istart)
						istart = 0;

					switch (base) {
					case 2:
						_multiply_base_two_vector((short unsigned) (istop - istart + 1), &(pbf->bsignificand[istart]), &(pp[j - istop]), product);
						if (product[2] != 0)
							_carry_propagate_two((unsigned long) product[2], &(pbf->bsignificand[j + 2]));
						if (product[1] != 0)
							_carry_propagate_two((unsigned long) product[1], &(pbf->bsignificand[j + 1]));
						break;
					case 10:
						_multiply_base_ten_vector((short unsigned) (istop - istart + 1), &(pbf->bsignificand[istart]), &(pp[j - istop]), product);
						if (product[2] != 0)
							_carry_propagate_ten((unsigned long) product[2], &(pbf->bsignificand[j + 2]));
						if (product[1] != 0)
							_carry_propagate_ten((unsigned long) product[1], &(pbf->bsignificand[j + 1]));
						break;
					}
					pbf->bsignificand[j] = product[0];
					lengthprod = lengthx + lengthp;
					if (pbf->bsignificand[lengthprod - 1] == 0)
						lengthprod--;
					if (lengthprod > needed_precision)
						excess = lengthprod - needed_precision;
					else
						excess = 0;
					if ((i == itlast) && ((j + 2) <= excess) && (pbf->bsignificand[j + 1] <= canquit)
					    && ((pbf->bsignificand[j + 1] | pbf->bsignificand[j]) != 0)) {
						/*
						 * On the last
						 * multiplication, it's not
						 * necessary to develop the
						 * entire product, if further
						 * digits can't possibly
						 * affect significant digits,
						 * unless there's a chance
						 * the product might be
						 * exact!
						 */
						/*
						 * Note that the product
						 * might be exact if the j
						 * and j+1 terms are zero; if
						 * they are non-zero, then it
						 * won't be after they're
						 * discarded.
						 */

						excess = j + 2;	/* Can discard j+1, j,
								 * ... 0 */
#ifdef DEBUG
						printf(" decided to quit early at j %d since s[j+1] is %d <= %d \n", j, pbf->bsignificand[j + 1], canquit);
						printf(" s[j+2..j] are %d %d %d \n", pbf->bsignificand[j + 2], pbf->bsignificand[j + 1], pbf->bsignificand[j]);
						printf(" requested precision %d needed_precision %d big digits out of %d \n", precision, needed_precision, lengthprod);
#endif
						if ((discard >= 0) && ((j + 2) > (discard - (int) sumlz))) {
#ifdef DEBUG
							printf(" early quit rejected because j+2 = %d > %d = discard \n", j + 2, discard);
#endif
							goto pastdiscard;
						}
						pbf->bsignificand[excess] |= 1;	/* Sticky bit on. */
#ifdef DEBUG
						printf(" discard %d digits - last gets %d \n", excess, pbf->bsignificand[excess]);
#endif
						trailing_zeros_to_delete = excess;
						goto donegeneral;
					}
			pastdiscard:	;
#ifdef DEBUG
					/*
					 * else { printf(" early termination
					 * rejected at j %d since s[j+1] =
					 * %d, canquit = %d \n", j,
					 * pbf->bsignificand[j + 1],
					 * canquit); printf(" s[j+2..j] are
					 * %d %d %d \n", pbf->bsignificand[j
					 * + 2], pbf->bsignificand[j + 1],
					 * pbf->bsignificand[j]); printf("
					 * requested precision %d
					 * needed_precision %d big digits out
					 * of %d \n", precision,
					 * needed_precision, lengthprod); }
					 */
#endif
				}
				trailing_zeros_to_delete = 0;
		donegeneral:
				pbf->blength = lengthx + lengthp;
				if (pbf->bsignificand[pbf->blength - 1] == 0)
					pbf->blength--;
				for (; pbf->bsignificand[trailing_zeros_to_delete] == 0; trailing_zeros_to_delete++);
				/*
				 * Look for additional trailing zeros to
				 * delete.
				 */
				if (trailing_zeros_to_delete != 0) {
#ifdef DEBUG
					printf(" %d trailing zeros deleted \n", trailing_zeros_to_delete);
#endif
					_copy_big_float_digits(pbf->bsignificand, &(pbf->bsignificand[trailing_zeros_to_delete]), pbf->blength - trailing_zeros_to_delete);
					pbf->blength -= trailing_zeros_to_delete;
					switch (base) {
					case 2:
						pbf->bexponent += 16 * trailing_zeros_to_delete;
						break;
					case 10:
						pbf->bexponent += 4 * trailing_zeros_to_delete;
						break;
					}
				}
			}
		}
	if ((pbfold != pbf) && (pbf->blength <= pbfold->bsize)) {	/* Don't need that huge
									 * buffer after all! */
#ifdef DEBUG
		printf(" free called from times_power because final length %d <= %d original size \n", pbf->blength, pbfold->bsize);
#endif

		/* Copy product to original buffer. */
		pbfold->blength = pbf->blength;
		pbfold->bexponent = pbf->bexponent;
		_copy_big_float_digits(pbfold->bsignificand, pbf->bsignificand, pbf->blength);
		_free_big_float(*pnewbf);	/* Free new buffer. */
		*pnewbf = pbfold;	/* New buffer pointer now agrees with
					 * original. */
	}
ret:
	return;
}