Home | History | Annotate | Line # | Download | only in gdtoa
misc.c revision 1.2 
      1 /* $NetBSD: misc.c,v 1.2 2006/01/25 15:27:42 kleink 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 (pmem_next - private_mem + len <= 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 		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 		ACQUIRE_DTOA_LOCK(0);
     95 		v->next = freelist[v->k];
     96 		freelist[v->k] = v;
     97 		FREE_DTOA_LOCK(0);
     98 		}
     99 	}
    100 
    101  int
    102 lo0bits
    103 #ifdef KR_headers
    104 	(y) ULong *y;
    105 #else
    106 	(ULong *y)
    107 #endif
    108 {
    109 	int k;
    110 	ULong x = *y;
    111 
    112 	if (x & 7) {
    113 		if (x & 1)
    114 			return 0;
    115 		if (x & 2) {
    116 			*y = x >> 1;
    117 			return 1;
    118 			}
    119 		*y = x >> 2;
    120 		return 2;
    121 		}
    122 	k = 0;
    123 	if (!(x & 0xffff)) {
    124 		k = 16;
    125 		x >>= 16;
    126 		}
    127 	if (!(x & 0xff)) {
    128 		k += 8;
    129 		x >>= 8;
    130 		}
    131 	if (!(x & 0xf)) {
    132 		k += 4;
    133 		x >>= 4;
    134 		}
    135 	if (!(x & 0x3)) {
    136 		k += 2;
    137 		x >>= 2;
    138 		}
    139 	if (!(x & 1)) {
    140 		k++;
    141 		x >>= 1;
    142 		if (!x)
    143 			return 32;
    144 		}
    145 	*y = x;
    146 	return k;
    147 	}
    148 
    149  Bigint *
    150 multadd
    151 #ifdef KR_headers
    152 	(b, m, a) Bigint *b; int m, a;
    153 #else
    154 	(Bigint *b, int m, int a)	/* multiply by m and add a */
    155 #endif
    156 {
    157 	int i, wds;
    158 #ifdef ULLong
    159 	ULong *x;
    160 	ULLong carry, y;
    161 #else
    162 	ULong carry, *x, y;
    163 #ifdef Pack_32
    164 	ULong xi, z;
    165 #endif
    166 #endif
    167 	Bigint *b1;
    168 
    169 	wds = b->wds;
    170 	x = b->x;
    171 	i = 0;
    172 	carry = a;
    173 	do {
    174 #ifdef ULLong
    175 		y = *x * (ULLong)m + carry;
    176 		carry = y >> 32;
    177 		/* LINTED conversion */
    178 		*x++ = y & 0xffffffffUL;
    179 #else
    180 #ifdef Pack_32
    181 		xi = *x;
    182 		y = (xi & 0xffff) * m + carry;
    183 		z = (xi >> 16) * m + (y >> 16);
    184 		carry = z >> 16;
    185 		*x++ = (z << 16) + (y & 0xffff);
    186 #else
    187 		y = *x * m + carry;
    188 		carry = y >> 16;
    189 		*x++ = y & 0xffff;
    190 #endif
    191 #endif
    192 		}
    193 		while(++i < wds);
    194 	if (carry) {
    195 		if (wds >= b->maxwds) {
    196 			b1 = Balloc(b->k+1);
    197 			Bcopy(b1, b);
    198 			Bfree(b);
    199 			b = b1;
    200 			}
    201 		/* LINTED conversion */
    202 		b->x[wds++] = carry;
    203 		b->wds = wds;
    204 		}
    205 	return b;
    206 	}
    207 
    208  int
    209 hi0bits_D2A
    210 #ifdef KR_headers
    211 	(x) ULong x;
    212 #else
    213 	(ULong x)
    214 #endif
    215 {
    216 	int k = 0;
    217 
    218 	if (!(x & 0xffff0000)) {
    219 		k = 16;
    220 		x <<= 16;
    221 		}
    222 	if (!(x & 0xff000000)) {
    223 		k += 8;
    224 		x <<= 8;
    225 		}
    226 	if (!(x & 0xf0000000)) {
    227 		k += 4;
    228 		x <<= 4;
    229 		}
    230 	if (!(x & 0xc0000000)) {
    231 		k += 2;
    232 		x <<= 2;
    233 		}
    234 	if (!(x & 0x80000000)) {
    235 		k++;
    236 		if (!(x & 0x40000000))
    237 			return 32;
    238 		}
    239 	return k;
    240 	}
    241 
    242  Bigint *
    243 i2b
    244 #ifdef KR_headers
    245 	(i) int i;
    246 #else
    247 	(int i)
    248 #endif
    249 {
    250 	Bigint *b;
    251 
    252 	b = Balloc(1);
    253 	b->x[0] = i;
    254 	b->wds = 1;
    255 	return b;
    256 	}
    257 
    258  Bigint *
    259 mult
    260 #ifdef KR_headers
    261 	(a, b) Bigint *a, *b;
    262 #else
    263 	(Bigint *a, Bigint *b)
    264 #endif
    265 {
    266 	Bigint *c;
    267 	int k, wa, wb, wc;
    268 	ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
    269 	ULong y;
    270 #ifdef ULLong
    271 	ULLong carry, z;
    272 #else
    273 	ULong carry, z;
    274 #ifdef Pack_32
    275 	ULong z2;
    276 #endif
    277 #endif
    278 
    279 	if (a->wds < b->wds) {
    280 		c = a;
    281 		a = b;
    282 		b = c;
    283 		}
    284 	k = a->k;
    285 	wa = a->wds;
    286 	wb = b->wds;
    287 	wc = wa + wb;
    288 	if (wc > a->maxwds)
    289 		k++;
    290 	c = Balloc(k);
    291 	for(x = c->x, xa = x + wc; x < xa; x++)
    292 		*x = 0;
    293 	xa = a->x;
    294 	xae = xa + wa;
    295 	xb = b->x;
    296 	xbe = xb + wb;
    297 	xc0 = c->x;
    298 #ifdef ULLong
    299 	for(; xb < xbe; xc0++) {
    300 		if ( (y = *xb++) !=0) {
    301 			x = xa;
    302 			xc = xc0;
    303 			carry = 0;
    304 			do {
    305 				z = *x++ * (ULLong)y + *xc + carry;
    306 				carry = z >> 32;
    307 				/* LINTED conversion */
    308 				*xc++ = z & 0xffffffffUL;
    309 				}
    310 				while(x < xae);
    311 			/* LINTED conversion */
    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 	CONST 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 = (unsigned int)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 = (unsigned int)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 = (unsigned int)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 		/* LINTED conversion */
    567 		*xc++ = y & 0xffffffffUL;
    568 		}
    569 		while(xb < xbe);
    570 	while(xa < xae) {
    571 		y = *xa++ - borrow;
    572 		borrow = y >> 32 & 1UL;
    573 		/* LINTED conversion */
    574 		*xc++ = y & 0xffffffffUL;
    575 		}
    576 #else
    577 #ifdef Pack_32
    578 	do {
    579 		y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
    580 		borrow = (y & 0x10000) >> 16;
    581 		z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
    582 		borrow = (z & 0x10000) >> 16;
    583 		Storeinc(xc, z, y);
    584 		}
    585 		while(xb < xbe);
    586 	while(xa < xae) {
    587 		y = (*xa & 0xffff) - borrow;
    588 		borrow = (y & 0x10000) >> 16;
    589 		z = (*xa++ >> 16) - borrow;
    590 		borrow = (z & 0x10000) >> 16;
    591 		Storeinc(xc, z, y);
    592 		}
    593 #else
    594 	do {
    595 		y = *xa++ - *xb++ - borrow;
    596 		borrow = (y & 0x10000) >> 16;
    597 		*xc++ = y & 0xffff;
    598 		}
    599 		while(xb < xbe);
    600 	while(xa < xae) {
    601 		y = *xa++ - borrow;
    602 		borrow = (y & 0x10000) >> 16;
    603 		*xc++ = y & 0xffff;
    604 		}
    605 #endif
    606 #endif
    607 	while(!*--xc)
    608 		wa--;
    609 	c->wds = wa;
    610 	return c;
    611 	}
    612 
    613  double
    614 b2d
    615 #ifdef KR_headers
    616 	(a, e) Bigint *a; int *e;
    617 #else
    618 	(Bigint *a, int *e)
    619 #endif
    620 {
    621 	ULong *xa, *xa0, w, y, z;
    622 	int k;
    623 	double d;
    624 #ifdef VAX
    625 	ULong d0, d1;
    626 #else
    627 #define d0 word0(d)
    628 #define d1 word1(d)
    629 #endif
    630 
    631 	xa0 = a->x;
    632 	xa = xa0 + a->wds;
    633 	y = *--xa;
    634 #ifdef DEBUG
    635 	if (!y) Bug("zero y in b2d");
    636 #endif
    637 	k = hi0bits(y);
    638 	*e = 32 - k;
    639 #ifdef Pack_32
    640 	if (k < Ebits) {
    641 		d0 = Exp_1 | y >> (Ebits - k);
    642 		w = xa > xa0 ? *--xa : 0;
    643 		d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
    644 		goto ret_d;
    645 		}
    646 	z = xa > xa0 ? *--xa : 0;
    647 	if (k -= Ebits) {
    648 		d0 = Exp_1 | y << k | z >> (32 - k);
    649 		y = xa > xa0 ? *--xa : 0;
    650 		d1 = z << k | y >> (32 - k);
    651 		}
    652 	else {
    653 		d0 = Exp_1 | y;
    654 		d1 = z;
    655 		}
    656 #else
    657 	if (k < Ebits + 16) {
    658 		z = xa > xa0 ? *--xa : 0;
    659 		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
    660 		w = xa > xa0 ? *--xa : 0;
    661 		y = xa > xa0 ? *--xa : 0;
    662 		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
    663 		goto ret_d;
    664 		}
    665 	z = xa > xa0 ? *--xa : 0;
    666 	w = xa > xa0 ? *--xa : 0;
    667 	k -= Ebits + 16;
    668 	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
    669 	y = xa > xa0 ? *--xa : 0;
    670 	d1 = w << k + 16 | y << k;
    671 #endif
    672  ret_d:
    673 #ifdef VAX
    674 	word0(d) = d0 >> 16 | d0 << 16;
    675 	word1(d) = d1 >> 16 | d1 << 16;
    676 #endif
    677 	return dval(d);
    678 	}
    679 #undef d0
    680 #undef d1
    681 
    682  Bigint *
    683 d2b
    684 #ifdef KR_headers
    685 	(d, e, bits) double d; int *e, *bits;
    686 #else
    687 	(double d, int *e, int *bits)
    688 #endif
    689 {
    690 	Bigint *b;
    691 #ifndef Sudden_Underflow
    692 	int i;
    693 #endif
    694 	int de, k;
    695 	ULong *x, y, z;
    696 #ifdef VAX
    697 	ULong d0, d1;
    698 	d0 = word0(d) >> 16 | word0(d) << 16;
    699 	d1 = word1(d) >> 16 | word1(d) << 16;
    700 #else
    701 #define d0 word0(d)
    702 #define d1 word1(d)
    703 #endif
    704 
    705 #ifdef Pack_32
    706 	b = Balloc(1);
    707 #else
    708 	b = Balloc(2);
    709 #endif
    710 	x = b->x;
    711 
    712 	z = d0 & Frac_mask;
    713 	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */
    714 #ifdef Sudden_Underflow
    715 	de = (int)(d0 >> Exp_shift);
    716 #ifndef IBM
    717 	z |= Exp_msk11;
    718 #endif
    719 #else
    720 	if ( (de = (int)(d0 >> Exp_shift)) !=0)
    721 		z |= Exp_msk1;
    722 #endif
    723 #ifdef Pack_32
    724 	if ( (y = d1) !=0) {
    725 		if ( (k = lo0bits(&y)) !=0) {
    726 			x[0] = y | z << (32 - k);
    727 			z >>= k;
    728 			}
    729 		else
    730 			x[0] = y;
    731 #ifndef Sudden_Underflow
    732 		i =
    733 #endif
    734 		     b->wds = (x[1] = z) !=0 ? 2 : 1;
    735 		}
    736 	else {
    737 #ifdef DEBUG
    738 		if (!z)
    739 			Bug("Zero passed to d2b");
    740 #endif
    741 		k = lo0bits(&z);
    742 		x[0] = z;
    743 #ifndef Sudden_Underflow
    744 		i =
    745 #endif
    746 		    b->wds = 1;
    747 		k += 32;
    748 		}
    749 #else
    750 	if ( (y = d1) !=0) {
    751 		if ( (k = lo0bits(&y)) !=0)
    752 			if (k >= 16) {
    753 				x[0] = y | z << 32 - k & 0xffff;
    754 				x[1] = z >> k - 16 & 0xffff;
    755 				x[2] = z >> k;
    756 				i = 2;
    757 				}
    758 			else {
    759 				x[0] = y & 0xffff;
    760 				x[1] = y >> 16 | z << 16 - k & 0xffff;
    761 				x[2] = z >> k & 0xffff;
    762 				x[3] = z >> k+16;
    763 				i = 3;
    764 				}
    765 		else {
    766 			x[0] = y & 0xffff;
    767 			x[1] = y >> 16;
    768 			x[2] = z & 0xffff;
    769 			x[3] = z >> 16;
    770 			i = 3;
    771 			}
    772 		}
    773 	else {
    774 #ifdef DEBUG
    775 		if (!z)
    776 			Bug("Zero passed to d2b");
    777 #endif
    778 		k = lo0bits(&z);
    779 		if (k >= 16) {
    780 			x[0] = z;
    781 			i = 0;
    782 			}
    783 		else {
    784 			x[0] = z & 0xffff;
    785 			x[1] = z >> 16;
    786 			i = 1;
    787 			}
    788 		k += 32;
    789 		}
    790 	while(!x[i])
    791 		--i;
    792 	b->wds = i + 1;
    793 #endif
    794 #ifndef Sudden_Underflow
    795 	if (de) {
    796 #endif
    797 #ifdef IBM
    798 		*e = (de - Bias - (P-1) << 2) + k;
    799 		*bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
    800 #else
    801 		*e = de - Bias - (P-1) + k;
    802 		*bits = P - k;
    803 #endif
    804 #ifndef Sudden_Underflow
    805 		}
    806 	else {
    807 		*e = de - Bias - (P-1) + 1 + k;
    808 #ifdef Pack_32
    809 		*bits = 32*i - hi0bits(x[i-1]);
    810 #else
    811 		*bits = (i+2)*16 - hi0bits(x[i]);
    812 #endif
    813 		}
    814 #endif
    815 	return b;
    816 	}
    817 #undef d0
    818 #undef d1
    819 
    820  CONST double
    821 #ifdef IEEE_Arith
    822 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
    823 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
    824 		};
    825 #else
    826 #ifdef IBM
    827 bigtens[] = { 1e16, 1e32, 1e64 };
    828 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
    829 #else
    830 bigtens[] = { 1e16, 1e32 };
    831 CONST double tinytens[] = { 1e-16, 1e-32 };
    832 #endif
    833 #endif
    834 
    835  CONST double
    836 tens[] = {
    837 		1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
    838 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
    839 		1e20, 1e21, 1e22
    840 #ifdef VAX
    841 		, 1e23, 1e24
    842 #endif
    843 		};
    844 
    845  char *
    846 #ifdef KR_headers
    847 strcp_D2A(a, b) char *a; char *b;
    848 #else
    849 strcp_D2A(char *a, CONST char *b)
    850 #endif
    851 {
    852 	while((*a = *b++))
    853 		a++;
    854 	return a;
    855 	}
    856 
    857 #ifdef NO_STRING_H
    858 
    859  Char *
    860 #ifdef KR_headers
    861 memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
    862 #else
    863 memcpy_D2A(void *a1, void *b1, size_t len)
    864 #endif
    865 {
    866 	char *a = (char*)a1, *ae = a + len;
    867 	char *b = (char*)b1, *a0 = a;
    868 	while(a < ae)
    869 		*a++ = *b++;
    870 	return a0;
    871 	}
    872 
    873 #endif /* NO_STRING_H */
    874