|
|
1.1 root 1: /* Common routines for _dtoa and strtod */
2:
3: #include "fconv.h"
4:
5: #ifdef DEBUG
6: #include <stdio.h>
7: #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
8: #endif
9:
10: double
11: _tens[] = {
12: 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
13: 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
14: 1e20, 1e21, 1e22
15: #ifdef VAX
16: , 1e23, 1e24
17: #endif
18: };
19:
20:
21: #ifdef IEEE_Arith
22: double _bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
23: double _tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
24: #else
25: #ifdef IBM
26: double _bigtens[] = { 1e16, 1e32, 1e64 };
27: double _tinytens[] = { 1e-16, 1e-32, 1e-64 };
28: #else
29: double _bigtens[] = { 1e16, 1e32 };
30: double _tinytens[] = { 1e-16, 1e-32 };
31: #endif
32: #endif
33:
34: static Bigint *freelist[Kmax+1];
35:
36: Bigint *
37: _Balloc(int k)
38: {
39: int x;
40: Bigint *rv;
41:
42: if (rv = freelist[k]) {
43: freelist[k] = rv->next;
44: }
45: else {
46: x = 1 << k;
47: rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long));
48: rv->k = k;
49: rv->maxwds = x;
50: }
51: rv->sign = rv->wds = 0;
52: return rv;
53: }
54:
55: void
56: _Bfree(Bigint *v)
57: {
58: if (v) {
59: v->next = freelist[v->k];
60: freelist[v->k] = v;
61: }
62: }
63:
64:
65: Bigint *
66: _multadd(Bigint *b, int m, int a) /* multiply by m and add a */
67: {
68: int i, wds;
69: unsigned long *x, y;
70: #ifdef Pack_32
71: unsigned long xi, z;
72: #endif
73: Bigint *b1;
74:
75: wds = b->wds;
76: x = b->x;
77: i = 0;
78: do {
79: #ifdef Pack_32
80: xi = *x;
81: y = (xi & 0xffff) * m + a;
82: z = (xi >> 16) * m + (y >> 16);
83: a = (int)(z >> 16);
84: *x++ = (z << 16) + (y & 0xffff);
85: #else
86: y = *x * m + a;
87: a = (int)(y >> 16);
88: *x++ = y & 0xffff;
89: #endif
90: }
91: while(++i < wds);
92: if (a) {
93: if (wds >= b->maxwds) {
94: b1 = Balloc(b->k+1);
95: Bcopy(b1, b);
96: Bfree(b);
97: b = b1;
98: }
99: b->x[wds++] = a;
100: b->wds = wds;
101: }
102: return b;
103: }
104:
105: int
106: _hi0bits(register unsigned long x)
107: {
108: register int k = 0;
109:
110: if (!(x & 0xffff0000)) {
111: k = 16;
112: x <<= 16;
113: }
114: if (!(x & 0xff000000)) {
115: k += 8;
116: x <<= 8;
117: }
118: if (!(x & 0xf0000000)) {
119: k += 4;
120: x <<= 4;
121: }
122: if (!(x & 0xc0000000)) {
123: k += 2;
124: x <<= 2;
125: }
126: if (!(x & 0x80000000)) {
127: k++;
128: if (!(x & 0x40000000))
129: return 32;
130: }
131: return k;
132: }
133:
134: static int
135: lo0bits(unsigned long *y)
136: {
137: register int k;
138: register unsigned long x = *y;
139:
140: if (x & 7) {
141: if (x & 1)
142: return 0;
143: if (x & 2) {
144: *y = x >> 1;
145: return 1;
146: }
147: *y = x >> 2;
148: return 2;
149: }
150: k = 0;
151: if (!(x & 0xffff)) {
152: k = 16;
153: x >>= 16;
154: }
155: if (!(x & 0xff)) {
156: k += 8;
157: x >>= 8;
158: }
159: if (!(x & 0xf)) {
160: k += 4;
161: x >>= 4;
162: }
163: if (!(x & 0x3)) {
164: k += 2;
165: x >>= 2;
166: }
167: if (!(x & 1)) {
168: k++;
169: x >>= 1;
170: if (!x & 1)
171: return 32;
172: }
173: *y = x;
174: return k;
175: }
176:
177: Bigint *
178: _i2b(int i)
179: {
180: Bigint *b;
181:
182: b = Balloc(1);
183: b->x[0] = i;
184: b->wds = 1;
185: return b;
186: }
187:
188: Bigint *
189: _mult(Bigint *a, Bigint *b)
190: {
191: Bigint *c;
192: int k, wa, wb, wc;
193: unsigned long carry, y, z;
194: unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
195: #ifdef Pack_32
196: unsigned long z2;
197: #endif
198:
199: if (a->wds < b->wds) {
200: c = a;
201: a = b;
202: b = c;
203: }
204: k = a->k;
205: wa = a->wds;
206: wb = b->wds;
207: wc = wa + wb;
208: if (wc > a->maxwds)
209: k++;
210: c = Balloc(k);
211: for(x = c->x, xa = x + wc; x < xa; x++)
212: *x = 0;
213: xa = a->x;
214: xae = xa + wa;
215: xb = b->x;
216: xbe = xb + wb;
217: xc0 = c->x;
218: #ifdef Pack_32
219: for(; xb < xbe; xb++, xc0++) {
220: if (y = *xb & 0xffff) {
221: x = xa;
222: xc = xc0;
223: carry = 0;
224: do {
225: z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
226: carry = z >> 16;
227: z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
228: carry = z2 >> 16;
229: Storeinc(xc, z2, z);
230: }
231: while(x < xae);
232: *xc = carry;
233: }
234: if (y = *xb >> 16) {
235: x = xa;
236: xc = xc0;
237: carry = 0;
238: z2 = *xc;
239: do {
240: z = (*x & 0xffff) * y + (*xc >> 16) + carry;
241: carry = z >> 16;
242: Storeinc(xc, z, z2);
243: z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
244: carry = z2 >> 16;
245: }
246: while(x < xae);
247: *xc = z2;
248: }
249: }
250: #else
251: for(; xb < xbe; xc0++) {
252: if (y = *xb++) {
253: x = xa;
254: xc = xc0;
255: carry = 0;
256: do {
257: z = *x++ * y + *xc + carry;
258: carry = z >> 16;
259: *xc++ = z & 0xffff;
260: }
261: while(x < xae);
262: *xc = carry;
263: }
264: }
265: #endif
266: for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
267: c->wds = wc;
268: return c;
269: }
270:
271: static Bigint *p5s;
272:
273: Bigint *
274: _pow5mult(Bigint *b, int k)
275: {
276: Bigint *b1, *p5, *p51;
277: int i;
278: static int p05[3] = { 5, 25, 125 };
279:
280: if (i = k & 3)
281: b = multadd(b, p05[i-1], 0);
282:
283: if (!(k >>= 2))
284: return b;
285: if (!(p5 = p5s)) {
286: /* first time */
287: p5 = p5s = i2b(625);
288: p5->next = 0;
289: }
290: for(;;) {
291: if (k & 1) {
292: b1 = mult(b, p5);
293: Bfree(b);
294: b = b1;
295: }
296: if (!(k >>= 1))
297: break;
298: if (!(p51 = p5->next)) {
299: p51 = p5->next = mult(p5,p5);
300: p51->next = 0;
301: }
302: p5 = p51;
303: }
304: return b;
305: }
306:
307: Bigint *
308: _lshift(Bigint *b, int k)
309: {
310: int i, k1, n, n1;
311: Bigint *b1;
312: unsigned long *x, *x1, *xe, z;
313:
314: #ifdef Pack_32
315: n = k >> 5;
316: #else
317: n = k >> 4;
318: #endif
319: k1 = b->k;
320: n1 = n + b->wds + 1;
321: for(i = b->maxwds; n1 > i; i <<= 1)
322: k1++;
323: b1 = Balloc(k1);
324: x1 = b1->x;
325: for(i = 0; i < n; i++)
326: *x1++ = 0;
327: x = b->x;
328: xe = x + b->wds;
329: #ifdef Pack_32
330: if (k &= 0x1f) {
331: k1 = 32 - k;
332: z = 0;
333: do {
334: *x1++ = *x << k | z;
335: z = *x++ >> k1;
336: }
337: while(x < xe);
338: if (*x1 = z)
339: ++n1;
340: }
341: #else
342: if (k &= 0xf) {
343: k1 = 16 - k;
344: z = 0;
345: do {
346: *x1++ = *x << k & 0xffff | z;
347: z = *x++ >> k1;
348: }
349: while(x < xe);
350: if (*x1 = z)
351: ++n1;
352: }
353: #endif
354: else do
355: *x1++ = *x++;
356: while(x < xe);
357: b1->wds = n1 - 1;
358: Bfree(b);
359: return b1;
360: }
361:
362: int
363: _cmp(Bigint *a, Bigint *b)
364: {
365: unsigned long *xa, *xa0, *xb, *xb0;
366: int i, j;
367:
368: i = a->wds;
369: j = b->wds;
370: #ifdef DEBUG
371: if (i > 1 && !a->x[i-1])
372: Bug("cmp called with a->x[a->wds-1] == 0");
373: if (j > 1 && !b->x[j-1])
374: Bug("cmp called with b->x[b->wds-1] == 0");
375: #endif
376: if (i -= j)
377: return i;
378: xa0 = a->x;
379: xa = xa0 + j;
380: xb0 = b->x;
381: xb = xb0 + j;
382: for(;;) {
383: if (*--xa != *--xb)
384: return *xa < *xb ? -1 : 1;
385: if (xa <= xa0)
386: break;
387: }
388: return 0;
389: }
390:
391: Bigint *
392: _diff(Bigint *a, Bigint *b)
393: {
394: Bigint *c;
395: int i, wa, wb;
396: long borrow, y; /* We need signed shifts here. */
397: unsigned long *xa, *xae, *xb, *xbe, *xc;
398: #ifdef Pack_32
399: long z;
400: #endif
401:
402: i = cmp(a,b);
403: if (!i) {
404: c = Balloc(0);
405: c->wds = 1;
406: c->x[0] = 0;
407: return c;
408: }
409: if (i < 0) {
410: c = a;
411: a = b;
412: b = c;
413: i = 1;
414: }
415: else
416: i = 0;
417: c = Balloc(a->k);
418: c->sign = i;
419: wa = a->wds;
420: xa = a->x;
421: xae = xa + wa;
422: wb = b->wds;
423: xb = b->x;
424: xbe = xb + wb;
425: xc = c->x;
426: borrow = 0;
427: #ifdef Pack_32
428: do {
429: y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
430: borrow = y >> 16;
431: Sign_Extend(borrow, y);
432: z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
433: borrow = z >> 16;
434: Sign_Extend(borrow, z);
435: Storeinc(xc, z, y);
436: }
437: while(xb < xbe);
438: while(xa < xae) {
439: y = (*xa & 0xffff) + borrow;
440: borrow = y >> 16;
441: Sign_Extend(borrow, y);
442: z = (*xa++ >> 16) + borrow;
443: borrow = z >> 16;
444: Sign_Extend(borrow, z);
445: Storeinc(xc, z, y);
446: }
447: #else
448: do {
449: y = *xa++ - *xb++ + borrow;
450: borrow = y >> 16;
451: Sign_Extend(borrow, y);
452: *xc++ = y & 0xffff;
453: }
454: while(xb < xbe);
455: while(xa < xae) {
456: y = *xa++ + borrow;
457: borrow = y >> 16;
458: Sign_Extend(borrow, y);
459: *xc++ = y & 0xffff;
460: }
461: #endif
462: while(!*--xc)
463: wa--;
464: c->wds = wa;
465: return c;
466: }
467:
468: Bigint *
469: _d2b(double darg, int *e, int *bits)
470: {
471: Bigint *b;
472: int de, i, k;
473: unsigned long *x, y, z;
474: Dul d;
475: #ifdef VAX
476: unsigned long d0, d1;
477: d.d = darg;
478: d0 = word0(d) >> 16 | word0(d) << 16;
479: d1 = word1(d) >> 16 | word1(d) << 16;
480: #else
481: d.d = darg;
482: #define d0 word0(d)
483: #define d1 word1(d)
484: #endif
485:
486: #ifdef Pack_32
487: b = Balloc(1);
488: #else
489: b = Balloc(2);
490: #endif
491: x = b->x;
492:
493: z = d0 & Frac_mask;
494: d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
495: #ifdef Sudden_Underflow
496: de = (int)(d0 >> Exp_shift);
497: #ifndef IBM
498: z |= Exp_msk11;
499: #endif
500: #else
501: if (de = (int)(d0 >> Exp_shift))
502: z |= Exp_msk1;
503: #endif
504: #ifdef Pack_32
505: if (y = d1) {
506: if (k = lo0bits(&y)) {
507: x[0] = y | z << 32 - k;
508: z >>= k;
509: }
510: else
511: x[0] = y;
512: i = b->wds = (x[1] = z) ? 2 : 1;
513: }
514: else {
515: #ifdef DEBUG
516: if (!z)
517: Bug("Zero passed to d2b");
518: #endif
519: k = lo0bits(&z);
520: x[0] = z;
521: i = b->wds = 1;
522: k += 32;
523: }
524: #else
525: if (y = d1) {
526: if (k = lo0bits(&y))
527: if (k >= 16) {
528: x[0] = y | z << 32 - k & 0xffff;
529: x[1] = z >> k - 16 & 0xffff;
530: x[2] = z >> k;
531: i = 2;
532: }
533: else {
534: x[0] = y & 0xffff;
535: x[1] = y >> 16 | z << 16 - k & 0xffff;
536: x[2] = z >> k & 0xffff;
537: x[3] = z >> k+16;
538: i = 3;
539: }
540: else {
541: x[0] = y & 0xffff;
542: x[1] = y >> 16;
543: x[2] = z & 0xffff;
544: x[3] = z >> 16;
545: i = 3;
546: }
547: }
548: else {
549: #ifdef DEBUG
550: if (!z)
551: Bug("Zero passed to d2b");
552: #endif
553: k = lo0bits(&z);
554: if (k >= 16) {
555: x[0] = z;
556: i = 0;
557: }
558: else {
559: x[0] = z & 0xffff;
560: x[1] = z >> 16;
561: i = 1;
562: }
563: k += 32;
564: }
565: while(!x[i])
566: --i;
567: b->wds = i + 1;
568: #endif
569: #ifndef Sudden_Underflow
570: if (de) {
571: #endif
572: #ifdef IBM
573: *e = (de - Bias - (P-1) << 2) + k;
574: *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
575: #else
576: *e = de - Bias - (P-1) + k;
577: *bits = P - k;
578: #endif
579: #ifndef Sudden_Underflow
580: }
581: else {
582: *e = de - Bias - (P-1) + 1 + k;
583: #ifdef Pack_32
584: *bits = 32*i - hi0bits(x[i-1]);
585: #else
586: *bits = (i+2)*16 - hi0bits(x[i]);
587: #endif
588: }
589: #endif
590: return b;
591: }
592: #undef d0
593: #undef d1
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.