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 (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 22 /* 23 * Copyright 2008 Sun Microsystems, Inc. All rights reserved. 24 * Use is subject to license terms. 25 */ 26 27 #include "lint.h" 28 #include "base_conversion.h" 29 #include <sys/types.h> 30 #include <malloc.h> 31 #include <memory.h> 32 #include <stdlib.h> 33 #include <errno.h> 34 35 /* 36 * Multiply a _big_float by a power of two or ten 37 */ 38 39 /* see comment in double_decim.c */ 40 static unsigned int 41 __quorem10000(unsigned int x, unsigned short *pr) 42 { 43 *pr = x % 10000; 44 return (x / 10000); 45 } 46 47 /* 48 * Multiply a base-2**16 significand by multiplier. Extend length as 49 * necessary to accommodate carries. 50 */ 51 static void 52 __multiply_base_two(_big_float *pbf, unsigned short multiplier) 53 { 54 unsigned int p, carry; 55 int j, length = pbf->blength; 56 57 carry = 0; 58 for (j = 0; j < length; j++) { 59 p = (unsigned int)pbf->bsignificand[j] * multiplier + carry; 60 pbf->bsignificand[j] = p & 0xffff; 61 carry = p >> 16; 62 } 63 if (carry != 0) 64 pbf->bsignificand[j++] = carry; 65 pbf->blength = j; 66 } 67 68 /* 69 * Multiply a base-10**4 significand by multiplier. Extend length as 70 * necessary to accommodate carries. 71 */ 72 static void 73 __multiply_base_ten(_big_float *pbf, unsigned short multiplier) 74 { 75 unsigned int p, carry; 76 int j, length = pbf->blength; 77 78 carry = 0; 79 for (j = 0; j < length; j++) { 80 p = (unsigned int)pbf->bsignificand[j] * multiplier + carry; 81 carry = __quorem10000(p, &pbf->bsignificand[j]); 82 } 83 while (carry != 0) { 84 carry = __quorem10000(carry, &pbf->bsignificand[j]); 85 j++; 86 } 87 pbf->blength = j; 88 } 89 90 /* 91 * Multiply a base-10**4 significand by 2**multiplier. Extend length 92 * as necessary to accommodate carries. 93 */ 94 static void 95 __multiply_base_ten_by_two(_big_float *pbf, unsigned short multiplier) 96 { 97 unsigned int p, carry; 98 int j, length = pbf->blength; 99 100 carry = 0; 101 for (j = 0; j < length; j++) { 102 p = ((unsigned int)pbf->bsignificand[j] << multiplier) + carry; 103 carry = __quorem10000(p, &pbf->bsignificand[j]); 104 } 105 while (carry != 0) { 106 carry = __quorem10000(carry, &pbf->bsignificand[j]); 107 j++; 108 } 109 pbf->blength = j; 110 } 111 112 /* 113 * Propagate carries in a base-2**16 significand. 114 */ 115 static void 116 __carry_propagate_two(unsigned int carry, unsigned short *psignificand) 117 { 118 unsigned int p; 119 int j; 120 121 j = 0; 122 while (carry != 0) { 123 p = psignificand[j] + carry; 124 psignificand[j++] = p & 0xffff; 125 carry = p >> 16; 126 } 127 } 128 129 /* 130 * Propagate carries in a base-10**4 significand. 131 */ 132 static void 133 __carry_propagate_ten(unsigned int carry, unsigned short *psignificand) 134 { 135 unsigned int p; 136 int j; 137 138 j = 0; 139 while (carry != 0) { 140 p = psignificand[j] + carry; 141 carry = __quorem10000(p, &psignificand[j]); 142 j++; 143 } 144 } 145 146 /* 147 * Given x[] and y[], base-2**16 vectors of length n, compute the 148 * dot product 149 * 150 * sum (i=0,n-1) of x[i]*y[n-1-i] 151 * 152 * The product may fill as many as three base-2**16 digits; product[0] 153 * is the least significant, product[2] the most. 154 */ 155 static void 156 __multiply_base_two_vector(unsigned short n, unsigned short *px, 157 unsigned short *py, unsigned short *product) 158 { 159 160 unsigned int acc, p; 161 unsigned short carry; 162 int i; 163 164 acc = 0; 165 carry = 0; 166 for (i = 0; i < (int)n; i++) { 167 p = px[i] * py[n - 1 - i] + acc; 168 if (p < acc) 169 carry++; 170 acc = p; 171 } 172 product[0] = acc & 0xffff; 173 product[1] = acc >> 16; 174 product[2] = carry; 175 } 176 177 /* 178 * Given x[] and y[], base-10**4 vectors of length n, compute the 179 * dot product 180 * 181 * sum (i=0,n-1) of x[i]*y[n-1-i] 182 * 183 * The product may fill as many as three base-10**4 digits; product[0] 184 * is the least significant, product[2] the most. 185 */ 186 #define ABASE 3000000000u /* base of accumulator */ 187 188 static void 189 __multiply_base_ten_vector(unsigned short n, unsigned short *px, 190 unsigned short *py, unsigned short *product) 191 { 192 193 unsigned int acc; 194 unsigned short carry; 195 int i; 196 197 acc = 0; 198 carry = 0; 199 for (i = 0; i < (int)n; i++) { 200 acc = px[i] * py[n - 1 - i] + acc; 201 if (acc >= ABASE) { 202 carry++; 203 acc -= ABASE; 204 } 205 } 206 product[0] = acc % 10000; 207 acc = acc / 10000; 208 product[1] = acc % 10000; 209 acc = acc / 10000; 210 product[2] = acc + (ABASE / 100000000) * carry; 211 } 212 213 /* 214 * Multiply *pbf by the n-th power of mult, which must be two or 215 * ten. If mult is two, *pbf is assumed to be base 10**4; if mult 216 * is ten, *pbf is assumed to be base 2**16. precision specifies 217 * the number of significant bits or decimal digits required in the 218 * result. (The product may have more or fewer digits than this, 219 * but it will be accurate to at least this many.) 220 * 221 * On exit, if the product is small enough, it overwrites *pbf, and 222 * *pnewbf is set to pbf. If the product is too large to fit in *pbf, 223 * this routine calls malloc(3M) to allocate storage and sets *pnewbf 224 * to point to this area; it is the caller's responsibility to free 225 * this storage when it is no longer needed. Note that *pbf may be 226 * modified even when the routine allocates new storage. 227 * 228 * If n is too large, we set errno to ERANGE and call abort(3C). 229 * If an attempted malloc fails, we set errno to ENOMEM and call 230 * abort(3C). 231 */ 232 void 233 __big_float_times_power(_big_float *pbf, int mult, int n, int precision, 234 _big_float **pnewbf) 235 { 236 int base, needed_precision, productsize; 237 int i, j, itlast, trailing_zeros_to_delete; 238 int tablepower[3], length[3]; 239 int lengthx, lengthp, istart, istop; 240 int excess_check; 241 unsigned short *pp, *table[3], canquit; 242 unsigned short multiplier, product[3]; 243 244 if (pbf->blength == 0) { 245 *pnewbf = pbf; 246 return; 247 } 248 249 if (mult == 2) { 250 /* 251 * Handle small n cases that don't require extra 252 * storage quickly. 253 */ 254 if (n <= 16 && pbf->blength + 2 < pbf->bsize) { 255 __multiply_base_ten_by_two(pbf, n); 256 *pnewbf = pbf; 257 return; 258 } 259 260 /* *pbf is base 10**4 */ 261 base = 10; 262 263 /* 264 * Convert precision (base ten digits) to needed_precision 265 * (base 10**4 digits), allowing an additional digit at 266 * each end. 267 */ 268 needed_precision = 2 + (precision >> 2); 269 270 /* 271 * Set up pointers to the table entries and compute their 272 * lengths. 273 */ 274 if (n < __TBL_2_SMALL_SIZE) { 275 itlast = 0; 276 tablepower[0] = n; 277 tablepower[1] = tablepower[2] = 0; 278 } else if (n < (__TBL_2_SMALL_SIZE * __TBL_2_BIG_SIZE)) { 279 itlast = 1; 280 tablepower[0] = n % __TBL_2_SMALL_SIZE; 281 tablepower[1] = n / __TBL_2_SMALL_SIZE; 282 tablepower[2] = 0; 283 } else if (n < (__TBL_2_SMALL_SIZE * __TBL_2_BIG_SIZE * 284 __TBL_2_HUGE_SIZE)) { 285 itlast = 2; 286 tablepower[0] = n % __TBL_2_SMALL_SIZE; 287 n /= __TBL_2_SMALL_SIZE; 288 tablepower[1] = n % __TBL_2_BIG_SIZE; 289 tablepower[2] = n / __TBL_2_BIG_SIZE; 290 } else { 291 errno = ERANGE; 292 abort(); 293 } 294 pp = (unsigned short *)__tbl_2_small_start + tablepower[0]; 295 table[0] = (unsigned short *)__tbl_2_small_digits + pp[0]; 296 length[0] = pp[1] - pp[0]; 297 pp = (unsigned short *)__tbl_2_big_start + tablepower[1]; 298 table[1] = (unsigned short *)__tbl_2_big_digits + pp[0]; 299 length[1] = pp[1] - pp[0]; 300 pp = (unsigned short *)__tbl_2_huge_start + tablepower[2]; 301 table[2] = (unsigned short *)__tbl_2_huge_digits + pp[0]; 302 length[2] = pp[1] - pp[0]; 303 } else { 304 if (n <= 4 && pbf->blength + 1 < pbf->bsize) { 305 pbf->bexponent += (short)n; 306 __multiply_base_two(pbf, __tbl_10_small_digits[n]); 307 *pnewbf = pbf; 308 return; 309 } 310 311 /* *pbf is base 2**16 */ 312 base = 2; 313 pbf->bexponent += (short)n; /* now need to multiply by 5**n */ 314 needed_precision = 2 + (precision >> 4); 315 if (n < __TBL_10_SMALL_SIZE) { 316 itlast = 0; 317 tablepower[0] = n; 318 tablepower[1] = tablepower[2] = 0; 319 } else if (n < (__TBL_10_SMALL_SIZE * __TBL_10_BIG_SIZE)) { 320 itlast = 1; 321 tablepower[0] = n % __TBL_10_SMALL_SIZE; 322 tablepower[1] = n / __TBL_10_SMALL_SIZE; 323 tablepower[2] = 0; 324 } else if (n < (__TBL_10_SMALL_SIZE * __TBL_10_BIG_SIZE * 325 __TBL_10_HUGE_SIZE)) { 326 itlast = 2; 327 tablepower[0] = n % __TBL_10_SMALL_SIZE; 328 n /= __TBL_10_SMALL_SIZE; 329 tablepower[1] = n % __TBL_10_BIG_SIZE; 330 tablepower[2] = n / __TBL_10_BIG_SIZE; 331 } else { 332 errno = ERANGE; 333 abort(); 334 } 335 pp = (unsigned short *)__tbl_10_small_start + tablepower[0]; 336 table[0] = (unsigned short *)__tbl_10_small_digits + pp[0]; 337 length[0] = pp[1] - pp[0]; 338 pp = (unsigned short *)__tbl_10_big_start + tablepower[1]; 339 table[1] = (unsigned short *)__tbl_10_big_digits + pp[0]; 340 length[1] = pp[1] - pp[0]; 341 pp = (unsigned short *)__tbl_10_huge_start + tablepower[2]; 342 table[2] = (unsigned short *)__tbl_10_huge_digits + pp[0]; 343 length[2] = pp[1] - pp[0]; 344 } 345 346 /* compute an upper bound on the size of the product */ 347 productsize = pbf->blength; 348 for (i = 0; i <= itlast; i++) 349 productsize += length[i]; 350 351 if (productsize < needed_precision) 352 needed_precision = productsize; 353 354 if (productsize <= pbf->bsize) { 355 *pnewbf = pbf; 356 } else { 357 i = sizeof (_big_float) + sizeof (unsigned short) * 358 (productsize - _BIG_FLOAT_SIZE); 359 if ((*pnewbf = malloc(i)) == NULL) { 360 errno = ENOMEM; 361 abort(); 362 } 363 (void) memcpy((*pnewbf)->bsignificand, pbf->bsignificand, 364 pbf->blength * sizeof (unsigned short)); 365 (*pnewbf)->blength = pbf->blength; 366 (*pnewbf)->bexponent = pbf->bexponent; 367 pbf = *pnewbf; 368 pbf->bsize = productsize; 369 } 370 371 /* 372 * Now pbf points to the input and the output. Step through 373 * each level of the tables. 374 */ 375 for (i = 0; i <= itlast; i++) { 376 if (tablepower[i] == 0) 377 continue; 378 379 lengthp = length[i]; 380 if (lengthp == 1) { 381 /* short multiplier (<= 10**4 or 2**13) */ 382 if (base == 10) { 383 /* left shift by tablepower[i] */ 384 __multiply_base_ten_by_two(pbf, tablepower[i]); 385 } else { 386 __multiply_base_two(pbf, (table[i])[0]); 387 } 388 continue; 389 } 390 391 lengthx = pbf->blength; 392 if (lengthx == 1) { 393 /* short multiplicand */ 394 multiplier = pbf->bsignificand[0]; 395 (void) memcpy(pbf->bsignificand, table[i], 396 lengthp * sizeof (unsigned short)); 397 pbf->blength = lengthp; 398 if (base == 10) 399 __multiply_base_ten(pbf, multiplier); 400 else 401 __multiply_base_two(pbf, multiplier); 402 continue; 403 } 404 405 /* keep track of trailing zeroes */ 406 trailing_zeros_to_delete = 0; 407 408 /* initialize for carry propagation */ 409 pbf->bsignificand[lengthx + lengthp - 1] = 0; 410 411 /* 412 * General case - the result will be accumulated in *pbf 413 * from most significant digit to least significant. 414 */ 415 for (j = lengthx + lengthp - 2; j >= 0; j--) { 416 istart = j - lengthp + 1; 417 if (istart < 0) 418 istart = 0; 419 420 istop = lengthx - 1; 421 if (istop > j) 422 istop = j; 423 424 pp = table[i]; 425 if (base == 2) { 426 __multiply_base_two_vector(istop - istart + 1, 427 &(pbf->bsignificand[istart]), 428 &(pp[j - istop]), product); 429 if (product[2] != 0) 430 __carry_propagate_two( 431 (unsigned int)product[2], 432 &(pbf->bsignificand[j + 2])); 433 if (product[1] != 0) 434 __carry_propagate_two( 435 (unsigned int)product[1], 436 &(pbf->bsignificand[j + 1])); 437 } else { 438 __multiply_base_ten_vector(istop - istart + 1, 439 &(pbf->bsignificand[istart]), 440 &(pp[j - istop]), product); 441 if (product[2] != 0) 442 __carry_propagate_ten( 443 (unsigned int)product[2], 444 &(pbf->bsignificand[j + 2])); 445 if (product[1] != 0) 446 __carry_propagate_ten( 447 (unsigned int)product[1], 448 &(pbf->bsignificand[j + 1])); 449 } 450 pbf->bsignificand[j] = product[0]; 451 if (i < itlast || j > lengthx + lengthp - 4 452 - needed_precision) 453 continue; 454 455 /* 456 * On the last multiplication, it's not necessary 457 * to develop the entire product if further digits 458 * can't possibly affect significant digits. But 459 * note that further digits can affect the product 460 * in one of two ways: (i) the sum of digits beyond 461 * the significant ones can cause a carry that would 462 * propagate into the significant digits, or (ii) no 463 * carry will occur, but there may be more nonzero 464 * digits that will need to be recorded in a sticky 465 * bit. 466 */ 467 excess_check = lengthx + lengthp; 468 if (pbf->bsignificand[excess_check - 1] == 0) 469 excess_check--; 470 excess_check -= needed_precision + 4; 471 canquit = ((base == 2)? 65535 : 9999) - 472 ((lengthx < lengthp)? lengthx : lengthp); 473 /* 474 * If j <= excess_check, then we have all the 475 * significant digits. If the (j + 1)-st digit 476 * is no larger than canquit, then the sum of the 477 * digits not yet computed can't carry into the 478 * significant digits. If the j-th and (j + 1)-st 479 * digits are not both zero, then we know we are 480 * discarding nonzero digits. (If both of these 481 * digits are zero, we need to keep forming more 482 * of the product to see whether or not there are 483 * any more nonzero digits.) 484 */ 485 if (j <= excess_check && 486 pbf->bsignificand[j + 1] <= canquit && 487 (pbf->bsignificand[j + 1] | pbf->bsignificand[j]) 488 != 0) { 489 /* can discard j+1, j, ... 0 */ 490 trailing_zeros_to_delete = j + 2; 491 492 /* set sticky bit */ 493 pbf->bsignificand[j + 2] |= 1; 494 break; 495 } 496 } 497 498 /* if the product didn't carry, delete the leading zero */ 499 pbf->blength = lengthx + lengthp; 500 if (pbf->bsignificand[pbf->blength - 1] == 0) 501 pbf->blength--; 502 503 /* look for additional trailing zeros to delete */ 504 for (; pbf->bsignificand[trailing_zeros_to_delete] == 0; 505 trailing_zeros_to_delete++) 506 continue; 507 508 if (trailing_zeros_to_delete > 0) { 509 for (j = 0; j < (int)pbf->blength - 510 trailing_zeros_to_delete; j++) { 511 pbf->bsignificand[j] = pbf->bsignificand[j 512 + trailing_zeros_to_delete]; 513 } 514 pbf->blength -= trailing_zeros_to_delete; 515 pbf->bexponent += trailing_zeros_to_delete << 516 ((base == 2)? 4 : 2); 517 } 518 } 519 } 520