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