xref: /netbsd-src/lib/libc/gdtoa/misc.c (revision b2707ef3ebc49c3a5af380e63ffb1753ac6def01)
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