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 *
_Balloc(int k)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
_Bfree(Bigint * v)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 *
_multadd(Bigint * b,int m,int a)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
_hi0bits(register unsigned long x)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
lo0bits(unsigned long * y)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 *
_i2b(int i)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 *
_mult(Bigint * a,Bigint * b)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 *
_pow5mult(Bigint * b,int k)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 *
_lshift(Bigint * b,int k)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
_cmp(Bigint * a,Bigint * b)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 *
_diff(Bigint * a,Bigint * b)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 *
_d2b(double darg,int * e,int * bits)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
594