Home | History | Annotate | Line # | Download | only in gdtoa
misc.c revision 1.7
      1 /* $NetBSD: misc.c,v 1.7 2011/03/21 04:52:09 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 				return NULL;
    420 			p5->next = 0;
    421 			}
    422 		FREE_DTOA_LOCK(1);
    423 #else
    424 		p5 = p5s = i2b(625);
    425 		if (p5 == NULL)
    426 			return NULL;
    427 		p5->next = 0;
    428 #endif
    429 		}
    430 	for(;;) {
    431 		if (k & 1) {
    432 			b1 = mult(b, p5);
    433 			if (b1 == NULL)
    434 				return NULL;
    435 			b = b1;
    436 			}
    437 		if (!(k = (unsigned int)k >> 1))
    438 			break;
    439 		if ((p51 = p5->next) == 0) {
    440 #ifdef MULTIPLE_THREADS
    441 			ACQUIRE_DTOA_LOCK(1);
    442 			if (!(p51 = p5->next)) {
    443 				p51 = p5->next = mult(p5,p5);
    444 				if (p51 == NULL)
    445 					return NULL;
    446 				p51->next = 0;
    447 				}
    448 			FREE_DTOA_LOCK(1);
    449 #else
    450 			p51 = p5->next = mult(p5,p5);
    451 			if (p51 == NULL)
    452 				return NULL;
    453 			p51->next = 0;
    454 #endif
    455 			}
    456 		p5 = p51;
    457 		}
    458 	return b;
    459 	}
    460 
    461  Bigint *
    462 lshift
    463 #ifdef KR_headers
    464 	(b, k) Bigint *b; int k;
    465 #else
    466 	(Bigint *b, int k)
    467 #endif
    468 {
    469 	int i, k1, n, n1;
    470 	Bigint *b1;
    471 	ULong *x, *x1, *xe, z;
    472 
    473 	n = (unsigned int)k >> kshift;
    474 	k1 = b->k;
    475 	n1 = n + b->wds + 1;
    476 	for(i = b->maxwds; n1 > i; i <<= 1)
    477 		k1++;
    478 	b1 = Balloc(k1);
    479 	if (b1 == NULL)
    480 		return NULL;
    481 	x1 = b1->x;
    482 	for(i = 0; i < n; i++)
    483 		*x1++ = 0;
    484 	x = b->x;
    485 	xe = x + b->wds;
    486 	if (k &= kmask) {
    487 #ifdef Pack_32
    488 		k1 = 32 - k;
    489 		z = 0;
    490 		do {
    491 			*x1++ = *x << k | z;
    492 			z = *x++ >> k1;
    493 			}
    494 			while(x < xe);
    495 		if ((*x1 = z) !=0)
    496 			++n1;
    497 #else
    498 		k1 = 16 - k;
    499 		z = 0;
    500 		do {
    501 			*x1++ = *x << k  & 0xffff | z;
    502 			z = *x++ >> k1;
    503 			}
    504 			while(x < xe);
    505 		if (*x1 = z)
    506 			++n1;
    507 #endif
    508 		}
    509 	else do
    510 		*x1++ = *x++;
    511 		while(x < xe);
    512 	b1->wds = n1 - 1;
    513 	Bfree(b);
    514 	return b1;
    515 	}
    516 
    517  int
    518 cmp
    519 #ifdef KR_headers
    520 	(a, b) Bigint *a, *b;
    521 #else
    522 	(Bigint *a, Bigint *b)
    523 #endif
    524 {
    525 	ULong *xa, *xa0, *xb, *xb0;
    526 	int i, j;
    527 
    528 	i = a->wds;
    529 	j = b->wds;
    530 #ifdef DEBUG
    531 	if (i > 1 && !a->x[i-1])
    532 		Bug("cmp called with a->x[a->wds-1] == 0");
    533 	if (j > 1 && !b->x[j-1])
    534 		Bug("cmp called with b->x[b->wds-1] == 0");
    535 #endif
    536 	if (i -= j)
    537 		return i;
    538 	xa0 = a->x;
    539 	xa = xa0 + j;
    540 	xb0 = b->x;
    541 	xb = xb0 + j;
    542 	for(;;) {
    543 		if (*--xa != *--xb)
    544 			return *xa < *xb ? -1 : 1;
    545 		if (xa <= xa0)
    546 			break;
    547 		}
    548 	return 0;
    549 	}
    550 
    551  Bigint *
    552 diff
    553 #ifdef KR_headers
    554 	(a, b) Bigint *a, *b;
    555 #else
    556 	(Bigint *a, Bigint *b)
    557 #endif
    558 {
    559 	Bigint *c;
    560 	int i, wa, wb;
    561 	ULong *xa, *xae, *xb, *xbe, *xc;
    562 #ifdef ULLong
    563 	ULLong borrow, y;
    564 #else
    565 	ULong borrow, y;
    566 #ifdef Pack_32
    567 	ULong z;
    568 #endif
    569 #endif
    570 
    571 	i = cmp(a,b);
    572 	if (!i) {
    573 		c = Balloc(0);
    574 		if (c == NULL)
    575 			return NULL;
    576 		c->wds = 1;
    577 		c->x[0] = 0;
    578 		return c;
    579 		}
    580 	if (i < 0) {
    581 		c = a;
    582 		a = b;
    583 		b = c;
    584 		i = 1;
    585 		}
    586 	else
    587 		i = 0;
    588 	c = Balloc(a->k);
    589 	if (c == NULL)
    590 		return NULL;
    591 	c->sign = i;
    592 	wa = a->wds;
    593 	xa = a->x;
    594 	xae = xa + wa;
    595 	wb = b->wds;
    596 	xb = b->x;
    597 	xbe = xb + wb;
    598 	xc = c->x;
    599 	borrow = 0;
    600 #ifdef ULLong
    601 	do {
    602 		y = (ULLong)*xa++ - *xb++ - borrow;
    603 		borrow = y >> 32 & 1UL;
    604 		/* LINTED conversion */
    605 		*xc++ = y & 0xffffffffUL;
    606 		}
    607 		while(xb < xbe);
    608 	while(xa < xae) {
    609 		y = *xa++ - borrow;
    610 		borrow = y >> 32 & 1UL;
    611 		/* LINTED conversion */
    612 		*xc++ = y & 0xffffffffUL;
    613 		}
    614 #else
    615 #ifdef Pack_32
    616 	do {
    617 		y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
    618 		borrow = (y & 0x10000) >> 16;
    619 		z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
    620 		borrow = (z & 0x10000) >> 16;
    621 		Storeinc(xc, z, y);
    622 		}
    623 		while(xb < xbe);
    624 	while(xa < xae) {
    625 		y = (*xa & 0xffff) - borrow;
    626 		borrow = (y & 0x10000) >> 16;
    627 		z = (*xa++ >> 16) - borrow;
    628 		borrow = (z & 0x10000) >> 16;
    629 		Storeinc(xc, z, y);
    630 		}
    631 #else
    632 	do {
    633 		y = *xa++ - *xb++ - borrow;
    634 		borrow = (y & 0x10000) >> 16;
    635 		*xc++ = y & 0xffff;
    636 		}
    637 		while(xb < xbe);
    638 	while(xa < xae) {
    639 		y = *xa++ - borrow;
    640 		borrow = (y & 0x10000) >> 16;
    641 		*xc++ = y & 0xffff;
    642 		}
    643 #endif
    644 #endif
    645 	while(!*--xc)
    646 		wa--;
    647 	c->wds = wa;
    648 	return c;
    649 	}
    650 
    651  double
    652 b2d
    653 #ifdef KR_headers
    654 	(a, e) Bigint *a; int *e;
    655 #else
    656 	(Bigint *a, int *e)
    657 #endif
    658 {
    659 	ULong *xa, *xa0, w, y, z;
    660 	int k;
    661 	U d;
    662 #ifdef VAX
    663 	ULong d0, d1;
    664 #else
    665 #define d0 word0(&d)
    666 #define d1 word1(&d)
    667 #endif
    668 
    669 	xa0 = a->x;
    670 	xa = xa0 + a->wds;
    671 	y = *--xa;
    672 #ifdef DEBUG
    673 	if (!y) Bug("zero y in b2d");
    674 #endif
    675 	k = hi0bits(y);
    676 	*e = 32 - k;
    677 #ifdef Pack_32
    678 	if (k < Ebits) {
    679 		d0 = Exp_1 | y >> (Ebits - k);
    680 		w = xa > xa0 ? *--xa : 0;
    681 		d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
    682 		goto ret_d;
    683 		}
    684 	z = xa > xa0 ? *--xa : 0;
    685 	if (k -= Ebits) {
    686 		d0 = Exp_1 | y << k | z >> (32 - k);
    687 		y = xa > xa0 ? *--xa : 0;
    688 		d1 = z << k | y >> (32 - k);
    689 		}
    690 	else {
    691 		d0 = Exp_1 | y;
    692 		d1 = z;
    693 		}
    694 #else
    695 	if (k < Ebits + 16) {
    696 		z = xa > xa0 ? *--xa : 0;
    697 		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
    698 		w = xa > xa0 ? *--xa : 0;
    699 		y = xa > xa0 ? *--xa : 0;
    700 		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
    701 		goto ret_d;
    702 		}
    703 	z = xa > xa0 ? *--xa : 0;
    704 	w = xa > xa0 ? *--xa : 0;
    705 	k -= Ebits + 16;
    706 	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
    707 	y = xa > xa0 ? *--xa : 0;
    708 	d1 = w << k + 16 | y << k;
    709 #endif
    710  ret_d:
    711 #ifdef VAX
    712 	word0(&d) = d0 >> 16 | d0 << 16;
    713 	word1(&d) = d1 >> 16 | d1 << 16;
    714 #endif
    715 	return dval(&d);
    716 	}
    717 #undef d0
    718 #undef d1
    719 
    720  Bigint *
    721 d2b
    722 #ifdef KR_headers
    723 	(dd, e, bits) double dd; int *e, *bits;
    724 #else
    725 	(double dd, int *e, int *bits)
    726 #endif
    727 {
    728 	Bigint *b;
    729 	U d;
    730 #ifndef Sudden_Underflow
    731 	int i;
    732 #endif
    733 	int de, k;
    734 	ULong *x, y, z;
    735 #ifdef VAX
    736 	ULong d0, d1;
    737 #else
    738 #define d0 word0(&d)
    739 #define d1 word1(&d)
    740 #endif
    741 	d.d = dd;
    742 #ifdef VAX
    743 	d0 = word0(&d) >> 16 | word0(&d) << 16;
    744 	d1 = word1(&d) >> 16 | word1(&d) << 16;
    745 #endif
    746 
    747 #ifdef Pack_32
    748 	b = Balloc(1);
    749 #else
    750 	b = Balloc(2);
    751 #endif
    752 	if (b == NULL)
    753 		return NULL;
    754 	x = b->x;
    755 
    756 	z = d0 & Frac_mask;
    757 	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */
    758 #ifdef Sudden_Underflow
    759 	de = (int)(d0 >> Exp_shift);
    760 #ifndef IBM
    761 	z |= Exp_msk11;
    762 #endif
    763 #else
    764 	if ( (de = (int)(d0 >> Exp_shift)) !=0)
    765 		z |= Exp_msk1;
    766 #endif
    767 #ifdef Pack_32
    768 	if ( (y = d1) !=0) {
    769 		if ( (k = lo0bits(&y)) !=0) {
    770 			x[0] = y | z << (32 - k);
    771 			z >>= k;
    772 			}
    773 		else
    774 			x[0] = y;
    775 #ifndef Sudden_Underflow
    776 		i =
    777 #endif
    778 		     b->wds = (x[1] = z) !=0 ? 2 : 1;
    779 		}
    780 	else {
    781 		k = lo0bits(&z);
    782 		x[0] = z;
    783 #ifndef Sudden_Underflow
    784 		i =
    785 #endif
    786 		    b->wds = 1;
    787 		k += 32;
    788 		}
    789 #else
    790 	if ( (y = d1) !=0) {
    791 		if ( (k = lo0bits(&y)) !=0)
    792 			if (k >= 16) {
    793 				x[0] = y | z << 32 - k & 0xffff;
    794 				x[1] = z >> k - 16 & 0xffff;
    795 				x[2] = z >> k;
    796 				i = 2;
    797 				}
    798 			else {
    799 				x[0] = y & 0xffff;
    800 				x[1] = y >> 16 | z << 16 - k & 0xffff;
    801 				x[2] = z >> k & 0xffff;
    802 				x[3] = z >> k+16;
    803 				i = 3;
    804 				}
    805 		else {
    806 			x[0] = y & 0xffff;
    807 			x[1] = y >> 16;
    808 			x[2] = z & 0xffff;
    809 			x[3] = z >> 16;
    810 			i = 3;
    811 			}
    812 		}
    813 	else {
    814 #ifdef DEBUG
    815 		if (!z)
    816 			Bug("Zero passed to d2b");
    817 #endif
    818 		k = lo0bits(&z);
    819 		if (k >= 16) {
    820 			x[0] = z;
    821 			i = 0;
    822 			}
    823 		else {
    824 			x[0] = z & 0xffff;
    825 			x[1] = z >> 16;
    826 			i = 1;
    827 			}
    828 		k += 32;
    829 		}
    830 	while(!x[i])
    831 		--i;
    832 	b->wds = i + 1;
    833 #endif
    834 #ifndef Sudden_Underflow
    835 	if (de) {
    836 #endif
    837 #ifdef IBM
    838 		*e = (de - Bias - (P-1) << 2) + k;
    839 		*bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
    840 #else
    841 		*e = de - Bias - (P-1) + k;
    842 		*bits = P - k;
    843 #endif
    844 #ifndef Sudden_Underflow
    845 		}
    846 	else {
    847 		*e = de - Bias - (P-1) + 1 + k;
    848 #ifdef Pack_32
    849 		*bits = 32*i - hi0bits(x[i-1]);
    850 #else
    851 		*bits = (i+2)*16 - hi0bits(x[i]);
    852 #endif
    853 		}
    854 #endif
    855 	return b;
    856 	}
    857 #undef d0
    858 #undef d1
    859 
    860  CONST double
    861 #ifdef IEEE_Arith
    862 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
    863 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
    864 		};
    865 #else
    866 #ifdef IBM
    867 bigtens[] = { 1e16, 1e32, 1e64 };
    868 CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
    869 #else
    870 bigtens[] = { 1e16, 1e32 };
    871 CONST double tinytens[] = { 1e-16, 1e-32 };
    872 #endif
    873 #endif
    874 
    875  CONST double
    876 tens[] = {
    877 		1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
    878 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
    879 		1e20, 1e21, 1e22
    880 #ifdef VAX
    881 		, 1e23, 1e24
    882 #endif
    883 		};
    884 
    885  char *
    886 #ifdef KR_headers
    887 strcp_D2A(a, b) char *a; char *b;
    888 #else
    889 strcp_D2A(char *a, CONST char *b)
    890 #endif
    891 {
    892 	while((*a = *b++))
    893 		a++;
    894 	return a;
    895 	}
    896 
    897 #ifdef NO_STRING_H
    898 
    899  Char *
    900 #ifdef KR_headers
    901 memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
    902 #else
    903 memcpy_D2A(void *a1, void *b1, size_t len)
    904 #endif
    905 {
    906 	char *a = (char*)a1, *ae = a + len;
    907 	char *b = (char*)b1, *a0 = a;
    908 	while(a < ae)
    909 		*a++ = *b++;
    910 	return a0;
    911 	}
    912 
    913 #endif /* NO_STRING_H */
    914