Home | History | Annotate | Download | only in gdtoa
      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 		if (rv == NULL)
     68 			return (NULL);
     69 #else
     70 		len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
     71 			/sizeof(double);
     72 		if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
     73 			rv = (Bigint*)pmem_next;
     74 			pmem_next += len;
     75 			}
     76 		else {
     77 			rv = (Bigint*)MALLOC(len*sizeof(double));
     78 			if (rv == NULL)
     79 				return (NULL);
     80 		}
     81 #endif
     82 		rv->k = k;
     83 		rv->maxwds = x;
     84 		}
     85 	FREE_DTOA_LOCK(0);
     86 	rv->sign = rv->wds = 0;
     87 	return rv;
     88 	}
     89 
     90  void
     91 Bfree
     92 #ifdef KR_headers
     93 	(v) Bigint *v;
     94 #else
     95 	(Bigint *v)
     96 #endif
     97 {
     98 	if (v) {
     99 		if (v->k > Kmax)
    100 #ifdef FREE
    101 			FREE((void*)v);
    102 #else
    103 			free((void*)v);
    104 #endif
    105 		else {
    106 			ACQUIRE_DTOA_LOCK(0);
    107 			v->next = freelist[v->k];
    108 			freelist[v->k] = v;
    109 			FREE_DTOA_LOCK(0);
    110 			}
    111 		}
    112 	}
    113 
    114  int
    115 lo0bits
    116 #ifdef KR_headers
    117 	(y) ULong *y;
    118 #else
    119 	(ULong *y)
    120 #endif
    121 {
    122 	int k;
    123 	ULong x = *y;
    124 
    125 	if (x & 7) {
    126 		if (x & 1)
    127 			return 0;
    128 		if (x & 2) {
    129 			*y = x >> 1;
    130 			return 1;
    131 			}
    132 		*y = x >> 2;
    133 		return 2;
    134 		}
    135 	k = 0;
    136 	if (!(x & 0xffff)) {
    137 		k = 16;
    138 		x >>= 16;
    139 		}
    140 	if (!(x & 0xff)) {
    141 		k += 8;
    142 		x >>= 8;
    143 		}
    144 	if (!(x & 0xf)) {
    145 		k += 4;
    146 		x >>= 4;
    147 		}
    148 	if (!(x & 0x3)) {
    149 		k += 2;
    150 		x >>= 2;
    151 		}
    152 	if (!(x & 1)) {
    153 		k++;
    154 		x >>= 1;
    155 		if (!x)
    156 			return 32;
    157 		}
    158 	*y = x;
    159 	return k;
    160 	}
    161 
    162  Bigint *
    163 multadd
    164 #ifdef KR_headers
    165 	(b, m, a) Bigint *b; int m, a;
    166 #else
    167 	(Bigint *b, int m, int a)	/* multiply by m and add a */
    168 #endif
    169 {
    170 	int i, wds;
    171 #ifdef ULLong
    172 	ULong *x;
    173 	ULLong carry, y;
    174 #else
    175 	ULong carry, *x, y;
    176 #ifdef Pack_32
    177 	ULong xi, z;
    178 #endif
    179 #endif
    180 	Bigint *b1;
    181 
    182 	wds = b->wds;
    183 	x = b->x;
    184 	i = 0;
    185 	carry = a;
    186 	do {
    187 #ifdef ULLong
    188 		y = *x * (ULLong)m + carry;
    189 		carry = y >> 32;
    190 		*x++ = y & 0xffffffffUL;
    191 #else
    192 #ifdef Pack_32
    193 		xi = *x;
    194 		y = (xi & 0xffff) * m + carry;
    195 		z = (xi >> 16) * m + (y >> 16);
    196 		carry = z >> 16;
    197 		*x++ = (z << 16) + (y & 0xffff);
    198 #else
    199 		y = *x * m + carry;
    200 		carry = y >> 16;
    201 		*x++ = y & 0xffff;
    202 #endif
    203 #endif
    204 		}
    205 		while(++i < wds);
    206 	if (carry) {
    207 		if (wds >= b->maxwds) {
    208 			b1 = Balloc(b->k+1);
    209 			if (b1 == NULL)
    210 				return (NULL);
    211 			Bcopy(b1, b);
    212 			Bfree(b);
    213 			b = b1;
    214 			}
    215 		b->x[wds++] = carry;
    216 		b->wds = wds;
    217 		}
    218 	return b;
    219 	}
    220 
    221  int
    222 hi0bits_D2A
    223 #ifdef KR_headers
    224 	(x) ULong x;
    225 #else
    226 	(ULong x)
    227 #endif
    228 {
    229 	int k = 0;
    230 
    231 	if (!(x & 0xffff0000)) {
    232 		k = 16;
    233 		x <<= 16;
    234 		}
    235 	if (!(x & 0xff000000)) {
    236 		k += 8;
    237 		x <<= 8;
    238 		}
    239 	if (!(x & 0xf0000000)) {
    240 		k += 4;
    241 		x <<= 4;
    242 		}
    243 	if (!(x & 0xc0000000)) {
    244 		k += 2;
    245 		x <<= 2;
    246 		}
    247 	if (!(x & 0x80000000)) {
    248 		k++;
    249 		if (!(x & 0x40000000))
    250 			return 32;
    251 		}
    252 	return k;
    253 	}
    254 
    255  Bigint *
    256 i2b
    257 #ifdef KR_headers
    258 	(i) int i;
    259 #else
    260 	(int i)
    261 #endif
    262 {
    263 	Bigint *b;
    264 
    265 	b = Balloc(1);
    266 	if (b == NULL)
    267 		return (NULL);
    268 	b->x[0] = i;
    269 	b->wds = 1;
    270 	return b;
    271 	}
    272 
    273  Bigint *
    274 mult
    275 #ifdef KR_headers
    276 	(a, b) Bigint *a, *b;
    277 #else
    278 	(Bigint *a, Bigint *b)
    279 #endif
    280 {
    281 	Bigint *c;
    282 	int k, wa, wb, wc;
    283 	ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
    284 	ULong y;
    285 #ifdef ULLong
    286 	ULLong carry, z;
    287 #else
    288 	ULong carry, z;
    289 #ifdef Pack_32
    290 	ULong z2;
    291 #endif
    292 #endif
    293 
    294 	if (a->wds < b->wds) {
    295 		c = a;
    296 		a = b;
    297 		b = c;
    298 		}
    299 	k = a->k;
    300 	wa = a->wds;
    301 	wb = b->wds;
    302 	wc = wa + wb;
    303 	if (wc > a->maxwds)
    304 		k++;
    305 	c = Balloc(k);
    306 	if (c == NULL)
    307 		return (NULL);
    308 	for(x = c->x, xa = x + wc; x < xa; x++)
    309 		*x = 0;
    310 	xa = a->x;
    311 	xae = xa + wa;
    312 	xb = b->x;
    313 	xbe = xb + wb;
    314 	xc0 = c->x;
    315 #ifdef ULLong
    316 	for(; xb < xbe; xc0++) {
    317 		if ( (y = *xb++) !=0) {
    318 			x = xa;
    319 			xc = xc0;
    320 			carry = 0;
    321 			do {
    322 				z = *x++ * (ULLong)y + *xc + carry;
    323 				carry = z >> 32;
    324 				*xc++ = z & 0xffffffffUL;
    325 				}
    326 				while(x < xae);
    327 			*xc = carry;
    328 			}
    329 		}
    330 #else
    331 #ifdef Pack_32
    332 	for(; xb < xbe; xb++, xc0++) {
    333 		if ( (y = *xb & 0xffff) !=0) {
    334 			x = xa;
    335 			xc = xc0;
    336 			carry = 0;
    337 			do {
    338 				z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
    339 				carry = z >> 16;
    340 				z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
    341 				carry = z2 >> 16;
    342 				Storeinc(xc, z2, z);
    343 				}
    344 				while(x < xae);
    345 			*xc = carry;
    346 			}
    347 		if ( (y = *xb >> 16) !=0) {
    348 			x = xa;
    349 			xc = xc0;
    350 			carry = 0;
    351 			z2 = *xc;
    352 			do {
    353 				z = (*x & 0xffff) * y + (*xc >> 16) + carry;
    354 				carry = z >> 16;
    355 				Storeinc(xc, z, z2);
    356 				z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
    357 				carry = z2 >> 16;
    358 				}
    359 				while(x < xae);
    360 			*xc = z2;
    361 			}
    362 		}
    363 #else
    364 	for(; xb < xbe; xc0++) {
    365 		if ( (y = *xb++) !=0) {
    366 			x = xa;
    367 			xc = xc0;
    368 			carry = 0;
    369 			do {
    370 				z = *x++ * y + *xc + carry;
    371 				carry = z >> 16;
    372 				*xc++ = z & 0xffff;
    373 				}
    374 				while(x < xae);
    375 			*xc = carry;
    376 			}
    377 		}
    378 #endif
    379 #endif
    380 	for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
    381 	c->wds = wc;
    382 	return c;
    383 	}
    384 
    385  static Bigint *p5s;
    386 
    387  Bigint *
    388 pow5mult
    389 #ifdef KR_headers
    390 	(b, k) Bigint *b; int k;
    391 #else
    392 	(Bigint *b, int k)
    393 #endif
    394 {
    395 	Bigint *b1, *p5, *p51;
    396 	int i;
    397 	static int p05[3] = { 5, 25, 125 };
    398 
    399 	if ( (i = k & 3) !=0) {
    400 		b = multadd(b, p05[i-1], 0);
    401 		if (b == NULL)
    402 			return (NULL);
    403 		}
    404 
    405 	if (!(k >>= 2))
    406 		return b;
    407 	if ((p5 = p5s) == 0) {
    408 		/* first time */
    409 #ifdef MULTIPLE_THREADS
    410 		ACQUIRE_DTOA_LOCK(1);
    411 		if (!(p5 = p5s)) {
    412 			p5 = p5s = i2b(625);
    413 			if (p5 == NULL)
    414 				return (NULL);
    415 			p5->next = 0;
    416 			}
    417 		FREE_DTOA_LOCK(1);
    418 #else
    419 		p5 = p5s = i2b(625);
    420 		if (p5 == NULL)
    421 			return (NULL);
    422 		p5->next = 0;
    423 #endif
    424 		}
    425 	for(;;) {
    426 		if (k & 1) {
    427 			b1 = mult(b, p5);
    428 			if (b1 == NULL)
    429 				return (NULL);
    430 			Bfree(b);
    431 			b = b1;
    432 			}
    433 		if (!(k >>= 1))
    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 = 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 		*xc++ = y & 0xffffffffUL;
    601 		}
    602 		while(xb < xbe);
    603 	while(xa < xae) {
    604 		y = *xa++ - borrow;
    605 		borrow = y >> 32 & 1UL;
    606 		*xc++ = y & 0xffffffffUL;
    607 		}
    608 #else
    609 #ifdef Pack_32
    610 	do {
    611 		y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
    612 		borrow = (y & 0x10000) >> 16;
    613 		z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
    614 		borrow = (z & 0x10000) >> 16;
    615 		Storeinc(xc, z, y);
    616 		}
    617 		while(xb < xbe);
    618 	while(xa < xae) {
    619 		y = (*xa & 0xffff) - borrow;
    620 		borrow = (y & 0x10000) >> 16;
    621 		z = (*xa++ >> 16) - borrow;
    622 		borrow = (z & 0x10000) >> 16;
    623 		Storeinc(xc, z, y);
    624 		}
    625 #else
    626 	do {
    627 		y = *xa++ - *xb++ - borrow;
    628 		borrow = (y & 0x10000) >> 16;
    629 		*xc++ = y & 0xffff;
    630 		}
    631 		while(xb < xbe);
    632 	while(xa < xae) {
    633 		y = *xa++ - borrow;
    634 		borrow = (y & 0x10000) >> 16;
    635 		*xc++ = y & 0xffff;
    636 		}
    637 #endif
    638 #endif
    639 	while(!*--xc)
    640 		wa--;
    641 	c->wds = wa;
    642 	return c;
    643 	}
    644 
    645  double
    646 b2d
    647 #ifdef KR_headers
    648 	(a, e) Bigint *a; int *e;
    649 #else
    650 	(Bigint *a, int *e)
    651 #endif
    652 {
    653 	ULong *xa, *xa0, w, y, z;
    654 	int k;
    655 	U d;
    656 #ifdef VAX
    657 	ULong d0, d1;
    658 #else
    659 #define d0 word0(&d)
    660 #define d1 word1(&d)
    661 #endif
    662 
    663 	xa0 = a->x;
    664 	xa = xa0 + a->wds;
    665 	y = *--xa;
    666 #ifdef DEBUG
    667 	if (!y) Bug("zero y in b2d");
    668 #endif
    669 	k = hi0bits(y);
    670 	*e = 32 - k;
    671 #ifdef Pack_32
    672 	if (k < Ebits) {
    673 		d0 = Exp_1 | y >> (Ebits - k);
    674 		w = xa > xa0 ? *--xa : 0;
    675 		d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
    676 		goto ret_d;
    677 		}
    678 	z = xa > xa0 ? *--xa : 0;
    679 	if (k -= Ebits) {
    680 		d0 = Exp_1 | y << k | z >> (32 - k);
    681 		y = xa > xa0 ? *--xa : 0;
    682 		d1 = z << k | y >> (32 - k);
    683 		}
    684 	else {
    685 		d0 = Exp_1 | y;
    686 		d1 = z;
    687 		}
    688 #else
    689 	if (k < Ebits + 16) {
    690 		z = xa > xa0 ? *--xa : 0;
    691 		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
    692 		w = xa > xa0 ? *--xa : 0;
    693 		y = xa > xa0 ? *--xa : 0;
    694 		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
    695 		goto ret_d;
    696 		}
    697 	z = xa > xa0 ? *--xa : 0;
    698 	w = xa > xa0 ? *--xa : 0;
    699 	k -= Ebits + 16;
    700 	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
    701 	y = xa > xa0 ? *--xa : 0;
    702 	d1 = w << k + 16 | y << k;
    703 #endif
    704  ret_d:
    705 #ifdef VAX
    706 	word0(&d) = d0 >> 16 | d0 << 16;
    707 	word1(&d) = d1 >> 16 | d1 << 16;
    708 #endif
    709 	return dval(&d);
    710 	}
    711 #undef d0
    712 #undef d1
    713 
    714  Bigint *
    715 d2b
    716 #ifdef KR_headers
    717 	(dd, e, bits) double dd; int *e, *bits;
    718 #else
    719 	(double dd, int *e, int *bits)
    720 #endif
    721 {
    722 	Bigint *b;
    723 	U d;
    724 #ifndef Sudden_Underflow
    725 	int i;
    726 #endif
    727 	int de, k;
    728 	ULong *x, y, z;
    729 #ifdef VAX
    730 	ULong d0, d1;
    731 #else
    732 #define d0 word0(&d)
    733 #define d1 word1(&d)
    734 #endif
    735 	d.d = dd;
    736 #ifdef VAX
    737 	d0 = word0(&d) >> 16 | word0(&d) << 16;
    738 	d1 = word1(&d) >> 16 | word1(&d) << 16;
    739 #endif
    740 
    741 #ifdef Pack_32
    742 	b = Balloc(1);
    743 #else
    744 	b = Balloc(2);
    745 #endif
    746 	if (b == NULL)
    747 		return (NULL);
    748 	x = b->x;
    749 
    750 	z = d0 & Frac_mask;
    751 	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */
    752 #ifdef Sudden_Underflow
    753 	de = (int)(d0 >> Exp_shift);
    754 #ifndef IBM
    755 	z |= Exp_msk11;
    756 #endif
    757 #else
    758 	if ( (de = (int)(d0 >> Exp_shift)) !=0)
    759 		z |= Exp_msk1;
    760 #endif
    761 #ifdef Pack_32
    762 	if ( (y = d1) !=0) {
    763 		if ( (k = lo0bits(&y)) !=0) {
    764 			x[0] = y | z << (32 - k);
    765 			z >>= k;
    766 			}
    767 		else
    768 			x[0] = y;
    769 #ifndef Sudden_Underflow
    770 		i =
    771 #endif
    772 		     b->wds = (x[1] = z) !=0 ? 2 : 1;
    773 		}
    774 	else {
    775 		k = lo0bits(&z);
    776 		x[0] = z;
    777 #ifndef Sudden_Underflow
    778 		i =
    779 #endif
    780 		    b->wds = 1;
    781 		k += 32;
    782 		}
    783 #else
    784 	if ( (y = d1) !=0) {
    785 		if ( (k = lo0bits(&y)) !=0)
    786 			if (k >= 16) {
    787 				x[0] = y | z << 32 - k & 0xffff;
    788 				x[1] = z >> k - 16 & 0xffff;
    789 				x[2] = z >> k;
    790 				i = 2;
    791 				}
    792 			else {
    793 				x[0] = y & 0xffff;
    794 				x[1] = y >> 16 | z << 16 - k & 0xffff;
    795 				x[2] = z >> k & 0xffff;
    796 				x[3] = z >> k+16;
    797 				i = 3;
    798 				}
    799 		else {
    800 			x[0] = y & 0xffff;
    801 			x[1] = y >> 16;
    802 			x[2] = z & 0xffff;
    803 			x[3] = z >> 16;
    804 			i = 3;
    805 			}
    806 		}
    807 	else {
    808 #ifdef DEBUG
    809 		if (!z)
    810 			Bug("Zero passed to d2b");
    811 #endif
    812 		k = lo0bits(&z);
    813 		if (k >= 16) {
    814 			x[0] = z;
    815 			i = 0;
    816 			}
    817 		else {
    818 			x[0] = z & 0xffff;
    819 			x[1] = z >> 16;
    820 			i = 1;
    821 			}
    822 		k += 32;
    823 		}
    824 	while(!x[i])
    825 		--i;
    826 	b->wds = i + 1;
    827 #endif
    828 #ifndef Sudden_Underflow
    829 	if (de) {
    830 #endif
    831 #ifdef IBM
    832 		*e = (de - Bias - (P-1) << 2) + k;
    833 		*bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
    834 #else
    835 		*e = de - Bias - (P-1) + k;
    836 		*bits = P - k;
    837 #endif
    838 #ifndef Sudden_Underflow
    839 		}
    840 	else {
    841 		*e = de - Bias - (P-1) + 1 + k;
    842 #ifdef Pack_32
    843 		*bits = 32*i - hi0bits(x[i-1]);
    844 #else
    845 		*bits = (i+2)*16 - hi0bits(x[i]);
    846 #endif
    847 		}
    848 #endif
    849 	return b;
    850 	}
    851 #undef d0
    852 #undef d1
    853 
    854  CONST double
    855 #ifdef IEEE_Arith
    856 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
    857 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
    858 		};
    859 #else
    860 #ifdef IBM
    861 bigtens[] = { 1e16, 1e32, 1e64 };
    862 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
    863 #else
    864 bigtens[] = { 1e16, 1e32 };
    865 CONST double tinytens[] = { 1e-16, 1e-32 };
    866 #endif
    867 #endif
    868 
    869  CONST double
    870 tens[] = {
    871 		1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
    872 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
    873 		1e20, 1e21, 1e22
    874 #ifdef VAX
    875 		, 1e23, 1e24
    876 #endif
    877 		};
    878 
    879  char *
    880 #ifdef KR_headers
    881 strcp_D2A(a, b) char *a; char *b;
    882 #else
    883 strcp_D2A(char *a, CONST char *b)
    884 #endif
    885 {
    886 	while((*a = *b++))
    887 		a++;
    888 	return a;
    889 	}
    890 
    891 #ifdef NO_STRING_H
    892 
    893  Char *
    894 #ifdef KR_headers
    895 memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
    896 #else
    897 memcpy_D2A(void *a1, void *b1, size_t len)
    898 #endif
    899 {
    900 	char *a = (char*)a1, *ae = a + len;
    901 	char *b = (char*)b1, *a0 = a;
    902 	while(a < ae)
    903 		*a++ = *b++;
    904 	return a0;
    905 	}
    906 
    907 #endif /* NO_STRING_H */
    908