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