misc.c revision 196916 
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			free((void*)v);
96		else {
97			ACQUIRE_DTOA_LOCK(0);
98			v->next = freelist[v->k];
99			freelist[v->k] = v;
100			FREE_DTOA_LOCK(0);
101			}
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)
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_D2A
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#ifndef Sudden_Underflow
690	int i;
691#endif
692	int de, k;
693	ULong *x, y, z;
694#ifdef VAX
695	ULong d0, d1;
696	d0 = word0(d) >> 16 | word0(d) << 16;
697	d1 = word1(d) >> 16 | word1(d) << 16;
698#else
699#define d0 word0(d)
700#define d1 word1(d)
701#endif
702
703#ifdef Pack_32
704	b = Balloc(1);
705#else
706	b = Balloc(2);
707#endif
708	x = b->x;
709
710	z = d0 & Frac_mask;
711	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */
712#ifdef Sudden_Underflow
713	de = (int)(d0 >> Exp_shift);
714#ifndef IBM
715	z |= Exp_msk11;
716#endif
717#else
718	if ( (de = (int)(d0 >> Exp_shift)) !=0)
719		z |= Exp_msk1;
720#endif
721#ifdef Pack_32
722	if ( (y = d1) !=0) {
723		if ( (k = lo0bits(&y)) !=0) {
724			x[0] = y | z << 32 - k;
725			z >>= k;
726			}
727		else
728			x[0] = y;
729#ifndef Sudden_Underflow
730		i =
731#endif
732		     b->wds = (x[1] = z) !=0 ? 2 : 1;
733		}
734	else {
735#ifdef DEBUG
736		if (!z)
737			Bug("Zero passed to d2b");
738#endif
739		k = lo0bits(&z);
740		x[0] = z;
741#ifndef Sudden_Underflow
742		i =
743#endif
744		    b->wds = 1;
745		k += 32;
746		}
747#else
748	if ( (y = d1) !=0) {
749		if ( (k = lo0bits(&y)) !=0)
750			if (k >= 16) {
751				x[0] = y | z << 32 - k & 0xffff;
752				x[1] = z >> k - 16 & 0xffff;
753				x[2] = z >> k;
754				i = 2;
755				}
756			else {
757				x[0] = y & 0xffff;
758				x[1] = y >> 16 | z << 16 - k & 0xffff;
759				x[2] = z >> k & 0xffff;
760				x[3] = z >> k+16;
761				i = 3;
762				}
763		else {
764			x[0] = y & 0xffff;
765			x[1] = y >> 16;
766			x[2] = z & 0xffff;
767			x[3] = z >> 16;
768			i = 3;
769			}
770		}
771	else {
772#ifdef DEBUG
773		if (!z)
774			Bug("Zero passed to d2b");
775#endif
776		k = lo0bits(&z);
777		if (k >= 16) {
778			x[0] = z;
779			i = 0;
780			}
781		else {
782			x[0] = z & 0xffff;
783			x[1] = z >> 16;
784			i = 1;
785			}
786		k += 32;
787		}
788	while(!x[i])
789		--i;
790	b->wds = i + 1;
791#endif
792#ifndef Sudden_Underflow
793	if (de) {
794#endif
795#ifdef IBM
796		*e = (de - Bias - (P-1) << 2) + k;
797		*bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
798#else
799		*e = de - Bias - (P-1) + k;
800		*bits = P - k;
801#endif
802#ifndef Sudden_Underflow
803		}
804	else {
805		*e = de - Bias - (P-1) + 1 + k;
806#ifdef Pack_32
807		*bits = 32*i - hi0bits(x[i-1]);
808#else
809		*bits = (i+2)*16 - hi0bits(x[i]);
810#endif
811		}
812#endif
813	return b;
814	}
815#undef d0
816#undef d1
817
818 CONST double
819#ifdef IEEE_Arith
820bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
821CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
822		};
823#else
824#ifdef IBM
825bigtens[] = { 1e16, 1e32, 1e64 };
826CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
827#else
828bigtens[] = { 1e16, 1e32 };
829CONST double tinytens[] = { 1e-16, 1e-32 };
830#endif
831#endif
832
833 CONST double
834tens[] = {
835		1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
836		1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
837		1e20, 1e21, 1e22
838#ifdef VAX
839		, 1e23, 1e24
840#endif
841		};
842
843 char *
844#ifdef KR_headers
845strcp_D2A(a, b) char *a; char *b;
846#else
847strcp_D2A(char *a, CONST char *b)
848#endif
849{
850	while(*a = *b++)
851		a++;
852	return a;
853	}
854
855#ifdef NO_STRING_H
856
857 Char *
858#ifdef KR_headers
859memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
860#else
861memcpy_D2A(void *a1, void *b1, size_t len)
862#endif
863{
864	register char *a = (char*)a1, *ae = a + len;
865	register char *b = (char*)b1, *a0 = a;
866	while(a < ae)
867		*a++ = *b++;
868	return a0;
869	}
870
871#endif /* NO_STRING_H */
872