_times_power.c revision 5d54f3d8999eac1762fe0a8c7177d20f1f201fae
/*
* CDDL HEADER START
*
* The contents of this file are subject to the terms of the
* Common Development and Distribution License, Version 1.0 only
* (the "License"). You may not use this file except in compliance
* with the License.
*
* You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
* See the License for the specific language governing permissions
* and limitations under the License.
*
* When distributing Covered Code, include this CDDL HEADER in each
* file and include the License file at usr/src/OPENSOLARIS.LICENSE.
* If applicable, add the following below this CDDL HEADER, with the
* fields enclosed by brackets "[]" replaced with your own identifying
* information: Portions Copyright [yyyy] [name of copyright owner]
*
* CDDL HEADER END
*/
/*
* Copyright 1995 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
#pragma ident "%Z%%M% %I% %E% SMI"
#include "base_conversion.h"
#include <malloc.h>
void
short unsigned n)
{ /* Copies p1[n] = p2[n] */
short unsigned i;
for (i = 0; i < n; i++)
}
void
{
/* Central routine to call free for base conversion. */
char *freearg = (char *) p;
#ifdef DEBUG
#endif
}
void
{
char pstring[160];
(void) sprintf(pstring, " libc base conversion file %s line %d: %s", __FILE__, __LINE__, bcastring);
abort();
}
/*
* function to multiply a big_float times a positive power of two or ten.
*
* Arguments
* pbf: Operand x, to be replaced by the product x * mult ** n.
* mult: if mult is two, x is base 10**4;
* if mult is ten, x is base 2**16
* n:
* 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.
* 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.
*
* if precision is < 0, then -pfb->bexponent/{4 or 16} digits are discarded
* on the last product.
*/
void
_big_float **pnewbf)
{
unsigned short productsize, trailing_zeros_to_delete, needed_precision, *pp, *table[3], max[3], *start[3], *lz[3], tablepower[3];
int i, j, itlast;
int discard;
if (precision >= 0)
discard = -1;
else {
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;
table[0] = _tiny_powers_two;
lz[0] = 0;
lz[1] = 0;
lz[2] = 0;
start[0] = _start_tiny_powers_two;
* 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;
table[0] = _tiny_powers_ten;
start[0] = _start_tiny_powers_ten;
* two; counts round and
* sticky. */
break;
}
for (i = 0; i < 3; i++) {
tablepower[i] = n % max[i];
n = n / max[i];
}
/* 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);
goto ret;
}
for (i = 0; i < 3; i++)
if (productsize < needed_precision)
} else { /* Need more significance than *pbf can hold. */
char *mallocresult;
int mallocarg;
#ifdef DEBUG
#endif
if (mallocresult == (char *) 0) { /* Not enough memory
* left, bail out. */
goto ret;
}
}
/* 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. */
/* Powers of 10**4 have leading zeros in base 2**16. */
if (discard >= 0)
switch (base) {
case 2:
break;
case 10:
break;
}
#ifdef DEBUG
{
long basexp;
int id;
if (base == 2)
if (base == 10)
printf("\n");
}
printf("\n");
}
#endif
if (base == 2) {
}
* <= 10**4 or 2**13 */
switch (base) {
case 10:
break;
case 2:
break;
}
#ifdef DEBUG
#endif
* multiplicand. */
_copy_big_float_digits(pbf->bsignificand, (unsigned short *) &((table[i])[(start[i])[tablepower[i]]]), lengthp);
switch (base) {
case 10:
break;
case 2:
break;
}
#ifdef DEBUG
#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 = (short unsigned)65536;
break;
case 10:
canquit = 10000;
break;
}
* carries. */
short unsigned lengthprod;
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)
if (product[1] != 0)
break;
case 10:
_multiply_base_ten_vector((short unsigned) (istop - istart + 1), &(pbf->bsignificand[istart]), &(pp[j - istop]), product);
if (product[2] != 0)
if (product[1] != 0)
break;
}
lengthprod--;
if (lengthprod > needed_precision)
else
excess = 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.
*/
* ... 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
#ifdef DEBUG
#endif
goto pastdiscard;
}
#ifdef DEBUG
#endif
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
}
/*
* Look for additional trailing zeros to
* delete.
*/
/*
* fix for bug 1070565; if too many trailing
* zeroes are deleted, we'll violate the
* assertion that bexponent is in [-3,+4]
*/
if (base == 10) {
if ((int)trailing_zeros_to_delete > deletelimit) {
#ifdef DEBUG
printf("\n __x_power trailing zeros delete count lowered from %d to
#endif
}
}
if (trailing_zeros_to_delete != 0) {
#ifdef DEBUG
#endif
_copy_big_float_digits(pbf->bsignificand, &(pbf->bsignificand[trailing_zeros_to_delete]), pbf->blength - trailing_zeros_to_delete);
switch (base) {
case 2:
break;
case 10:
break;
}
}
}
}
* 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. */
* original. */
}
ret:
return;
}