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