1 /* Common routines for _dtoa and strtod */ 2 3 #include "fconv.h" 4 5 #ifdef DEBUG 6 #include <stdio.h> 7 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 8 #endif 9 10 double 11 _tens[] = { 12 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 13 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 14 1e20, 1e21, 1e22 15 #ifdef VAX 16 , 1e23, 1e24 17 #endif 18 }; 19 20 21 #ifdef IEEE_Arith 22 double _bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 23 double _tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 }; 24 #else 25 #ifdef IBM 26 double _bigtens[] = { 1e16, 1e32, 1e64 }; 27 double _tinytens[] = { 1e-16, 1e-32, 1e-64 }; 28 #else 29 double _bigtens[] = { 1e16, 1e32 }; 30 double _tinytens[] = { 1e-16, 1e-32 }; 31 #endif 32 #endif 33 34 static Bigint *freelist[Kmax+1]; 35 36 Bigint * 37 _Balloc(int k) 38 { 39 int x; 40 Bigint *rv; 41 42 if (rv = freelist[k]) { 43 freelist[k] = rv->next; 44 } 45 else { 46 x = 1 << k; 47 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long)); 48 rv->k = k; 49 rv->maxwds = x; 50 } 51 rv->sign = rv->wds = 0; 52 return rv; 53 } 54 55 void 56 _Bfree(Bigint *v) 57 { 58 if (v) { 59 v->next = freelist[v->k]; 60 freelist[v->k] = v; 61 } 62 } 63 64 65 Bigint * 66 _multadd(Bigint *b, int m, int a) /* multiply by m and add a */ 67 { 68 int i, wds; 69 unsigned long *x, y; 70 #ifdef Pack_32 71 unsigned long xi, z; 72 #endif 73 Bigint *b1; 74 75 wds = b->wds; 76 x = b->x; 77 i = 0; 78 do { 79 #ifdef Pack_32 80 xi = *x; 81 y = (xi & 0xffff) * m + a; 82 z = (xi >> 16) * m + (y >> 16); 83 a = (int)(z >> 16); 84 *x++ = (z << 16) + (y & 0xffff); 85 #else 86 y = *x * m + a; 87 a = (int)(y >> 16); 88 *x++ = y & 0xffff; 89 #endif 90 } 91 while(++i < wds); 92 if (a) { 93 if (wds >= b->maxwds) { 94 b1 = Balloc(b->k+1); 95 Bcopy(b1, b); 96 Bfree(b); 97 b = b1; 98 } 99 b->x[wds++] = a; 100 b->wds = wds; 101 } 102 return b; 103 } 104 105 int 106 _hi0bits(register unsigned long x) 107 { 108 register int k = 0; 109 110 if (!(x & 0xffff0000)) { 111 k = 16; 112 x <<= 16; 113 } 114 if (!(x & 0xff000000)) { 115 k += 8; 116 x <<= 8; 117 } 118 if (!(x & 0xf0000000)) { 119 k += 4; 120 x <<= 4; 121 } 122 if (!(x & 0xc0000000)) { 123 k += 2; 124 x <<= 2; 125 } 126 if (!(x & 0x80000000)) { 127 k++; 128 if (!(x & 0x40000000)) 129 return 32; 130 } 131 return k; 132 } 133 134 static int 135 lo0bits(unsigned long *y) 136 { 137 register int k; 138 register unsigned long x = *y; 139 140 if (x & 7) { 141 if (x & 1) 142 return 0; 143 if (x & 2) { 144 *y = x >> 1; 145 return 1; 146 } 147 *y = x >> 2; 148 return 2; 149 } 150 k = 0; 151 if (!(x & 0xffff)) { 152 k = 16; 153 x >>= 16; 154 } 155 if (!(x & 0xff)) { 156 k += 8; 157 x >>= 8; 158 } 159 if (!(x & 0xf)) { 160 k += 4; 161 x >>= 4; 162 } 163 if (!(x & 0x3)) { 164 k += 2; 165 x >>= 2; 166 } 167 if (!(x & 1)) { 168 k++; 169 x >>= 1; 170 if (!x & 1) 171 return 32; 172 } 173 *y = x; 174 return k; 175 } 176 177 Bigint * 178 _i2b(int i) 179 { 180 Bigint *b; 181 182 b = Balloc(1); 183 b->x[0] = i; 184 b->wds = 1; 185 return b; 186 } 187 188 Bigint * 189 _mult(Bigint *a, Bigint *b) 190 { 191 Bigint *c; 192 int k, wa, wb, wc; 193 unsigned long carry, y, z; 194 unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 195 #ifdef Pack_32 196 unsigned long z2; 197 #endif 198 199 if (a->wds < b->wds) { 200 c = a; 201 a = b; 202 b = c; 203 } 204 k = a->k; 205 wa = a->wds; 206 wb = b->wds; 207 wc = wa + wb; 208 if (wc > a->maxwds) 209 k++; 210 c = Balloc(k); 211 for(x = c->x, xa = x + wc; x < xa; x++) 212 *x = 0; 213 xa = a->x; 214 xae = xa + wa; 215 xb = b->x; 216 xbe = xb + wb; 217 xc0 = c->x; 218 #ifdef Pack_32 219 for(; xb < xbe; xb++, xc0++) { 220 if (y = *xb & 0xffff) { 221 x = xa; 222 xc = xc0; 223 carry = 0; 224 do { 225 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 226 carry = z >> 16; 227 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 228 carry = z2 >> 16; 229 Storeinc(xc, z2, z); 230 } 231 while(x < xae); 232 *xc = carry; 233 } 234 if (y = *xb >> 16) { 235 x = xa; 236 xc = xc0; 237 carry = 0; 238 z2 = *xc; 239 do { 240 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 241 carry = z >> 16; 242 Storeinc(xc, z, z2); 243 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 244 carry = z2 >> 16; 245 } 246 while(x < xae); 247 *xc = z2; 248 } 249 } 250 #else 251 for(; xb < xbe; xc0++) { 252 if (y = *xb++) { 253 x = xa; 254 xc = xc0; 255 carry = 0; 256 do { 257 z = *x++ * y + *xc + carry; 258 carry = z >> 16; 259 *xc++ = z & 0xffff; 260 } 261 while(x < xae); 262 *xc = carry; 263 } 264 } 265 #endif 266 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 267 c->wds = wc; 268 return c; 269 } 270 271 static Bigint *p5s; 272 273 Bigint * 274 _pow5mult(Bigint *b, int k) 275 { 276 Bigint *b1, *p5, *p51; 277 int i; 278 static int p05[3] = { 5, 25, 125 }; 279 280 if (i = k & 3) 281 b = multadd(b, p05[i-1], 0); 282 283 if (!(k >>= 2)) 284 return b; 285 if (!(p5 = p5s)) { 286 /* first time */ 287 p5 = p5s = i2b(625); 288 p5->next = 0; 289 } 290 for(;;) { 291 if (k & 1) { 292 b1 = mult(b, p5); 293 Bfree(b); 294 b = b1; 295 } 296 if (!(k >>= 1)) 297 break; 298 if (!(p51 = p5->next)) { 299 p51 = p5->next = mult(p5,p5); 300 p51->next = 0; 301 } 302 p5 = p51; 303 } 304 return b; 305 } 306 307 Bigint * 308 _lshift(Bigint *b, int k) 309 { 310 int i, k1, n, n1; 311 Bigint *b1; 312 unsigned long *x, *x1, *xe, z; 313 314 #ifdef Pack_32 315 n = k >> 5; 316 #else 317 n = k >> 4; 318 #endif 319 k1 = b->k; 320 n1 = n + b->wds + 1; 321 for(i = b->maxwds; n1 > i; i <<= 1) 322 k1++; 323 b1 = Balloc(k1); 324 x1 = b1->x; 325 for(i = 0; i < n; i++) 326 *x1++ = 0; 327 x = b->x; 328 xe = x + b->wds; 329 #ifdef Pack_32 330 if (k &= 0x1f) { 331 k1 = 32 - k; 332 z = 0; 333 do { 334 *x1++ = *x << k | z; 335 z = *x++ >> k1; 336 } 337 while(x < xe); 338 if (*x1 = z) 339 ++n1; 340 } 341 #else 342 if (k &= 0xf) { 343 k1 = 16 - k; 344 z = 0; 345 do { 346 *x1++ = *x << k & 0xffff | z; 347 z = *x++ >> k1; 348 } 349 while(x < xe); 350 if (*x1 = z) 351 ++n1; 352 } 353 #endif 354 else do 355 *x1++ = *x++; 356 while(x < xe); 357 b1->wds = n1 - 1; 358 Bfree(b); 359 return b1; 360 } 361 362 int 363 _cmp(Bigint *a, Bigint *b) 364 { 365 unsigned long *xa, *xa0, *xb, *xb0; 366 int i, j; 367 368 i = a->wds; 369 j = b->wds; 370 #ifdef DEBUG 371 if (i > 1 && !a->x[i-1]) 372 Bug("cmp called with a->x[a->wds-1] == 0"); 373 if (j > 1 && !b->x[j-1]) 374 Bug("cmp called with b->x[b->wds-1] == 0"); 375 #endif 376 if (i -= j) 377 return i; 378 xa0 = a->x; 379 xa = xa0 + j; 380 xb0 = b->x; 381 xb = xb0 + j; 382 for(;;) { 383 if (*--xa != *--xb) 384 return *xa < *xb ? -1 : 1; 385 if (xa <= xa0) 386 break; 387 } 388 return 0; 389 } 390 391 Bigint * 392 _diff(Bigint *a, Bigint *b) 393 { 394 Bigint *c; 395 int i, wa, wb; 396 long borrow, y; /* We need signed shifts here. */ 397 unsigned long *xa, *xae, *xb, *xbe, *xc; 398 #ifdef Pack_32 399 long z; 400 #endif 401 402 i = cmp(a,b); 403 if (!i) { 404 c = Balloc(0); 405 c->wds = 1; 406 c->x[0] = 0; 407 return c; 408 } 409 if (i < 0) { 410 c = a; 411 a = b; 412 b = c; 413 i = 1; 414 } 415 else 416 i = 0; 417 c = Balloc(a->k); 418 c->sign = i; 419 wa = a->wds; 420 xa = a->x; 421 xae = xa + wa; 422 wb = b->wds; 423 xb = b->x; 424 xbe = xb + wb; 425 xc = c->x; 426 borrow = 0; 427 #ifdef Pack_32 428 do { 429 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; 430 borrow = y >> 16; 431 Sign_Extend(borrow, y); 432 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; 433 borrow = z >> 16; 434 Sign_Extend(borrow, z); 435 Storeinc(xc, z, y); 436 } 437 while(xb < xbe); 438 while(xa < xae) { 439 y = (*xa & 0xffff) + borrow; 440 borrow = y >> 16; 441 Sign_Extend(borrow, y); 442 z = (*xa++ >> 16) + borrow; 443 borrow = z >> 16; 444 Sign_Extend(borrow, z); 445 Storeinc(xc, z, y); 446 } 447 #else 448 do { 449 y = *xa++ - *xb++ + borrow; 450 borrow = y >> 16; 451 Sign_Extend(borrow, y); 452 *xc++ = y & 0xffff; 453 } 454 while(xb < xbe); 455 while(xa < xae) { 456 y = *xa++ + borrow; 457 borrow = y >> 16; 458 Sign_Extend(borrow, y); 459 *xc++ = y & 0xffff; 460 } 461 #endif 462 while(!*--xc) 463 wa--; 464 c->wds = wa; 465 return c; 466 } 467 468 Bigint * 469 _d2b(double darg, int *e, int *bits) 470 { 471 Bigint *b; 472 int de, i, k; 473 unsigned long *x, y, z; 474 Dul d; 475 #ifdef VAX 476 unsigned long d0, d1; 477 d.d = darg; 478 d0 = word0(d) >> 16 | word0(d) << 16; 479 d1 = word1(d) >> 16 | word1(d) << 16; 480 #else 481 d.d = darg; 482 #define d0 word0(d) 483 #define d1 word1(d) 484 #endif 485 486 #ifdef Pack_32 487 b = Balloc(1); 488 #else 489 b = Balloc(2); 490 #endif 491 x = b->x; 492 493 z = d0 & Frac_mask; 494 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 495 #ifdef Sudden_Underflow 496 de = (int)(d0 >> Exp_shift); 497 #ifndef IBM 498 z |= Exp_msk11; 499 #endif 500 #else 501 if (de = (int)(d0 >> Exp_shift)) 502 z |= Exp_msk1; 503 #endif 504 #ifdef Pack_32 505 if (y = d1) { 506 if (k = lo0bits(&y)) { 507 x[0] = y | z << 32 - k; 508 z >>= k; 509 } 510 else 511 x[0] = y; 512 i = b->wds = (x[1] = z) ? 2 : 1; 513 } 514 else { 515 #ifdef DEBUG 516 if (!z) 517 Bug("Zero passed to d2b"); 518 #endif 519 k = lo0bits(&z); 520 x[0] = z; 521 i = b->wds = 1; 522 k += 32; 523 } 524 #else 525 if (y = d1) { 526 if (k = lo0bits(&y)) 527 if (k >= 16) { 528 x[0] = y | z << 32 - k & 0xffff; 529 x[1] = z >> k - 16 & 0xffff; 530 x[2] = z >> k; 531 i = 2; 532 } 533 else { 534 x[0] = y & 0xffff; 535 x[1] = y >> 16 | z << 16 - k & 0xffff; 536 x[2] = z >> k & 0xffff; 537 x[3] = z >> k+16; 538 i = 3; 539 } 540 else { 541 x[0] = y & 0xffff; 542 x[1] = y >> 16; 543 x[2] = z & 0xffff; 544 x[3] = z >> 16; 545 i = 3; 546 } 547 } 548 else { 549 #ifdef DEBUG 550 if (!z) 551 Bug("Zero passed to d2b"); 552 #endif 553 k = lo0bits(&z); 554 if (k >= 16) { 555 x[0] = z; 556 i = 0; 557 } 558 else { 559 x[0] = z & 0xffff; 560 x[1] = z >> 16; 561 i = 1; 562 } 563 k += 32; 564 } 565 while(!x[i]) 566 --i; 567 b->wds = i + 1; 568 #endif 569 #ifndef Sudden_Underflow 570 if (de) { 571 #endif 572 #ifdef IBM 573 *e = (de - Bias - (P-1) << 2) + k; 574 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 575 #else 576 *e = de - Bias - (P-1) + k; 577 *bits = P - k; 578 #endif 579 #ifndef Sudden_Underflow 580 } 581 else { 582 *e = de - Bias - (P-1) + 1 + k; 583 #ifdef Pack_32 584 *bits = 32*i - hi0bits(x[i-1]); 585 #else 586 *bits = (i+2)*16 - hi0bits(x[i]); 587 #endif 588 } 589 #endif 590 return b; 591 } 592 #undef d0 593 #undef d1 594