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