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