StdLib/LibC/gdtoa/misc.c

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