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