contrib/gdtoa/misc.c

   1 /****************************************************************
   2
   3 The author of this software is David M. Gay.
   4
   5 Copyright (C) 1998, 1999 by Lucent Technologies
   6 All Rights Reserved
   7
   8 Permission to use, copy, modify, and distribute this software and
   9 its documentation for any purpose and without fee is hereby
  10 granted, provided that the above copyright notice appear in all
  11 copies and that both that the copyright notice and this
  12 permission notice and warranty disclaimer appear in supporting
  13 documentation, and that the name of Lucent or any of its entities
  14 not be used in advertising or publicity pertaining to
  15 distribution of the software without specific, written prior
  16 permission.
  17
  18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  25 THIS SOFTWARE.
  26
  27 ****************************************************************/
  28
  29 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  30  * with " at " changed at "@" and " dot " changed to ".").      */
  31
  32 #include "gdtoaimp.h"
  33
  34  static Bigint *freelist[Kmax+1];
  35 #ifndef Omit_Private_Memory
  36 #ifndef PRIVATE_MEM
  37 #define PRIVATE_MEM 2304
  38 #endif
  39 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
  40 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
  41 #endif
  42
  43  Bigint *
  44 Balloc
  45 #ifdef KR_headers
  46         (k) int k;
  47 #else
  48         (int k)
  49 #endif
  50 {
  51         int x;
  52         Bigint *rv;
  53 #ifndef Omit_Private_Memory
  54         unsigned int len;
  55 #endif
  56
  57         ACQUIRE_DTOA_LOCK(0);
  58         /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
  59         /* but this case seems very unlikely. */
  60         if (k <= Kmax && (rv = freelist[k]) !=0) {
  61                 freelist[k] = rv->next;
  62                 }
  63         else {
  64                 x = 1 << k;
  65 #ifdef Omit_Private_Memory
  66                 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
  67 #else
  68                 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
  69                         /sizeof(double);
  70                 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
  71                         rv = (Bigint*)pmem_next;
  72                         pmem_next += len;
  73                         }
  74                 else
  75                         rv = (Bigint*)MALLOC(len*sizeof(double));
  76 #endif
  77                 rv->k = k;
  78                 rv->maxwds = x;
  79                 }
  80         FREE_DTOA_LOCK(0);
  81         rv->sign = rv->wds = 0;
  82         return rv;
  83         }
  84
  85  void
  86 Bfree
  87 #ifdef KR_headers
  88         (v) Bigint *v;
  89 #else
  90         (Bigint *v)
  91 #endif
  92 {
  93         if (v) {
  94                 if (v->k > Kmax)
  95 #ifdef FREE
  96                         FREE((void*)v);
  97 #else
  98                         free((void*)v);
  99 #endif
 100                 else {
 101                         ACQUIRE_DTOA_LOCK(0);
 102                         v->next = freelist[v->k];
 103                         freelist[v->k] = v;
 104                         FREE_DTOA_LOCK(0);
 105                         }
 106                 }
 107         }
 108
 109  int
 110 lo0bits
 111 #ifdef KR_headers
 112         (y) ULong *y;
 113 #else
 114         (ULong *y)
 115 #endif
 116 {
 117         int k;
 118         ULong x = *y;
 119
 120         if (x & 7) {
 121                 if (x & 1)
 122                         return 0;
 123                 if (x & 2) {
 124                         *y = x >> 1;
 125                         return 1;
 126                         }
 127                 *y = x >> 2;
 128                 return 2;
 129                 }
 130         k = 0;
 131         if (!(x & 0xffff)) {
 132                 k = 16;
 133                 x >>= 16;
 134                 }
 135         if (!(x & 0xff)) {
 136                 k += 8;
 137                 x >>= 8;
 138                 }
 139         if (!(x & 0xf)) {
 140                 k += 4;
 141                 x >>= 4;
 142                 }
 143         if (!(x & 0x3)) {
 144                 k += 2;
 145                 x >>= 2;
 146                 }
 147         if (!(x & 1)) {
 148                 k++;
 149                 x >>= 1;
 150                 if (!x)
 151                         return 32;
 152                 }
 153         *y = x;
 154         return k;
 155         }
 156
 157  Bigint *
 158 multadd
 159 #ifdef KR_headers
 160         (b, m, a) Bigint *b; int m, a;
 161 #else
 162         (Bigint *b, int m, int a)       /* multiply by m and add a */
 163 #endif
 164 {
 165         int i, wds;
 166 #ifdef ULLong
 167         ULong *x;
 168         ULLong carry, y;
 169 #else
 170         ULong carry, *x, y;
 171 #ifdef Pack_32
 172         ULong xi, z;
 173 #endif
 174 #endif
 175         Bigint *b1;
 176
 177         wds = b->wds;
 178         x = b->x;
 179         i = 0;
 180         carry = a;
 181         do {
 182 #ifdef ULLong
 183                 y = *x * (ULLong)m + carry;
 184                 carry = y >> 32;
 185                 *x++ = y & 0xffffffffUL;
 186 #else
 187 #ifdef Pack_32
 188                 xi = *x;
 189                 y = (xi & 0xffff) * m + carry;
 190                 z = (xi >> 16) * m + (y >> 16);
 191                 carry = z >> 16;
 192                 *x++ = (z << 16) + (y & 0xffff);
 193 #else
 194                 y = *x * m + carry;
 195                 carry = y >> 16;
 196                 *x++ = y & 0xffff;
 197 #endif
 198 #endif
 199                 }
 200                 while(++i < wds);
 201         if (carry) {
 202                 if (wds >= b->maxwds) {
 203                         b1 = Balloc(b->k+1);
 204                         Bcopy(b1, b);
 205                         Bfree(b);
 206                         b = b1;
 207                         }
 208                 b->x[wds++] = carry;
 209                 b->wds = wds;
 210                 }
 211         return b;
 212         }
 213
 214  int
 215 hi0bits_D2A
 216 #ifdef KR_headers
 217         (x) ULong x;
 218 #else
 219         (ULong x)
 220 #endif
 221 {
 222         int k = 0;
 223
 224         if (!(x & 0xffff0000)) {
 225                 k = 16;
 226                 x <<= 16;
 227                 }
 228         if (!(x & 0xff000000)) {
 229                 k += 8;
 230                 x <<= 8;
 231                 }
 232         if (!(x & 0xf0000000)) {
 233                 k += 4;
 234                 x <<= 4;
 235                 }
 236         if (!(x & 0xc0000000)) {
 237                 k += 2;
 238                 x <<= 2;
 239                 }
 240         if (!(x & 0x80000000)) {
 241                 k++;
 242                 if (!(x & 0x40000000))
 243                         return 32;
 244                 }
 245         return k;
 246         }
 247
 248  Bigint *
 249 i2b
 250 #ifdef KR_headers
 251         (i) int i;
 252 #else
 253         (int i)
 254 #endif
 255 {
 256         Bigint *b;
 257
 258         b = Balloc(1);
 259         b->x[0] = i;
 260         b->wds = 1;
 261         return b;
 262         }
 263
 264  Bigint *
 265 mult
 266 #ifdef KR_headers
 267         (a, b) Bigint *a, *b;
 268 #else
 269         (Bigint *a, Bigint *b)
 270 #endif
 271 {
 272         Bigint *c;
 273         int k, wa, wb, wc;
 274         ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
 275         ULong y;
 276 #ifdef ULLong
 277         ULLong carry, z;
 278 #else
 279         ULong carry, z;
 280 #ifdef Pack_32
 281         ULong z2;
 282 #endif
 283 #endif
 284
 285         if (a->wds < b->wds) {
 286                 c = a;
 287                 a = b;
 288                 b = c;
 289                 }
 290         k = a->k;
 291         wa = a->wds;
 292         wb = b->wds;
 293         wc = wa + wb;
 294         if (wc > a->maxwds)
 295                 k++;
 296         c = Balloc(k);
 297         for(x = c->x, xa = x + wc; x < xa; x++)
 298                 *x = 0;
 299         xa = a->x;
 300         xae = xa + wa;
 301         xb = b->x;
 302         xbe = xb + wb;
 303         xc0 = c->x;
 304 #ifdef ULLong
 305         for(; xb < xbe; xc0++) {
 306                 if ( (y = *xb++) !=0) {
 307                         x = xa;
 308                         xc = xc0;
 309                         carry = 0;
 310                         do {
 311                                 z = *x++ * (ULLong)y + *xc + carry;
 312                                 carry = z >> 32;
 313                                 *xc++ = z & 0xffffffffUL;
 314                                 }
 315                                 while(x < xae);
 316                         *xc = carry;
 317                         }
 318                 }
 319 #else
 320 #ifdef Pack_32
 321         for(; xb < xbe; xb++, xc0++) {
 322                 if ( (y = *xb & 0xffff) !=0) {
 323                         x = xa;
 324                         xc = xc0;
 325                         carry = 0;
 326                         do {
 327                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
 328                                 carry = z >> 16;
 329                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
 330                                 carry = z2 >> 16;
 331                                 Storeinc(xc, z2, z);
 332                                 }
 333                                 while(x < xae);
 334                         *xc = carry;
 335                         }
 336                 if ( (y = *xb >> 16) !=0) {
 337                         x = xa;
 338                         xc = xc0;
 339                         carry = 0;
 340                         z2 = *xc;
 341                         do {
 342                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
 343                                 carry = z >> 16;
 344                                 Storeinc(xc, z, z2);
 345                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
 346                                 carry = z2 >> 16;
 347                                 }
 348                                 while(x < xae);
 349                         *xc = z2;
 350                         }
 351                 }
 352 #else
 353         for(; xb < xbe; xc0++) {
 354                 if ( (y = *xb++) !=0) {
 355                         x = xa;
 356                         xc = xc0;
 357                         carry = 0;
 358                         do {
 359                                 z = *x++ * y + *xc + carry;
 360                                 carry = z >> 16;
 361                                 *xc++ = z & 0xffff;
 362                                 }
 363                                 while(x < xae);
 364                         *xc = carry;
 365                         }
 366                 }
 367 #endif
 368 #endif
 369         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
 370         c->wds = wc;
 371         return c;
 372         }
 373
 374  static Bigint *p5s;
 375
 376  Bigint *
 377 pow5mult
 378 #ifdef KR_headers
 379         (b, k) Bigint *b; int k;
 380 #else
 381         (Bigint *b, int k)
 382 #endif
 383 {
 384         Bigint *b1, *p5, *p51;
 385         int i;
 386         static int p05[3] = { 5, 25, 125 };
 387
 388         if ( (i = k & 3) !=0)
 389                 b = multadd(b, p05[i-1], 0);
 390
 391         if (!(k >>= 2))
 392                 return b;
 393         if ((p5 = p5s) == 0) {
 394                 /* first time */
 395 #ifdef MULTIPLE_THREADS
 396                 ACQUIRE_DTOA_LOCK(1);
 397                 if (!(p5 = p5s)) {
 398                         p5 = p5s = i2b(625);
 399                         p5->next = 0;
 400                         }
 401                 FREE_DTOA_LOCK(1);
 402 #else
 403                 p5 = p5s = i2b(625);
 404                 p5->next = 0;
 405 #endif
 406                 }
 407         for(;;) {
 408                 if (k & 1) {
 409                         b1 = mult(b, p5);
 410                         Bfree(b);
 411                         b = b1;
 412                         }
 413                 if (!(k >>= 1))
 414                         break;
 415                 if ((p51 = p5->next) == 0) {
 416 #ifdef MULTIPLE_THREADS
 417                         ACQUIRE_DTOA_LOCK(1);
 418                         if (!(p51 = p5->next)) {
 419                                 p51 = p5->next = mult(p5,p5);
 420                                 p51->next = 0;
 421                                 }
 422                         FREE_DTOA_LOCK(1);
 423 #else
 424                         p51 = p5->next = mult(p5,p5);
 425                         p51->next = 0;
 426 #endif
 427                         }
 428                 p5 = p51;
 429                 }
 430         return b;
 431         }
 432
 433  Bigint *
 434 lshift
 435 #ifdef KR_headers
 436         (b, k) Bigint *b; int k;
 437 #else
 438         (Bigint *b, int k)
 439 #endif
 440 {
 441         int i, k1, n, n1;
 442         Bigint *b1;
 443         ULong *x, *x1, *xe, z;
 444
 445         n = k >> kshift;
 446         k1 = b->k;
 447         n1 = n + b->wds + 1;
 448         for(i = b->maxwds; n1 > i; i <<= 1)
 449                 k1++;
 450         b1 = Balloc(k1);
 451         x1 = b1->x;
 452         for(i = 0; i < n; i++)
 453                 *x1++ = 0;
 454         x = b->x;
 455         xe = x + b->wds;
 456         if (k &= kmask) {
 457 #ifdef Pack_32
 458                 k1 = 32 - k;
 459                 z = 0;
 460                 do {
 461                         *x1++ = *x << k | z;
 462                         z = *x++ >> k1;
 463                         }
 464                         while(x < xe);
 465                 if ((*x1 = z) !=0)
 466                         ++n1;
 467 #else
 468                 k1 = 16 - k;
 469                 z = 0;
 470                 do {
 471                         *x1++ = *x << k  & 0xffff | z;
 472                         z = *x++ >> k1;
 473                         }
 474                         while(x < xe);
 475                 if (*x1 = z)
 476                         ++n1;
 477 #endif
 478                 }
 479         else do
 480                 *x1++ = *x++;
 481                 while(x < xe);
 482         b1->wds = n1 - 1;
 483         Bfree(b);
 484         return b1;
 485         }
 486
 487  int
 488 cmp
 489 #ifdef KR_headers
 490         (a, b) Bigint *a, *b;
 491 #else
 492         (Bigint *a, Bigint *b)
 493 #endif
 494 {
 495         ULong *xa, *xa0, *xb, *xb0;
 496         int i, j;
 497
 498         i = a->wds;
 499         j = b->wds;
 500 #ifdef DEBUG
 501         if (i > 1 && !a->x[i-1])
 502                 Bug("cmp called with a->x[a->wds-1] == 0");
 503         if (j > 1 && !b->x[j-1])
 504                 Bug("cmp called with b->x[b->wds-1] == 0");
 505 #endif
 506         if (i -= j)
 507                 return i;
 508         xa0 = a->x;
 509         xa = xa0 + j;
 510         xb0 = b->x;
 511         xb = xb0 + j;
 512         for(;;) {
 513                 if (*--xa != *--xb)
 514                         return *xa < *xb ? -1 : 1;
 515                 if (xa <= xa0)
 516                         break;
 517                 }
 518         return 0;
 519         }
 520
 521  Bigint *
 522 diff
 523 #ifdef KR_headers
 524         (a, b) Bigint *a, *b;
 525 #else
 526         (Bigint *a, Bigint *b)
 527 #endif
 528 {
 529         Bigint *c;
 530         int i, wa, wb;
 531         ULong *xa, *xae, *xb, *xbe, *xc;
 532 #ifdef ULLong
 533         ULLong borrow, y;
 534 #else
 535         ULong borrow, y;
 536 #ifdef Pack_32
 537         ULong z;
 538 #endif
 539 #endif
 540
 541         i = cmp(a,b);
 542         if (!i) {
 543                 c = Balloc(0);
 544                 c->wds = 1;
 545                 c->x[0] = 0;
 546                 return c;
 547                 }
 548         if (i < 0) {
 549                 c = a;
 550                 a = b;
 551                 b = c;
 552                 i = 1;
 553                 }
 554         else
 555                 i = 0;
 556         c = Balloc(a->k);
 557         c->sign = i;
 558         wa = a->wds;
 559         xa = a->x;
 560         xae = xa + wa;
 561         wb = b->wds;
 562         xb = b->x;
 563         xbe = xb + wb;
 564         xc = c->x;
 565         borrow = 0;
 566 #ifdef ULLong
 567         do {
 568                 y = (ULLong)*xa++ - *xb++ - borrow;
 569                 borrow = y >> 32 & 1UL;
 570                 *xc++ = y & 0xffffffffUL;
 571                 }
 572                 while(xb < xbe);
 573         while(xa < xae) {
 574                 y = *xa++ - borrow;
 575                 borrow = y >> 32 & 1UL;
 576                 *xc++ = y & 0xffffffffUL;
 577                 }
 578 #else
 579 #ifdef Pack_32
 580         do {
 581                 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
 582                 borrow = (y & 0x10000) >> 16;
 583                 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
 584                 borrow = (z & 0x10000) >> 16;
 585                 Storeinc(xc, z, y);
 586                 }
 587                 while(xb < xbe);
 588         while(xa < xae) {
 589                 y = (*xa & 0xffff) - borrow;
 590                 borrow = (y & 0x10000) >> 16;
 591                 z = (*xa++ >> 16) - borrow;
 592                 borrow = (z & 0x10000) >> 16;
 593                 Storeinc(xc, z, y);
 594                 }
 595 #else
 596         do {
 597                 y = *xa++ - *xb++ - borrow;
 598                 borrow = (y & 0x10000) >> 16;
 599                 *xc++ = y & 0xffff;
 600                 }
 601                 while(xb < xbe);
 602         while(xa < xae) {
 603                 y = *xa++ - borrow;
 604                 borrow = (y & 0x10000) >> 16;
 605                 *xc++ = y & 0xffff;
 606                 }
 607 #endif
 608 #endif
 609         while(!*--xc)
 610                 wa--;
 611         c->wds = wa;
 612         return c;
 613         }
 614
 615  double
 616 b2d
 617 #ifdef KR_headers
 618         (a, e) Bigint *a; int *e;
 619 #else
 620         (Bigint *a, int *e)
 621 #endif
 622 {
 623         ULong *xa, *xa0, w, y, z;
 624         int k;
 625         U d;
 626 #ifdef VAX
 627         ULong d0, d1;
 628 #else
 629 #define d0 word0(&d)
 630 #define d1 word1(&d)
 631 #endif
 632
 633         xa0 = a->x;
 634         xa = xa0 + a->wds;
 635         y = *--xa;
 636 #ifdef DEBUG
 637         if (!y) Bug("zero y in b2d");
 638 #endif
 639         k = hi0bits(y);
 640         *e = 32 - k;
 641 #ifdef Pack_32
 642         if (k < Ebits) {
 643                 d0 = Exp_1 | y >> (Ebits - k);
 644                 w = xa > xa0 ? *--xa : 0;
 645                 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
 646                 goto ret_d;
 647                 }
 648         z = xa > xa0 ? *--xa : 0;
 649         if (k -= Ebits) {
 650                 d0 = Exp_1 | y << k | z >> (32 - k);
 651                 y = xa > xa0 ? *--xa : 0;
 652                 d1 = z << k | y >> (32 - k);
 653                 }
 654         else {
 655                 d0 = Exp_1 | y;
 656                 d1 = z;
 657                 }
 658 #else
 659         if (k < Ebits + 16) {
 660                 z = xa > xa0 ? *--xa : 0;
 661                 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
 662                 w = xa > xa0 ? *--xa : 0;
 663                 y = xa > xa0 ? *--xa : 0;
 664                 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
 665                 goto ret_d;
 666                 }
 667         z = xa > xa0 ? *--xa : 0;
 668         w = xa > xa0 ? *--xa : 0;
 669         k -= Ebits + 16;
 670         d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
 671         y = xa > xa0 ? *--xa : 0;
 672         d1 = w << k + 16 | y << k;
 673 #endif
 674  ret_d:
 675 #ifdef VAX
 676         word0(&d) = d0 >> 16 | d0 << 16;
 677         word1(&d) = d1 >> 16 | d1 << 16;
 678 #endif
 679         return dval(&d);
 680         }
 681 #undef d0
 682 #undef d1
 683
 684  Bigint *
 685 d2b
 686 #ifdef KR_headers
 687         (dd, e, bits) double dd; int *e, *bits;
 688 #else
 689         (double dd, int *e, int *bits)
 690 #endif
 691 {
 692         Bigint *b;
 693         U d;
 694 #ifndef Sudden_Underflow
 695         int i;
 696 #endif
 697         int de, k;
 698         ULong *x, y, z;
 699 #ifdef VAX
 700         ULong d0, d1;
 701 #else
 702 #define d0 word0(&d)
 703 #define d1 word1(&d)
 704 #endif
 705         d.d = dd;
 706 #ifdef VAX
 707         d0 = word0(&d) >> 16 | word0(&d) << 16;
 708         d1 = word1(&d) >> 16 | word1(&d) << 16;
 709 #endif
 710
 711 #ifdef Pack_32
 712         b = Balloc(1);
 713 #else
 714         b = Balloc(2);
 715 #endif
 716         x = b->x;
 717
 718         z = d0 & Frac_mask;
 719         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
 720 #ifdef Sudden_Underflow
 721         de = (int)(d0 >> Exp_shift);
 722 #ifndef IBM
 723         z |= Exp_msk11;
 724 #endif
 725 #else
 726         if ( (de = (int)(d0 >> Exp_shift)) !=0)
 727                 z |= Exp_msk1;
 728 #endif
 729 #ifdef Pack_32
 730         if ( (y = d1) !=0) {
 731                 if ( (k = lo0bits(&y)) !=0) {
 732                         x[0] = y | z << (32 - k);
 733                         z >>= k;
 734                         }
 735                 else
 736                         x[0] = y;
 737 #ifndef Sudden_Underflow
 738                 i =
 739 #endif
 740                      b->wds = (x[1] = z) !=0 ? 2 : 1;
 741                 }
 742         else {
 743                 k = lo0bits(&z);
 744                 x[0] = z;
 745 #ifndef Sudden_Underflow
 746                 i =
 747 #endif
 748                     b->wds = 1;
 749                 k += 32;
 750                 }
 751 #else
 752         if ( (y = d1) !=0) {
 753                 if ( (k = lo0bits(&y)) !=0)
 754                         if (k >= 16) {
 755                                 x[0] = y | z << 32 - k & 0xffff;
 756                                 x[1] = z >> k - 16 & 0xffff;
 757                                 x[2] = z >> k;
 758                                 i = 2;
 759                                 }
 760                         else {
 761                                 x[0] = y & 0xffff;
 762                                 x[1] = y >> 16 | z << 16 - k & 0xffff;
 763                                 x[2] = z >> k & 0xffff;
 764                                 x[3] = z >> k+16;
 765                                 i = 3;
 766                                 }
 767                 else {
 768                         x[0] = y & 0xffff;
 769                         x[1] = y >> 16;
 770                         x[2] = z & 0xffff;
 771                         x[3] = z >> 16;
 772                         i = 3;
 773                         }
 774                 }
 775         else {
 776 #ifdef DEBUG
 777                 if (!z)
 778                         Bug("Zero passed to d2b");
 779 #endif
 780                 k = lo0bits(&z);
 781                 if (k >= 16) {
 782                         x[0] = z;
 783                         i = 0;
 784                         }
 785                 else {
 786                         x[0] = z & 0xffff;
 787                         x[1] = z >> 16;
 788                         i = 1;
 789                         }
 790                 k += 32;
 791                 }
 792         while(!x[i])
 793                 --i;
 794         b->wds = i + 1;
 795 #endif
 796 #ifndef Sudden_Underflow
 797         if (de) {
 798 #endif
 799 #ifdef IBM
 800                 *e = (de - Bias - (P-1) << 2) + k;
 801                 *bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
 802 #else
 803                 *e = de - Bias - (P-1) + k;
 804                 *bits = P - k;
 805 #endif
 806 #ifndef Sudden_Underflow
 807                 }
 808         else {
 809                 *e = de - Bias - (P-1) + 1 + k;
 810 #ifdef Pack_32
 811                 *bits = 32*i - hi0bits(x[i-1]);
 812 #else
 813                 *bits = (i+2)*16 - hi0bits(x[i]);
 814 #endif
 815                 }
 816 #endif
 817         return b;
 818         }
 819 #undef d0
 820 #undef d1
 821
 822  CONST double
 823 #ifdef IEEE_Arith
 824 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
 825 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
 826                 };
 827 #else
 828 #ifdef IBM
 829 bigtens[] = { 1e16, 1e32, 1e64 };
 830 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
 831 #else
 832 bigtens[] = { 1e16, 1e32 };
 833 CONST double tinytens[] = { 1e-16, 1e-32 };
 834 #endif
 835 #endif
 836
 837  CONST double
 838 tens[] = {
 839                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
 840                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
 841                 1e20, 1e21, 1e22
 842 #ifdef VAX
 843                 , 1e23, 1e24
 844 #endif
 845                 };
 846
 847  char *
 848 #ifdef KR_headers
 849 strcp_D2A(a, b) char *a; char *b;
 850 #else
 851 strcp_D2A(char *a, CONST char *b)
 852 #endif
 853 {
 854         while((*a = *b++))
 855                 a++;
 856         return a;
 857         }
 858
 859 #ifdef NO_STRING_H
 860
 861  Char *
 862 #ifdef KR_headers
 863 memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
 864 #else
 865 memcpy_D2A(void *a1, void *b1, size_t len)
 866 #endif
 867 {
 868         char *a = (char*)a1, *ae = a + len;
 869         char *b = (char*)b1, *a0 = a;
 870         while(a < ae)
 871                 *a++ = *b++;
 872         return a0;
 873         }
 874
 875 #endif /* NO_STRING_H */