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