1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License, Version 1.0 only 6 * (the "License"). You may not use this file except in compliance 7 * with the License. 8 * 9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 10 * or http://www.opensolaris.org/os/licensing. 11 * See the License for the specific language governing permissions 12 * and limitations under the License. 13 * 14 * When distributing Covered Code, include this CDDL HEADER in each 15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 16 * If applicable, add the following below this CDDL HEADER, with the 17 * fields enclosed by brackets "[]" replaced with your own identifying 18 * information: Portions Copyright [yyyy] [name of copyright owner] 19 * 20 * CDDL HEADER END 21 */ 22 /* 23 * Copyright 1995 Sun Microsystems, Inc. All rights reserved. 24 * Use is subject to license terms. 25 */ 26 27 #include "base_conversion.h" 28 #include <malloc.h> 29 30 void 31 _copy_big_float_digits(_BIG_FLOAT_DIGIT *p1, _BIG_FLOAT_DIGIT *p2, 32 short unsigned n) 33 { /* Copies p1[n] = p2[n] */ 34 short unsigned i; 35 36 for (i = 0; i < n; i++) 37 *p1++ = *p2++; 38 } 39 40 void 41 _free_big_float(_big_float *p) 42 { 43 /* Central routine to call free for base conversion. */ 44 45 char *freearg = (char *) p; 46 47 (void) free(freearg); 48 #ifdef DEBUG 49 printf(" free called with %X \n", freearg); 50 #endif 51 } 52 53 void 54 _base_conversion_abort(int ern, char *bcastring) 55 { 56 char pstring[160]; 57 58 errno = ern; 59 (void) sprintf(pstring, " libc base conversion file %s line %d: %s", __FILE__, __LINE__, bcastring); 60 perror(pstring); 61 abort(); 62 } 63 64 /* 65 * function to multiply a big_float times a positive power of two or ten. 66 * 67 * Arguments 68 * pbf: Operand x, to be replaced by the product x * mult ** n. 69 * mult: if mult is two, x is base 10**4; 70 * if mult is ten, x is base 2**16 71 * n: 72 * precision: Number of bits of precision ultimately required 73 * (mult=10) or number of digits of precision ultimately 74 * required (mult=2). 75 * Extra bits are allowed internally to permit correct 76 * rounding. 77 * pnewbf: Return result *pnewbf is set to: pbf if uneventful 78 * BIG_FLOAT_TIMES_TOOBIG if n is bigger than the tables 79 * permit ; 80 * BIG_FLOAT_TIMES_NOMEM if pbf->blength was 81 * insufficient to hold the product, and malloc failed to 82 * produce a new block ; 83 * &newbf if pbf->blength was 84 * insufficient to hold the product, and a new _big_float 85 * was allocated by malloc. newbf holds the product. 86 * It's the caller's responsibility to free this space 87 * when no longer needed. 88 * 89 * if precision is < 0, then -pfb->bexponent/{4 or 16} digits are discarded 90 * on the last product. 91 */ 92 void 93 _big_float_times_power(_big_float *pbf, int mult, int n, int precision, 94 _big_float **pnewbf) 95 { 96 short unsigned base, sumlz = 0; 97 unsigned short productsize, trailing_zeros_to_delete, needed_precision, *pp, *table[3], max[3], *start[3], *lz[3], tablepower[3]; 98 int i, j, itlast; 99 _big_float *pbfold = pbf; 100 int discard; 101 102 if (precision >= 0) 103 discard = -1; 104 else { 105 precision = -precision; 106 discard = 0; 107 } 108 switch (mult) { 109 case 2: /* *pbf is in base 10**4 so multiply by a 110 * power of two */ 111 base = 10; 112 max[0] = _max_tiny_powers_two; 113 max[1] = _max_small_powers_two; 114 max[2] = _max_big_powers_two; 115 table[0] = _tiny_powers_two; 116 table[1] = _small_powers_two; 117 table[2] = _big_powers_two; 118 lz[0] = 0; 119 lz[1] = 0; 120 lz[2] = 0; 121 start[0] = _start_tiny_powers_two; 122 start[1] = _start_small_powers_two; 123 start[2] = _start_big_powers_two; 124 needed_precision = 2 + (precision + 1) / 4; /* Precision is in base 125 * ten; counts round and 126 * sticky. */ 127 break; 128 case 10: /* *pbf is in base 2**16 so multiply by a 129 * power of ten */ 130 base = 2; 131 max[0] = _max_tiny_powers_ten; 132 max[1] = _max_small_powers_ten; 133 max[2] = _max_big_powers_ten; 134 table[0] = _tiny_powers_ten; 135 table[1] = _small_powers_ten; 136 table[2] = _big_powers_ten; 137 start[0] = _start_tiny_powers_ten; 138 start[1] = _start_small_powers_ten; 139 start[2] = _start_big_powers_ten; 140 lz[0] = _leading_zeros_tiny_powers_ten; 141 lz[1] = _leading_zeros_small_powers_ten; 142 lz[2] = _leading_zeros_big_powers_ten; 143 needed_precision = 2 + (precision + 1) / 16; /* Precision is in base 144 * two; counts round and 145 * sticky. */ 146 break; 147 } 148 for (i = 0; i < 3; i++) { 149 tablepower[i] = n % max[i]; 150 n = n / max[i]; 151 } 152 for (itlast = 2; (itlast >= 0) && (tablepower[itlast] == 0); itlast--); 153 /* Determine last i; could be 0, 1, or 2. */ 154 if (n > 0) { /* The tables aren't big enough to accomodate 155 * mult**n, but it doesn't matter since the 156 * result would undoubtedly overflow even 157 * binary quadruple precision format. Return 158 * an error code. */ 159 (void) printf("\n _times_power failed due to exponent %d %d %d leftover: %d \n", tablepower[0], tablepower[1], tablepower[2], n); 160 *pnewbf = BIG_FLOAT_TIMES_TOOBIG; 161 goto ret; 162 } 163 productsize = pbf->blength; 164 for (i = 0; i < 3; i++) 165 productsize += (start[i])[tablepower[i] + 1] - (start[i])[tablepower[i]]; 166 167 if (productsize < needed_precision) 168 needed_precision = productsize; 169 170 if (productsize <= pbf->bsize) { 171 *pnewbf = pbf; /* Work with *pnewbf from now on. */ 172 } else { /* Need more significance than *pbf can hold. */ 173 char *mallocresult; 174 int mallocarg; 175 176 mallocarg = sizeof(_big_float) + sizeof(_BIG_FLOAT_DIGIT) * (productsize - _BIG_FLOAT_SIZE); 177 mallocresult = malloc(mallocarg); 178 #ifdef DEBUG 179 printf(" malloc arg %X result %X \n", mallocarg, (int) mallocresult); 180 #endif 181 if (mallocresult == (char *) 0) { /* Not enough memory 182 * left, bail out. */ 183 *pnewbf = BIG_FLOAT_TIMES_NOMEM; 184 goto ret; 185 } 186 *pnewbf = (_big_float *) mallocresult; 187 _copy_big_float_digits((*pnewbf)->bsignificand, pbf->bsignificand, pbf->blength); 188 (*pnewbf)->blength = pbf->blength; 189 (*pnewbf)->bexponent = pbf->bexponent; 190 pbf = *pnewbf; 191 pbf->bsize = productsize; 192 } 193 194 /* pbf now points to the input and the output big_floats. */ 195 196 for (i = 0; i <= itlast; i++) 197 if (tablepower[i] != 0) { /* Step through each of the 198 * tables. */ 199 unsigned lengthx, lengthp; 200 201 /* Powers of 10**4 have leading zeros in base 2**16. */ 202 lengthp = (start[i])[tablepower[i] + 1] - (start[i])[tablepower[i]]; 203 lengthx = pbf->blength; 204 205 if (discard >= 0) 206 switch (base) { 207 case 2: 208 discard = (-pbf->bexponent) / 16; 209 break; 210 case 10: 211 discard = (-pbf->bexponent) / 4; 212 break; 213 } 214 215 #ifdef DEBUG 216 { 217 long basexp; 218 int id; 219 220 printf(" step %d x operand length %d \n", i, lengthx); 221 _display_big_float(pbf, base); 222 printf(" step %d p operand length %d power %d \n", i, lengthp, tablepower[i]); 223 basexp = (base == 2) ? (lz[i])[tablepower[i]] : 0; 224 for (id = 0; id < lengthp; id++) { 225 printf("+ %d * ", (table[i])[id + (start[i])[tablepower[i]]]); 226 if (base == 2) 227 printf("2**%d", 16 * (basexp + id)); 228 if (base == 10) 229 printf("10**%d", 4 * (basexp + id)); 230 if ((id % 4) == 3) 231 printf("\n"); 232 } 233 printf("\n"); 234 } 235 if ((i == itlast) && (discard >= 0)) 236 printf(" alternative discard %d digits \n", discard); 237 #endif 238 239 if (base == 2) { 240 sumlz += (lz[i])[tablepower[i]]; 241 pbf->bexponent += 16 * (lz[i])[tablepower[i]]; 242 } 243 if (lengthp == 1) { /* Special case - multiply by 244 * <= 10**4 or 2**13 */ 245 switch (base) { 246 case 10: 247 _multiply_base_ten_by_two(pbf, tablepower[i]); 248 break; 249 case 2: 250 _multiply_base_two(pbf, (_BIG_FLOAT_DIGIT) ((table[i])[tablepower[i]]), (unsigned long) 0); 251 break; 252 } 253 #ifdef DEBUG 254 assert(pbf->blength <= pbf->bsize); 255 #endif 256 } else if (lengthx == 1) { /* Special case of short 257 * multiplicand. */ 258 _BIG_FLOAT_DIGIT multiplier = pbf->bsignificand[0]; 259 260 _copy_big_float_digits(pbf->bsignificand, (unsigned short *) &((table[i])[(start[i])[tablepower[i]]]), lengthp); 261 pbf->blength = lengthp; 262 switch (base) { 263 case 10: 264 _multiply_base_ten(pbf, multiplier); 265 break; 266 case 2: 267 _multiply_base_two(pbf, multiplier, (unsigned long) 0); 268 break; 269 } 270 #ifdef DEBUG 271 assert(pbf->blength <= pbf->bsize); 272 #endif 273 } else {/* General case. */ 274 short unsigned canquit; 275 short unsigned excess; 276 277 /* 278 * The result will be accumulated in *pbf 279 * from most significant to least 280 * significant. 281 */ 282 283 /* Generate criterion for early termination. */ 284 switch (base) { 285 case 2: 286 canquit = (short unsigned)65536; 287 break; 288 case 10: 289 canquit = 10000; 290 break; 291 } 292 canquit -= 3 + ((lengthx < lengthp) ? lengthx : lengthp); 293 294 pbf->bsignificand[lengthx + lengthp - 1] = 0; /* Only gets filled by 295 * carries. */ 296 for (j = lengthx + lengthp - 2; j >= 0; j--) { 297 int istart = j - lengthp + 1, istop = lengthx - 1; 298 short unsigned lengthprod; 299 _BIG_FLOAT_DIGIT product[3]; 300 301 pp = (unsigned short *) &((table[i])[(start[i])[tablepower[i]]]); 302 if (j < istop) 303 istop = j; 304 if (0 > istart) 305 istart = 0; 306 307 switch (base) { 308 case 2: 309 _multiply_base_two_vector((short unsigned) (istop - istart + 1), &(pbf->bsignificand[istart]), &(pp[j - istop]), product); 310 if (product[2] != 0) 311 _carry_propagate_two((unsigned long) product[2], &(pbf->bsignificand[j + 2])); 312 if (product[1] != 0) 313 _carry_propagate_two((unsigned long) product[1], &(pbf->bsignificand[j + 1])); 314 break; 315 case 10: 316 _multiply_base_ten_vector((short unsigned) (istop - istart + 1), &(pbf->bsignificand[istart]), &(pp[j - istop]), product); 317 if (product[2] != 0) 318 _carry_propagate_ten((unsigned long) product[2], &(pbf->bsignificand[j + 2])); 319 if (product[1] != 0) 320 _carry_propagate_ten((unsigned long) product[1], &(pbf->bsignificand[j + 1])); 321 break; 322 } 323 pbf->bsignificand[j] = product[0]; 324 lengthprod = lengthx + lengthp; 325 if (pbf->bsignificand[lengthprod - 1] == 0) 326 lengthprod--; 327 if (lengthprod > needed_precision) 328 excess = lengthprod - needed_precision; 329 else 330 excess = 0; 331 if ((i == itlast) && ((j + 2) <= excess) && (pbf->bsignificand[j + 1] <= canquit) 332 && ((pbf->bsignificand[j + 1] | pbf->bsignificand[j]) != 0)) { 333 /* 334 * On the last 335 * multiplication, it's not 336 * necessary to develop the 337 * entire product, if further 338 * digits can't possibly 339 * affect significant digits, 340 * unless there's a chance 341 * the product might be 342 * exact! 343 */ 344 /* 345 * Note that the product 346 * might be exact if the j 347 * and j+1 terms are zero; if 348 * they are non-zero, then it 349 * won't be after they're 350 * discarded. 351 */ 352 353 excess = j + 2; /* Can discard j+1, j, 354 * ... 0 */ 355 #ifdef DEBUG 356 printf(" decided to quit early at j %d since s[j+1] is %d <= %d \n", j, pbf->bsignificand[j + 1], canquit); 357 printf(" s[j+2..j] are %d %d %d \n", pbf->bsignificand[j + 2], pbf->bsignificand[j + 1], pbf->bsignificand[j]); 358 printf(" requested precision %d needed_precision %d big digits out of %d \n", precision, needed_precision, lengthprod); 359 #endif 360 if ((discard >= 0) && ((j + 2) > (discard - (int) sumlz))) { 361 #ifdef DEBUG 362 printf(" early quit rejected because j+2 = %d > %d = discard \n", j + 2, discard); 363 #endif 364 goto pastdiscard; 365 } 366 pbf->bsignificand[excess] |= 1; /* Sticky bit on. */ 367 #ifdef DEBUG 368 printf(" discard %d digits - last gets %d \n", excess, pbf->bsignificand[excess]); 369 #endif 370 trailing_zeros_to_delete = excess; 371 goto donegeneral; 372 } 373 pastdiscard: ; 374 #ifdef DEBUG 375 /* 376 * else { printf(" early termination 377 * rejected at j %d since s[j+1] = 378 * %d, canquit = %d \n", j, 379 * pbf->bsignificand[j + 1], 380 * canquit); printf(" s[j+2..j] are 381 * %d %d %d \n", pbf->bsignificand[j 382 * + 2], pbf->bsignificand[j + 1], 383 * pbf->bsignificand[j]); printf(" 384 * requested precision %d 385 * needed_precision %d big digits out 386 * of %d \n", precision, 387 * needed_precision, lengthprod); } 388 */ 389 #endif 390 } 391 trailing_zeros_to_delete = 0; 392 donegeneral: 393 pbf->blength = lengthx + lengthp; 394 if (pbf->bsignificand[pbf->blength - 1] == 0) 395 pbf->blength--; 396 for (; pbf->bsignificand[trailing_zeros_to_delete] == 0; trailing_zeros_to_delete++); 397 /* 398 * Look for additional trailing zeros to 399 * delete. 400 */ 401 402 /* 403 * fix for bug 1070565; if too many trailing 404 * zeroes are deleted, we'll violate the 405 * assertion that bexponent is in [-3,+4] 406 */ 407 if (base == 10) { 408 int deletelimit=(1-((pbf->bexponent+3)/4)); 409 410 if ((int)trailing_zeros_to_delete > deletelimit) { 411 #ifdef DEBUG 412 printf("\n __x_power trailing zeros delete count lowered from %d to " 413 "%d \n", trailing_zeros_to_delete,deletelimit); 414 #endif 415 416 trailing_zeros_to_delete = deletelimit; 417 } 418 } 419 420 421 if (trailing_zeros_to_delete != 0) { 422 #ifdef DEBUG 423 printf(" %d trailing zeros deleted \n", trailing_zeros_to_delete); 424 #endif 425 _copy_big_float_digits(pbf->bsignificand, &(pbf->bsignificand[trailing_zeros_to_delete]), pbf->blength - trailing_zeros_to_delete); 426 pbf->blength -= trailing_zeros_to_delete; 427 switch (base) { 428 case 2: 429 pbf->bexponent += 16 * trailing_zeros_to_delete; 430 break; 431 case 10: 432 pbf->bexponent += 4 * trailing_zeros_to_delete; 433 break; 434 } 435 } 436 } 437 } 438 if ((pbfold != pbf) && (pbf->blength <= pbfold->bsize)) { /* Don't need that huge 439 * buffer after all! */ 440 #ifdef DEBUG 441 printf(" free called from times_power because final length %d <= %d original size \n", pbf->blength, pbfold->bsize); 442 #endif 443 444 /* Copy product to original buffer. */ 445 pbfold->blength = pbf->blength; 446 pbfold->bexponent = pbf->bexponent; 447 _copy_big_float_digits(pbfold->bsignificand, pbf->bsignificand, pbf->blength); 448 _free_big_float(*pnewbf); /* Free new buffer. */ 449 *pnewbf = pbfold; /* New buffer pointer now agrees with 450 * original. */ 451 } 452 ret: 453 return; 454 } 455