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