lib/libc/gdtoa/misc.c

   1 /* $NetBSD: misc.c,v 1.5 2009/01/30 23:35:35 lukem 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
  34 #include "gdtoaimp.h"
  35
  36  static Bigint *freelist[Kmax+1];
  37 #ifndef Omit_Private_Memory
  38 #ifndef PRIVATE_MEM
  39 #define PRIVATE_MEM 2304
  40 #endif
  41 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
  42 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
  43 #endif
  44
  45  Bigint *
  46 Balloc
  47 #ifdef KR_headers
  48         (k) int k;
  49 #else
  50         (int k)
  51 #endif
  52 {
  53         int x;
  54         Bigint *rv;
  55 #ifndef Omit_Private_Memory
  56         unsigned int len;
  57 #endif
  58
  59         ACQUIRE_DTOA_LOCK(0);
  60         if ( (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 ((double *)(pmem_next - private_mem + len) <= (double *)PRIVATE_mem) {
  71                         rv = (Bigint*)(void *)pmem_next;
  72                         pmem_next += len;
  73                         }
  74                 else
  75                         rv = (Bigint*)MALLOC(len*sizeof(double));
  76 #endif
  77                 if (rv == NULL)
  78                         return NULL;
  79                 rv->k = k;
  80                 rv->maxwds = x;
  81                 }
  82         FREE_DTOA_LOCK(0);
  83         rv->sign = rv->wds = 0;
  84         return rv;
  85         }
  86
  87  void
  88 Bfree
  89 #ifdef KR_headers
  90         (v) Bigint *v;
  91 #else
  92         (Bigint *v)
  93 #endif
  94 {
  95         if (v) {
  96                 ACQUIRE_DTOA_LOCK(0);
  97                 v->next = freelist[v->k];
  98                 freelist[v->k] = v;
  99                 FREE_DTOA_LOCK(0);
 100                 }
 101         }
 102
 103  int
 104 lo0bits
 105 #ifdef KR_headers
 106         (y) ULong *y;
 107 #else
 108         (ULong *y)
 109 #endif
 110 {
 111         int k;
 112         ULong x = *y;
 113
 114         if (x & 7) {
 115                 if (x & 1)
 116                         return 0;
 117                 if (x & 2) {
 118                         *y = x >> 1;
 119                         return 1;
 120                         }
 121                 *y = x >> 2;
 122                 return 2;
 123                 }
 124         k = 0;
 125         if (!(x & 0xffff)) {
 126                 k = 16;
 127                 x >>= 16;
 128                 }
 129         if (!(x & 0xff)) {
 130                 k += 8;
 131                 x >>= 8;
 132                 }
 133         if (!(x & 0xf)) {
 134                 k += 4;
 135                 x >>= 4;
 136                 }
 137         if (!(x & 0x3)) {
 138                 k += 2;
 139                 x >>= 2;
 140                 }
 141         if (!(x & 1)) {
 142                 k++;
 143                 x >>= 1;
 144                 if (!x)
 145                         return 32;
 146                 }
 147         *y = x;
 148         return k;
 149         }
 150
 151  Bigint *
 152 multadd
 153 #ifdef KR_headers
 154         (b, m, a) Bigint *b; int m, a;
 155 #else
 156         (Bigint *b, int m, int a)       /* multiply by m and add a */
 157 #endif
 158 {
 159         int i, wds;
 160 #ifdef ULLong
 161         ULong *x;
 162         ULLong carry, y;
 163 #else
 164         ULong carry, *x, y;
 165 #ifdef Pack_32
 166         ULong xi, z;
 167 #endif
 168 #endif
 169         Bigint *b1;
 170
 171         wds = b->wds;
 172         x = b->x;
 173         i = 0;
 174         carry = a;
 175         do {
 176 #ifdef ULLong
 177                 y = *x * (ULLong)m + carry;
 178                 carry = y >> 32;
 179                 /* LINTED conversion */
 180                 *x++ = y & 0xffffffffUL;
 181 #else
 182 #ifdef Pack_32
 183                 xi = *x;
 184                 y = (xi & 0xffff) * m + carry;
 185                 z = (xi >> 16) * m + (y >> 16);
 186                 carry = z >> 16;
 187                 *x++ = (z << 16) + (y & 0xffff);
 188 #else
 189                 y = *x * m + carry;
 190                 carry = y >> 16;
 191                 *x++ = y & 0xffff;
 192 #endif
 193 #endif
 194                 }
 195                 while(++i < wds);
 196         if (carry) {
 197                 if (wds >= b->maxwds) {
 198                         b1 = Balloc(b->k+1);
 199                         if (b1 == NULL) {
 200                                 Bfree(b);
 201                                 return NULL;
 202                                 }
 203                         Bcopy(b1, b);
 204                         Bfree(b);
 205                         b = b1;
 206                         }
 207                 /* LINTED conversion */
 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         if (b == NULL)
 260                 return NULL;
 261         b->x[0] = i;
 262         b->wds = 1;
 263         return b;
 264         }
 265
 266  Bigint *
 267 mult
 268 #ifdef KR_headers
 269         (a, b) Bigint *a, *b;
 270 #else
 271         (Bigint *a, Bigint *b)
 272 #endif
 273 {
 274         Bigint *c;
 275         int k, wa, wb, wc;
 276         ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
 277         ULong y;
 278 #ifdef ULLong
 279         ULLong carry, z;
 280 #else
 281         ULong carry, z;
 282 #ifdef Pack_32
 283         ULong z2;
 284 #endif
 285 #endif
 286
 287         if (a->wds < b->wds) {
 288                 c = a;
 289                 a = b;
 290                 b = c;
 291                 }
 292         k = a->k;
 293         wa = a->wds;
 294         wb = b->wds;
 295         wc = wa + wb;
 296         if (wc > a->maxwds)
 297                 k++;
 298         c = Balloc(k);
 299         if (c == NULL)
 300                 return NULL;
 301         for(x = c->x, xa = x + wc; x < xa; x++)
 302                 *x = 0;
 303         xa = a->x;
 304         xae = xa + wa;
 305         xb = b->x;
 306         xbe = xb + wb;
 307         xc0 = c->x;
 308 #ifdef ULLong
 309         for(; xb < xbe; xc0++) {
 310                 if ( (y = *xb++) !=0) {
 311                         x = xa;
 312                         xc = xc0;
 313                         carry = 0;
 314                         do {
 315                                 z = *x++ * (ULLong)y + *xc + carry;
 316                                 carry = z >> 32;
 317                                 /* LINTED conversion */
 318                                 *xc++ = z & 0xffffffffUL;
 319                                 }
 320                                 while(x < xae);
 321                         /* LINTED conversion */
 322                         *xc = carry;
 323                         }
 324                 }
 325 #else
 326 #ifdef Pack_32
 327         for(; xb < xbe; xb++, xc0++) {
 328                 if ( (y = *xb & 0xffff) !=0) {
 329                         x = xa;
 330                         xc = xc0;
 331                         carry = 0;
 332                         do {
 333                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
 334                                 carry = z >> 16;
 335                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
 336                                 carry = z2 >> 16;
 337                                 Storeinc(xc, z2, z);
 338                                 }
 339                                 while(x < xae);
 340                         *xc = carry;
 341                         }
 342                 if ( (y = *xb >> 16) !=0) {
 343                         x = xa;
 344                         xc = xc0;
 345                         carry = 0;
 346                         z2 = *xc;
 347                         do {
 348                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
 349                                 carry = z >> 16;
 350                                 Storeinc(xc, z, z2);
 351                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
 352                                 carry = z2 >> 16;
 353                                 }
 354                                 while(x < xae);
 355                         *xc = z2;
 356                         }
 357                 }
 358 #else
 359         for(; xb < xbe; xc0++) {
 360                 if ( (y = *xb++) !=0) {
 361                         x = xa;
 362                         xc = xc0;
 363                         carry = 0;
 364                         do {
 365                                 z = *x++ * y + *xc + carry;
 366                                 carry = z >> 16;
 367                                 *xc++ = z & 0xffff;
 368                                 }
 369                                 while(x < xae);
 370                         *xc = carry;
 371                         }
 372                 }
 373 #endif
 374 #endif
 375         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
 376         c->wds = wc;
 377         return c;
 378         }
 379
 380  static Bigint *p5s;
 381
 382  Bigint *
 383 pow5mult
 384 #ifdef KR_headers
 385         (b, k) Bigint *b; int k;
 386 #else
 387         (Bigint *b, int k)
 388 #endif
 389 {
 390         Bigint *b1, *p5, *p51;
 391         int i;
 392         static CONST int p05[3] = { 5, 25, 125 };
 393
 394         if ( (i = k & 3) !=0) {
 395                 b = multadd(b, p05[i-1], 0);
 396                 if (b == NULL)
 397                         return NULL;
 398                 }
 399
 400         if (!(k = (unsigned int)k >> 2))
 401                 return b;
 402         if ((p5 = p5s) == 0) {
 403                 /* first time */
 404 #ifdef MULTIPLE_THREADS
 405                 ACQUIRE_DTOA_LOCK(1);
 406                 if (!(p5 = p5s)) {
 407                         p5 = p5s = i2b(625);
 408                         if (p5 == NULL)
 409                                 return NULL;
 410                         p5->next = 0;
 411                         }
 412                 FREE_DTOA_LOCK(1);
 413 #else
 414                 p5 = p5s = i2b(625);
 415                 if (p5 == NULL)
 416                         return NULL;
 417                 p5->next = 0;
 418 #endif
 419                 }
 420         for(;;) {
 421                 if (k & 1) {
 422                         b1 = mult(b, p5);
 423                         if (b1 == NULL)
 424                                 return NULL;
 425                         b = b1;
 426                         }
 427                 if (!(k = (unsigned int)k >> 1))
 428                         break;
 429                 if ((p51 = p5->next) == 0) {
 430 #ifdef MULTIPLE_THREADS
 431                         ACQUIRE_DTOA_LOCK(1);
 432                         if (!(p51 = p5->next)) {
 433                                 p51 = p5->next = mult(p5,p5);
 434                                 if (p51 == NULL)
 435                                         return NULL;
 436                                 p51->next = 0;
 437                                 }
 438                         FREE_DTOA_LOCK(1);
 439 #else
 440                         p51 = p5->next = mult(p5,p5);
 441                         if (p51 == NULL)
 442                                 return NULL;
 443                         p51->next = 0;
 444 #endif
 445                         }
 446                 p5 = p51;
 447                 }
 448         return b;
 449         }
 450
 451  Bigint *
 452 lshift
 453 #ifdef KR_headers
 454         (b, k) Bigint *b; int k;
 455 #else
 456         (Bigint *b, int k)
 457 #endif
 458 {
 459         int i, k1, n, n1;
 460         Bigint *b1;
 461         ULong *x, *x1, *xe, z;
 462
 463         n = (unsigned int)k >> kshift;
 464         k1 = b->k;
 465         n1 = n + b->wds + 1;
 466         for(i = b->maxwds; n1 > i; i <<= 1)
 467                 k1++;
 468         b1 = Balloc(k1);
 469         if (b1 == NULL)
 470                 return NULL;
 471         x1 = b1->x;
 472         for(i = 0; i < n; i++)
 473                 *x1++ = 0;
 474         x = b->x;
 475         xe = x + b->wds;
 476         if (k &= kmask) {
 477 #ifdef Pack_32
 478                 k1 = 32 - k;
 479                 z = 0;
 480                 do {
 481                         *x1++ = *x << k | z;
 482                         z = *x++ >> k1;
 483                         }
 484                         while(x < xe);
 485                 if ((*x1 = z) !=0)
 486                         ++n1;
 487 #else
 488                 k1 = 16 - k;
 489                 z = 0;
 490                 do {
 491                         *x1++ = *x << k  & 0xffff | z;
 492                         z = *x++ >> k1;
 493                         }
 494                         while(x < xe);
 495                 if (*x1 = z)
 496                         ++n1;
 497 #endif
 498                 }
 499         else do
 500                 *x1++ = *x++;
 501                 while(x < xe);
 502         b1->wds = n1 - 1;
 503         Bfree(b);
 504         return b1;
 505         }
 506
 507  int
 508 cmp
 509 #ifdef KR_headers
 510         (a, b) Bigint *a, *b;
 511 #else
 512         (Bigint *a, Bigint *b)
 513 #endif
 514 {
 515         ULong *xa, *xa0, *xb, *xb0;
 516         int i, j;
 517
 518         i = a->wds;
 519         j = b->wds;
 520 #ifdef DEBUG
 521         if (i > 1 && !a->x[i-1])
 522                 Bug("cmp called with a->x[a->wds-1] == 0");
 523         if (j > 1 && !b->x[j-1])
 524                 Bug("cmp called with b->x[b->wds-1] == 0");
 525 #endif
 526         if (i -= j)
 527                 return i;
 528         xa0 = a->x;
 529         xa = xa0 + j;
 530         xb0 = b->x;
 531         xb = xb0 + j;
 532         for(;;) {
 533                 if (*--xa != *--xb)
 534                         return *xa < *xb ? -1 : 1;
 535                 if (xa <= xa0)
 536                         break;
 537                 }
 538         return 0;
 539         }
 540
 541  Bigint *
 542 diff
 543 #ifdef KR_headers
 544         (a, b) Bigint *a, *b;
 545 #else
 546         (Bigint *a, Bigint *b)
 547 #endif
 548 {
 549         Bigint *c;
 550         int i, wa, wb;
 551         ULong *xa, *xae, *xb, *xbe, *xc;
 552 #ifdef ULLong
 553         ULLong borrow, y;
 554 #else
 555         ULong borrow, y;
 556 #ifdef Pack_32
 557         ULong z;
 558 #endif
 559 #endif
 560
 561         i = cmp(a,b);
 562         if (!i) {
 563                 c = Balloc(0);
 564                 if (c == NULL)
 565                         return NULL;
 566                 c->wds = 1;
 567                 c->x[0] = 0;
 568                 return c;
 569                 }
 570         if (i < 0) {
 571                 c = a;
 572                 a = b;
 573                 b = c;
 574                 i = 1;
 575                 }
 576         else
 577                 i = 0;
 578         c = Balloc(a->k);
 579         if (c == NULL)
 580                 return NULL;
 581         c->sign = i;
 582         wa = a->wds;
 583         xa = a->x;
 584         xae = xa + wa;
 585         wb = b->wds;
 586         xb = b->x;
 587         xbe = xb + wb;
 588         xc = c->x;
 589         borrow = 0;
 590 #ifdef ULLong
 591         do {
 592                 y = (ULLong)*xa++ - *xb++ - borrow;
 593                 borrow = y >> 32 & 1UL;
 594                 /* LINTED conversion */
 595                 *xc++ = y & 0xffffffffUL;
 596                 }
 597                 while(xb < xbe);
 598         while(xa < xae) {
 599                 y = *xa++ - borrow;
 600                 borrow = y >> 32 & 1UL;
 601                 /* LINTED conversion */
 602                 *xc++ = y & 0xffffffffUL;
 603                 }
 604 #else
 605 #ifdef Pack_32
 606         do {
 607                 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
 608                 borrow = (y & 0x10000) >> 16;
 609                 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
 610                 borrow = (z & 0x10000) >> 16;
 611                 Storeinc(xc, z, y);
 612                 }
 613                 while(xb < xbe);
 614         while(xa < xae) {
 615                 y = (*xa & 0xffff) - borrow;
 616                 borrow = (y & 0x10000) >> 16;
 617                 z = (*xa++ >> 16) - borrow;
 618                 borrow = (z & 0x10000) >> 16;
 619                 Storeinc(xc, z, y);
 620                 }
 621 #else
 622         do {
 623                 y = *xa++ - *xb++ - borrow;
 624                 borrow = (y & 0x10000) >> 16;
 625                 *xc++ = y & 0xffff;
 626                 }
 627                 while(xb < xbe);
 628         while(xa < xae) {
 629                 y = *xa++ - borrow;
 630                 borrow = (y & 0x10000) >> 16;
 631                 *xc++ = y & 0xffff;
 632                 }
 633 #endif
 634 #endif
 635         while(!*--xc)
 636                 wa--;
 637         c->wds = wa;
 638         return c;
 639         }
 640
 641  double
 642 b2d
 643 #ifdef KR_headers
 644         (a, e) Bigint *a; int *e;
 645 #else
 646         (Bigint *a, int *e)
 647 #endif
 648 {
 649         ULong *xa, *xa0, w, y, z;
 650         int k;
 651         double d;
 652 #ifdef VAX
 653         ULong d0, d1;
 654 #else
 655 #define d0 word0(d)
 656 #define d1 word1(d)
 657 #endif
 658
 659         xa0 = a->x;
 660         xa = xa0 + a->wds;
 661         y = *--xa;
 662 #ifdef DEBUG
 663         if (!y) Bug("zero y in b2d");
 664 #endif
 665         k = hi0bits(y);
 666         *e = 32 - k;
 667 #ifdef Pack_32
 668         if (k < Ebits) {
 669                 d0 = Exp_1 | y >> (Ebits - k);
 670                 w = xa > xa0 ? *--xa : 0;
 671                 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
 672                 goto ret_d;
 673                 }
 674         z = xa > xa0 ? *--xa : 0;
 675         if (k -= Ebits) {
 676                 d0 = Exp_1 | y << k | z >> (32 - k);
 677                 y = xa > xa0 ? *--xa : 0;
 678                 d1 = z << k | y >> (32 - k);
 679                 }
 680         else {
 681                 d0 = Exp_1 | y;
 682                 d1 = z;
 683                 }
 684 #else
 685         if (k < Ebits + 16) {
 686                 z = xa > xa0 ? *--xa : 0;
 687                 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
 688                 w = xa > xa0 ? *--xa : 0;
 689                 y = xa > xa0 ? *--xa : 0;
 690                 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
 691                 goto ret_d;
 692                 }
 693         z = xa > xa0 ? *--xa : 0;
 694         w = xa > xa0 ? *--xa : 0;
 695         k -= Ebits + 16;
 696         d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
 697         y = xa > xa0 ? *--xa : 0;
 698         d1 = w << k + 16 | y << k;
 699 #endif
 700  ret_d:
 701 #ifdef VAX
 702         word0(d) = d0 >> 16 | d0 << 16;
 703         word1(d) = d1 >> 16 | d1 << 16;
 704 #endif
 705         return dval(d);
 706         }
 707 #undef d0
 708 #undef d1
 709
 710  Bigint *
 711 d2b
 712 #ifdef KR_headers
 713         (d, e, bits) double d; int *e, *bits;
 714 #else
 715         (double d, int *e, int *bits)
 716 #endif
 717 {
 718         Bigint *b;
 719 #ifndef Sudden_Underflow
 720         int i;
 721 #endif
 722         int de, k;
 723         ULong *x, y, z;
 724 #ifdef VAX
 725         ULong d0, d1;
 726         d0 = word0(d) >> 16 | word0(d) << 16;
 727         d1 = word1(d) >> 16 | word1(d) << 16;
 728 #else
 729 #define d0 word0(d)
 730 #define d1 word1(d)
 731 #endif
 732
 733 #ifdef Pack_32
 734         b = Balloc(1);
 735 #else
 736         b = Balloc(2);
 737 #endif
 738         if (b == NULL)
 739                 return NULL;
 740         x = b->x;
 741
 742         z = d0 & Frac_mask;
 743         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
 744 #ifdef Sudden_Underflow
 745         de = (int)(d0 >> Exp_shift);
 746 #ifndef IBM
 747         z |= Exp_msk11;
 748 #endif
 749 #else
 750         if ( (de = (int)(d0 >> Exp_shift)) !=0)
 751                 z |= Exp_msk1;
 752 #endif
 753 #ifdef Pack_32
 754         if ( (y = d1) !=0) {
 755                 if ( (k = lo0bits(&y)) !=0) {
 756                         x[0] = y | z << (32 - k);
 757                         z >>= k;
 758                         }
 759                 else
 760                         x[0] = y;
 761 #ifndef Sudden_Underflow
 762                 i =
 763 #endif
 764                      b->wds = (x[1] = z) !=0 ? 2 : 1;
 765                 }
 766         else {
 767 #ifdef DEBUG
 768                 if (!z)
 769                         Bug("Zero passed to d2b");
 770 #endif
 771                 k = lo0bits(&z);
 772                 x[0] = z;
 773 #ifndef Sudden_Underflow
 774                 i =
 775 #endif
 776                     b->wds = 1;
 777                 k += 32;
 778                 }
 779 #else
 780         if ( (y = d1) !=0) {
 781                 if ( (k = lo0bits(&y)) !=0)
 782                         if (k >= 16) {
 783                                 x[0] = y | z << 32 - k & 0xffff;
 784                                 x[1] = z >> k - 16 & 0xffff;
 785                                 x[2] = z >> k;
 786                                 i = 2;
 787                                 }
 788                         else {
 789                                 x[0] = y & 0xffff;
 790                                 x[1] = y >> 16 | z << 16 - k & 0xffff;
 791                                 x[2] = z >> k & 0xffff;
 792                                 x[3] = z >> k+16;
 793                                 i = 3;
 794                                 }
 795                 else {
 796                         x[0] = y & 0xffff;
 797                         x[1] = y >> 16;
 798                         x[2] = z & 0xffff;
 799                         x[3] = z >> 16;
 800                         i = 3;
 801                         }
 802                 }
 803         else {
 804 #ifdef DEBUG
 805                 if (!z)
 806                         Bug("Zero passed to d2b");
 807 #endif
 808                 k = lo0bits(&z);
 809                 if (k >= 16) {
 810                         x[0] = z;
 811                         i = 0;
 812                         }
 813                 else {
 814                         x[0] = z & 0xffff;
 815                         x[1] = z >> 16;
 816                         i = 1;
 817                         }
 818                 k += 32;
 819                 }
 820         while(!x[i])
 821                 --i;
 822         b->wds = i + 1;
 823 #endif
 824 #ifndef Sudden_Underflow
 825         if (de) {
 826 #endif
 827 #ifdef IBM
 828                 *e = (de - Bias - (P-1) << 2) + k;
 829                 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
 830 #else
 831                 *e = de - Bias - (P-1) + k;
 832                 *bits = P - k;
 833 #endif
 834 #ifndef Sudden_Underflow
 835                 }
 836         else {
 837                 *e = de - Bias - (P-1) + 1 + k;
 838 #ifdef Pack_32
 839                 *bits = 32*i - hi0bits(x[i-1]);
 840 #else
 841                 *bits = (i+2)*16 - hi0bits(x[i]);
 842 #endif
 843                 }
 844 #endif
 845         return b;
 846         }
 847 #undef d0
 848 #undef d1
 849
 850  CONST double
 851 #ifdef IEEE_Arith
 852 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
 853 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
 854                 };
 855 #else
 856 #ifdef IBM
 857 bigtens[] = { 1e16, 1e32, 1e64 };
 858 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
 859 #else
 860 bigtens[] = { 1e16, 1e32 };
 861 CONST double tinytens[] = { 1e-16, 1e-32 };
 862 #endif
 863 #endif
 864
 865  CONST double
 866 tens[] = {
 867                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
 868                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
 869                 1e20, 1e21, 1e22
 870 #ifdef VAX
 871                 , 1e23, 1e24
 872 #endif
 873                 };
 874
 875  char *
 876 #ifdef KR_headers
 877 strcp_D2A(a, b) char *a; char *b;
 878 #else
 879 strcp_D2A(char *a, CONST char *b)
 880 #endif
 881 {
 882         while((*a = *b++))
 883                 a++;
 884         return a;
 885         }
 886
 887 #ifdef NO_STRING_H
 888
 889  Char *
 890 #ifdef KR_headers
 891 memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
 892 #else
 893 memcpy_D2A(void *a1, void *b1, size_t len)
 894 #endif
 895 {
 896         char *a = (char*)a1, *ae = a + len;
 897         char *b = (char*)b1, *a0 = a;
 898         while(a < ae)
 899                 *a++ = *b++;
 900         return a0;
 901         }
 902
 903 #endif /* NO_STRING_H */