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