1 | /**************************************************************** |
---|
2 | * |
---|
3 | * The author of this software is David M. Gay. |
---|
4 | * |
---|
5 | * Copyright (c) 1991 by AT&T. |
---|
6 | * |
---|
7 | * Permission to use, copy, modify, and distribute this software for any |
---|
8 | * purpose without fee is hereby granted, provided that this entire notice |
---|
9 | * is included in all copies of any software which is or includes a copy |
---|
10 | * or modification of this software and in all copies of the supporting |
---|
11 | * documentation for such software. |
---|
12 | * |
---|
13 | * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED |
---|
14 | * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY |
---|
15 | * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY |
---|
16 | * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. |
---|
17 | * |
---|
18 | ***************************************************************/ |
---|
19 | |
---|
20 | /* Please send bug reports to |
---|
21 | David M. Gay |
---|
22 | AT&T Bell Laboratories, Room 2C-463 |
---|
23 | 600 Mountain Avenue |
---|
24 | Murray Hill, NJ 07974-2070 |
---|
25 | U.S.A. |
---|
26 | dmg@research.att.com or research!dmg |
---|
27 | */ |
---|
28 | |
---|
29 | /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. |
---|
30 | * |
---|
31 | * This strtod returns a nearest machine number to the input decimal |
---|
32 | * string (or sets errno to ERANGE). With IEEE arithmetic, ties are |
---|
33 | * broken by the IEEE round-even rule. Otherwise ties are broken by |
---|
34 | * biased rounding (add half and chop). |
---|
35 | * |
---|
36 | * Inspired loosely by William D. Clinger's paper "How to Read Floating |
---|
37 | * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. |
---|
38 | * |
---|
39 | * Modifications: |
---|
40 | * |
---|
41 | * 1. We only require IEEE, IBM, or VAX double-precision |
---|
42 | * arithmetic (not IEEE double-extended). |
---|
43 | * 2. We get by with floating-point arithmetic in a case that |
---|
44 | * Clinger missed -- when we're computing d * 10^n |
---|
45 | * for a small integer d and the integer n is not too |
---|
46 | * much larger than 22 (the maximum integer k for which |
---|
47 | * we can represent 10^k exactly), we may be able to |
---|
48 | * compute (d*10^k) * 10^(e-k) with just one roundoff. |
---|
49 | * 3. Rather than a bit-at-a-time adjustment of the binary |
---|
50 | * result in the hard case, we use floating-point |
---|
51 | * arithmetic to determine the adjustment to within |
---|
52 | * one bit; only in really hard cases do we need to |
---|
53 | * compute a second residual. |
---|
54 | * 4. Because of 3., we don't need a large table of powers of 10 |
---|
55 | * for ten-to-e (just some small tables, e.g. of 10^k |
---|
56 | * for 0 <= k <= 22). |
---|
57 | */ |
---|
58 | |
---|
59 | /* |
---|
60 | * #define IEEE_8087 for IEEE-arithmetic machines where the least |
---|
61 | * significant byte has the lowest address. |
---|
62 | * #define IEEE_MC68k for IEEE-arithmetic machines where the most |
---|
63 | * significant byte has the lowest address. |
---|
64 | * #define Sudden_Underflow for IEEE-format machines without gradual |
---|
65 | * underflow (i.e., that flush to zero on underflow). |
---|
66 | * #define IBM for IBM mainframe-style floating-point arithmetic. |
---|
67 | * #define VAX for VAX-style floating-point arithmetic. |
---|
68 | * #define Unsigned_Shifts if >> does treats its left operand as unsigned. |
---|
69 | * #define No_leftright to omit left-right logic in fast floating-point |
---|
70 | * computation of dtoa. |
---|
71 | * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3. |
---|
72 | * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines |
---|
73 | * that use extended-precision instructions to compute rounded |
---|
74 | * products and quotients) with IBM. |
---|
75 | * #define ROUND_BIASED for IEEE-format with biased rounding. |
---|
76 | * #define Inaccurate_Divide for IEEE-format with correctly rounded |
---|
77 | * products but inaccurate quotients, e.g., for Intel i860. |
---|
78 | * #define Just_16 to store 16 bits per 32-bit long when doing high-precision |
---|
79 | * integer arithmetic. Whether this speeds things up or slows things |
---|
80 | * down depends on the machine and the number being converted. |
---|
81 | */ |
---|
82 | |
---|
83 | #include <_ansi.h> |
---|
84 | #include <stdlib.h> |
---|
85 | #include <string.h> |
---|
86 | #include <reent.h> |
---|
87 | #include "mprec.h" |
---|
88 | |
---|
89 | /* This is defined in sys/reent.h as (sizeof (size_t) << 3) now, as in NetBSD. |
---|
90 | The old value of 15 was wrong and made newlib vulnerable against buffer |
---|
91 | overrun attacks (CVE-2009-0689), same as other implementations of gdtoa |
---|
92 | based on BSD code. |
---|
93 | #define _Kmax 15 |
---|
94 | */ |
---|
95 | |
---|
96 | _Bigint * |
---|
97 | Balloc (struct _reent *ptr, int k) |
---|
98 | { |
---|
99 | int x; |
---|
100 | _Bigint *rv ; |
---|
101 | |
---|
102 | _REENT_CHECK_MP(ptr); |
---|
103 | if (_REENT_MP_FREELIST(ptr) == NULL) |
---|
104 | { |
---|
105 | /* Allocate a list of pointers to the mprec objects */ |
---|
106 | _REENT_MP_FREELIST(ptr) = (struct _Bigint **) _calloc_r (ptr, |
---|
107 | sizeof (struct _Bigint *), |
---|
108 | _Kmax + 1); |
---|
109 | if (_REENT_MP_FREELIST(ptr) == NULL) |
---|
110 | { |
---|
111 | return NULL; |
---|
112 | } |
---|
113 | } |
---|
114 | |
---|
115 | if ((rv = _REENT_MP_FREELIST(ptr)[k]) != 0) |
---|
116 | { |
---|
117 | _REENT_MP_FREELIST(ptr)[k] = rv->_next; |
---|
118 | } |
---|
119 | else |
---|
120 | { |
---|
121 | x = 1 << k; |
---|
122 | /* Allocate an mprec Bigint and stick in in the freelist */ |
---|
123 | rv = (_Bigint *) _calloc_r (ptr, |
---|
124 | 1, |
---|
125 | sizeof (_Bigint) + |
---|
126 | (x-1) * sizeof(rv->_x)); |
---|
127 | if (rv == NULL) return NULL; |
---|
128 | rv->_k = k; |
---|
129 | rv->_maxwds = x; |
---|
130 | } |
---|
131 | rv->_sign = rv->_wds = 0; |
---|
132 | return rv; |
---|
133 | } |
---|
134 | |
---|
135 | void |
---|
136 | Bfree (struct _reent *ptr, _Bigint * v) |
---|
137 | { |
---|
138 | _REENT_CHECK_MP(ptr); |
---|
139 | if (v) |
---|
140 | { |
---|
141 | v->_next = _REENT_MP_FREELIST(ptr)[v->_k]; |
---|
142 | _REENT_MP_FREELIST(ptr)[v->_k] = v; |
---|
143 | } |
---|
144 | } |
---|
145 | |
---|
146 | _Bigint * |
---|
147 | multadd (struct _reent *ptr, |
---|
148 | _Bigint * b, |
---|
149 | int m, |
---|
150 | int a) |
---|
151 | { |
---|
152 | int i, wds; |
---|
153 | __ULong *x, y; |
---|
154 | #ifdef Pack_32 |
---|
155 | __ULong xi, z; |
---|
156 | #endif |
---|
157 | _Bigint *b1; |
---|
158 | |
---|
159 | wds = b->_wds; |
---|
160 | x = b->_x; |
---|
161 | i = 0; |
---|
162 | do |
---|
163 | { |
---|
164 | #ifdef Pack_32 |
---|
165 | xi = *x; |
---|
166 | y = (xi & 0xffff) * m + a; |
---|
167 | z = (xi >> 16) * m + (y >> 16); |
---|
168 | a = (int) (z >> 16); |
---|
169 | *x++ = (z << 16) + (y & 0xffff); |
---|
170 | #else |
---|
171 | y = *x * m + a; |
---|
172 | a = (int) (y >> 16); |
---|
173 | *x++ = y & 0xffff; |
---|
174 | #endif |
---|
175 | } |
---|
176 | while (++i < wds); |
---|
177 | if (a) |
---|
178 | { |
---|
179 | if (wds >= b->_maxwds) |
---|
180 | { |
---|
181 | b1 = Balloc (ptr, b->_k + 1); |
---|
182 | Bcopy (b1, b); |
---|
183 | Bfree (ptr, b); |
---|
184 | b = b1; |
---|
185 | } |
---|
186 | b->_x[wds++] = a; |
---|
187 | b->_wds = wds; |
---|
188 | } |
---|
189 | return b; |
---|
190 | } |
---|
191 | |
---|
192 | _Bigint * |
---|
193 | s2b (struct _reent * ptr, |
---|
194 | const char *s, |
---|
195 | int nd0, |
---|
196 | int nd, |
---|
197 | __ULong y9) |
---|
198 | { |
---|
199 | _Bigint *b; |
---|
200 | int i, k; |
---|
201 | __Long x, y; |
---|
202 | |
---|
203 | x = (nd + 8) / 9; |
---|
204 | for (k = 0, y = 1; x > y; y <<= 1, k++); |
---|
205 | #ifdef Pack_32 |
---|
206 | b = Balloc (ptr, k); |
---|
207 | b->_x[0] = y9; |
---|
208 | b->_wds = 1; |
---|
209 | #else |
---|
210 | b = Balloc (ptr, k + 1); |
---|
211 | b->_x[0] = y9 & 0xffff; |
---|
212 | b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1; |
---|
213 | #endif |
---|
214 | |
---|
215 | i = 9; |
---|
216 | if (9 < nd0) |
---|
217 | { |
---|
218 | s += 9; |
---|
219 | do |
---|
220 | b = multadd (ptr, b, 10, *s++ - '0'); |
---|
221 | while (++i < nd0); |
---|
222 | s++; |
---|
223 | } |
---|
224 | else |
---|
225 | s += 10; |
---|
226 | for (; i < nd; i++) |
---|
227 | b = multadd (ptr, b, 10, *s++ - '0'); |
---|
228 | return b; |
---|
229 | } |
---|
230 | |
---|
231 | int |
---|
232 | hi0bits (register __ULong x) |
---|
233 | { |
---|
234 | register int k = 0; |
---|
235 | |
---|
236 | if (!(x & 0xffff0000)) |
---|
237 | { |
---|
238 | k = 16; |
---|
239 | x <<= 16; |
---|
240 | } |
---|
241 | if (!(x & 0xff000000)) |
---|
242 | { |
---|
243 | k += 8; |
---|
244 | x <<= 8; |
---|
245 | } |
---|
246 | if (!(x & 0xf0000000)) |
---|
247 | { |
---|
248 | k += 4; |
---|
249 | x <<= 4; |
---|
250 | } |
---|
251 | if (!(x & 0xc0000000)) |
---|
252 | { |
---|
253 | k += 2; |
---|
254 | x <<= 2; |
---|
255 | } |
---|
256 | if (!(x & 0x80000000)) |
---|
257 | { |
---|
258 | k++; |
---|
259 | if (!(x & 0x40000000)) |
---|
260 | return 32; |
---|
261 | } |
---|
262 | return k; |
---|
263 | } |
---|
264 | |
---|
265 | int |
---|
266 | lo0bits (__ULong *y) |
---|
267 | { |
---|
268 | register int k; |
---|
269 | register __ULong x = *y; |
---|
270 | |
---|
271 | if (x & 7) |
---|
272 | { |
---|
273 | if (x & 1) |
---|
274 | return 0; |
---|
275 | if (x & 2) |
---|
276 | { |
---|
277 | *y = x >> 1; |
---|
278 | return 1; |
---|
279 | } |
---|
280 | *y = x >> 2; |
---|
281 | return 2; |
---|
282 | } |
---|
283 | k = 0; |
---|
284 | if (!(x & 0xffff)) |
---|
285 | { |
---|
286 | k = 16; |
---|
287 | x >>= 16; |
---|
288 | } |
---|
289 | if (!(x & 0xff)) |
---|
290 | { |
---|
291 | k += 8; |
---|
292 | x >>= 8; |
---|
293 | } |
---|
294 | if (!(x & 0xf)) |
---|
295 | { |
---|
296 | k += 4; |
---|
297 | x >>= 4; |
---|
298 | } |
---|
299 | if (!(x & 0x3)) |
---|
300 | { |
---|
301 | k += 2; |
---|
302 | x >>= 2; |
---|
303 | } |
---|
304 | if (!(x & 1)) |
---|
305 | { |
---|
306 | k++; |
---|
307 | x >>= 1; |
---|
308 | if (!x & 1) |
---|
309 | return 32; |
---|
310 | } |
---|
311 | *y = x; |
---|
312 | return k; |
---|
313 | } |
---|
314 | |
---|
315 | _Bigint * |
---|
316 | i2b (struct _reent * ptr, int i) |
---|
317 | { |
---|
318 | _Bigint *b; |
---|
319 | |
---|
320 | b = Balloc (ptr, 1); |
---|
321 | b->_x[0] = i; |
---|
322 | b->_wds = 1; |
---|
323 | return b; |
---|
324 | } |
---|
325 | |
---|
326 | _Bigint * |
---|
327 | mult (struct _reent * ptr, _Bigint * a, _Bigint * b) |
---|
328 | { |
---|
329 | _Bigint *c; |
---|
330 | int k, wa, wb, wc; |
---|
331 | __ULong carry, y, z; |
---|
332 | __ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; |
---|
333 | #ifdef Pack_32 |
---|
334 | __ULong z2; |
---|
335 | #endif |
---|
336 | |
---|
337 | if (a->_wds < b->_wds) |
---|
338 | { |
---|
339 | c = a; |
---|
340 | a = b; |
---|
341 | b = c; |
---|
342 | } |
---|
343 | k = a->_k; |
---|
344 | wa = a->_wds; |
---|
345 | wb = b->_wds; |
---|
346 | wc = wa + wb; |
---|
347 | if (wc > a->_maxwds) |
---|
348 | k++; |
---|
349 | c = Balloc (ptr, k); |
---|
350 | for (x = c->_x, xa = x + wc; x < xa; x++) |
---|
351 | *x = 0; |
---|
352 | xa = a->_x; |
---|
353 | xae = xa + wa; |
---|
354 | xb = b->_x; |
---|
355 | xbe = xb + wb; |
---|
356 | xc0 = c->_x; |
---|
357 | #ifdef Pack_32 |
---|
358 | for (; xb < xbe; xb++, xc0++) |
---|
359 | { |
---|
360 | if ((y = *xb & 0xffff) != 0) |
---|
361 | { |
---|
362 | x = xa; |
---|
363 | xc = xc0; |
---|
364 | carry = 0; |
---|
365 | do |
---|
366 | { |
---|
367 | z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; |
---|
368 | carry = z >> 16; |
---|
369 | z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; |
---|
370 | carry = z2 >> 16; |
---|
371 | Storeinc (xc, z2, z); |
---|
372 | } |
---|
373 | while (x < xae); |
---|
374 | *xc = carry; |
---|
375 | } |
---|
376 | if ((y = *xb >> 16) != 0) |
---|
377 | { |
---|
378 | x = xa; |
---|
379 | xc = xc0; |
---|
380 | carry = 0; |
---|
381 | z2 = *xc; |
---|
382 | do |
---|
383 | { |
---|
384 | z = (*x & 0xffff) * y + (*xc >> 16) + carry; |
---|
385 | carry = z >> 16; |
---|
386 | Storeinc (xc, z, z2); |
---|
387 | z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; |
---|
388 | carry = z2 >> 16; |
---|
389 | } |
---|
390 | while (x < xae); |
---|
391 | *xc = z2; |
---|
392 | } |
---|
393 | } |
---|
394 | #else |
---|
395 | for (; xb < xbe; xc0++) |
---|
396 | { |
---|
397 | if (y = *xb++) |
---|
398 | { |
---|
399 | x = xa; |
---|
400 | xc = xc0; |
---|
401 | carry = 0; |
---|
402 | do |
---|
403 | { |
---|
404 | z = *x++ * y + *xc + carry; |
---|
405 | carry = z >> 16; |
---|
406 | *xc++ = z & 0xffff; |
---|
407 | } |
---|
408 | while (x < xae); |
---|
409 | *xc = carry; |
---|
410 | } |
---|
411 | } |
---|
412 | #endif |
---|
413 | for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc); |
---|
414 | c->_wds = wc; |
---|
415 | return c; |
---|
416 | } |
---|
417 | |
---|
418 | _Bigint * |
---|
419 | pow5mult (struct _reent * ptr, _Bigint * b, int k) |
---|
420 | { |
---|
421 | _Bigint *b1, *p5, *p51; |
---|
422 | int i; |
---|
423 | static const int p05[3] = {5, 25, 125}; |
---|
424 | |
---|
425 | if ((i = k & 3) != 0) |
---|
426 | b = multadd (ptr, b, p05[i - 1], 0); |
---|
427 | |
---|
428 | if (!(k >>= 2)) |
---|
429 | return b; |
---|
430 | _REENT_CHECK_MP(ptr); |
---|
431 | if (!(p5 = _REENT_MP_P5S(ptr))) |
---|
432 | { |
---|
433 | /* first time */ |
---|
434 | p5 = _REENT_MP_P5S(ptr) = i2b (ptr, 625); |
---|
435 | p5->_next = 0; |
---|
436 | } |
---|
437 | for (;;) |
---|
438 | { |
---|
439 | if (k & 1) |
---|
440 | { |
---|
441 | b1 = mult (ptr, b, p5); |
---|
442 | Bfree (ptr, b); |
---|
443 | b = b1; |
---|
444 | } |
---|
445 | if (!(k >>= 1)) |
---|
446 | break; |
---|
447 | if (!(p51 = p5->_next)) |
---|
448 | { |
---|
449 | p51 = p5->_next = mult (ptr, p5, p5); |
---|
450 | p51->_next = 0; |
---|
451 | } |
---|
452 | p5 = p51; |
---|
453 | } |
---|
454 | return b; |
---|
455 | } |
---|
456 | |
---|
457 | _Bigint * |
---|
458 | lshift (struct _reent * ptr, _Bigint * b, int k) |
---|
459 | { |
---|
460 | int i, k1, n, n1; |
---|
461 | _Bigint *b1; |
---|
462 | __ULong *x, *x1, *xe, z; |
---|
463 | |
---|
464 | #ifdef Pack_32 |
---|
465 | n = k >> 5; |
---|
466 | #else |
---|
467 | n = k >> 4; |
---|
468 | #endif |
---|
469 | k1 = b->_k; |
---|
470 | n1 = n + b->_wds + 1; |
---|
471 | for (i = b->_maxwds; n1 > i; i <<= 1) |
---|
472 | k1++; |
---|
473 | b1 = Balloc (ptr, k1); |
---|
474 | x1 = b1->_x; |
---|
475 | for (i = 0; i < n; i++) |
---|
476 | *x1++ = 0; |
---|
477 | x = b->_x; |
---|
478 | xe = x + b->_wds; |
---|
479 | #ifdef Pack_32 |
---|
480 | if (k &= 0x1f) |
---|
481 | { |
---|
482 | k1 = 32 - k; |
---|
483 | z = 0; |
---|
484 | do |
---|
485 | { |
---|
486 | *x1++ = *x << k | z; |
---|
487 | z = *x++ >> k1; |
---|
488 | } |
---|
489 | while (x < xe); |
---|
490 | if ((*x1 = z) != 0) |
---|
491 | ++n1; |
---|
492 | } |
---|
493 | #else |
---|
494 | if (k &= 0xf) |
---|
495 | { |
---|
496 | k1 = 16 - k; |
---|
497 | z = 0; |
---|
498 | do |
---|
499 | { |
---|
500 | *x1++ = *x << k & 0xffff | z; |
---|
501 | z = *x++ >> k1; |
---|
502 | } |
---|
503 | while (x < xe); |
---|
504 | if (*x1 = z) |
---|
505 | ++n1; |
---|
506 | } |
---|
507 | #endif |
---|
508 | else |
---|
509 | do |
---|
510 | *x1++ = *x++; |
---|
511 | while (x < xe); |
---|
512 | b1->_wds = n1 - 1; |
---|
513 | Bfree (ptr, b); |
---|
514 | return b1; |
---|
515 | } |
---|
516 | |
---|
517 | int |
---|
518 | cmp (_Bigint * a, _Bigint * b) |
---|
519 | { |
---|
520 | __ULong *xa, *xa0, *xb, *xb0; |
---|
521 | int i, j; |
---|
522 | |
---|
523 | i = a->_wds; |
---|
524 | j = b->_wds; |
---|
525 | #ifdef DEBUG |
---|
526 | if (i > 1 && !a->_x[i - 1]) |
---|
527 | Bug ("cmp called with a->_x[a->_wds-1] == 0"); |
---|
528 | if (j > 1 && !b->_x[j - 1]) |
---|
529 | Bug ("cmp called with b->_x[b->_wds-1] == 0"); |
---|
530 | #endif |
---|
531 | if (i -= j) |
---|
532 | return i; |
---|
533 | xa0 = a->_x; |
---|
534 | xa = xa0 + j; |
---|
535 | xb0 = b->_x; |
---|
536 | xb = xb0 + j; |
---|
537 | for (;;) |
---|
538 | { |
---|
539 | if (*--xa != *--xb) |
---|
540 | return *xa < *xb ? -1 : 1; |
---|
541 | if (xa <= xa0) |
---|
542 | break; |
---|
543 | } |
---|
544 | return 0; |
---|
545 | } |
---|
546 | |
---|
547 | _Bigint * |
---|
548 | diff (struct _reent * ptr, |
---|
549 | _Bigint * a, _Bigint * b) |
---|
550 | { |
---|
551 | _Bigint *c; |
---|
552 | int i, wa, wb; |
---|
553 | __Long borrow, y; /* We need signed shifts here. */ |
---|
554 | __ULong *xa, *xae, *xb, *xbe, *xc; |
---|
555 | #ifdef Pack_32 |
---|
556 | __Long z; |
---|
557 | #endif |
---|
558 | |
---|
559 | i = cmp (a, b); |
---|
560 | if (!i) |
---|
561 | { |
---|
562 | c = Balloc (ptr, 0); |
---|
563 | c->_wds = 1; |
---|
564 | c->_x[0] = 0; |
---|
565 | return c; |
---|
566 | } |
---|
567 | if (i < 0) |
---|
568 | { |
---|
569 | c = a; |
---|
570 | a = b; |
---|
571 | b = c; |
---|
572 | i = 1; |
---|
573 | } |
---|
574 | else |
---|
575 | i = 0; |
---|
576 | c = Balloc (ptr, a->_k); |
---|
577 | c->_sign = i; |
---|
578 | wa = a->_wds; |
---|
579 | xa = a->_x; |
---|
580 | xae = xa + wa; |
---|
581 | wb = b->_wds; |
---|
582 | xb = b->_x; |
---|
583 | xbe = xb + wb; |
---|
584 | xc = c->_x; |
---|
585 | borrow = 0; |
---|
586 | #ifdef Pack_32 |
---|
587 | do |
---|
588 | { |
---|
589 | y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; |
---|
590 | borrow = y >> 16; |
---|
591 | Sign_Extend (borrow, y); |
---|
592 | z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; |
---|
593 | borrow = z >> 16; |
---|
594 | Sign_Extend (borrow, z); |
---|
595 | Storeinc (xc, z, y); |
---|
596 | } |
---|
597 | while (xb < xbe); |
---|
598 | while (xa < xae) |
---|
599 | { |
---|
600 | y = (*xa & 0xffff) + borrow; |
---|
601 | borrow = y >> 16; |
---|
602 | Sign_Extend (borrow, y); |
---|
603 | z = (*xa++ >> 16) + borrow; |
---|
604 | borrow = z >> 16; |
---|
605 | Sign_Extend (borrow, z); |
---|
606 | Storeinc (xc, z, y); |
---|
607 | } |
---|
608 | #else |
---|
609 | do |
---|
610 | { |
---|
611 | y = *xa++ - *xb++ + borrow; |
---|
612 | borrow = y >> 16; |
---|
613 | Sign_Extend (borrow, y); |
---|
614 | *xc++ = y & 0xffff; |
---|
615 | } |
---|
616 | while (xb < xbe); |
---|
617 | while (xa < xae) |
---|
618 | { |
---|
619 | y = *xa++ + borrow; |
---|
620 | borrow = y >> 16; |
---|
621 | Sign_Extend (borrow, y); |
---|
622 | *xc++ = y & 0xffff; |
---|
623 | } |
---|
624 | #endif |
---|
625 | while (!*--xc) |
---|
626 | wa--; |
---|
627 | c->_wds = wa; |
---|
628 | return c; |
---|
629 | } |
---|
630 | |
---|
631 | double |
---|
632 | ulp (double _x) |
---|
633 | { |
---|
634 | union double_union x, a; |
---|
635 | register __Long L; |
---|
636 | |
---|
637 | x.d = _x; |
---|
638 | |
---|
639 | L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1; |
---|
640 | #ifndef Sudden_Underflow |
---|
641 | if (L > 0) |
---|
642 | { |
---|
643 | #endif |
---|
644 | #ifdef IBM |
---|
645 | L |= Exp_msk1 >> 4; |
---|
646 | #endif |
---|
647 | word0 (a) = L; |
---|
648 | #ifndef _DOUBLE_IS_32BITS |
---|
649 | word1 (a) = 0; |
---|
650 | #endif |
---|
651 | |
---|
652 | #ifndef Sudden_Underflow |
---|
653 | } |
---|
654 | else |
---|
655 | { |
---|
656 | L = -L >> Exp_shift; |
---|
657 | if (L < Exp_shift) |
---|
658 | { |
---|
659 | word0 (a) = 0x80000 >> L; |
---|
660 | #ifndef _DOUBLE_IS_32BITS |
---|
661 | word1 (a) = 0; |
---|
662 | #endif |
---|
663 | } |
---|
664 | else |
---|
665 | { |
---|
666 | word0 (a) = 0; |
---|
667 | L -= Exp_shift; |
---|
668 | #ifndef _DOUBLE_IS_32BITS |
---|
669 | word1 (a) = L >= 31 ? 1 : 1 << (31 - L); |
---|
670 | #endif |
---|
671 | } |
---|
672 | } |
---|
673 | #endif |
---|
674 | return a.d; |
---|
675 | } |
---|
676 | |
---|
677 | double |
---|
678 | b2d (_Bigint * a, int *e) |
---|
679 | { |
---|
680 | __ULong *xa, *xa0, w, y, z; |
---|
681 | int k; |
---|
682 | union double_union d; |
---|
683 | #ifdef VAX |
---|
684 | __ULong d0, d1; |
---|
685 | #else |
---|
686 | #define d0 word0(d) |
---|
687 | #define d1 word1(d) |
---|
688 | #endif |
---|
689 | |
---|
690 | xa0 = a->_x; |
---|
691 | xa = xa0 + a->_wds; |
---|
692 | y = *--xa; |
---|
693 | #ifdef DEBUG |
---|
694 | if (!y) |
---|
695 | Bug ("zero y in b2d"); |
---|
696 | #endif |
---|
697 | k = hi0bits (y); |
---|
698 | *e = 32 - k; |
---|
699 | #ifdef Pack_32 |
---|
700 | if (k < Ebits) |
---|
701 | { |
---|
702 | d0 = Exp_1 | y >> (Ebits - k); |
---|
703 | w = xa > xa0 ? *--xa : 0; |
---|
704 | #ifndef _DOUBLE_IS_32BITS |
---|
705 | d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k); |
---|
706 | #endif |
---|
707 | goto ret_d; |
---|
708 | } |
---|
709 | z = xa > xa0 ? *--xa : 0; |
---|
710 | if (k -= Ebits) |
---|
711 | { |
---|
712 | d0 = Exp_1 | y << k | z >> (32 - k); |
---|
713 | y = xa > xa0 ? *--xa : 0; |
---|
714 | #ifndef _DOUBLE_IS_32BITS |
---|
715 | d1 = z << k | y >> (32 - k); |
---|
716 | #endif |
---|
717 | } |
---|
718 | else |
---|
719 | { |
---|
720 | d0 = Exp_1 | y; |
---|
721 | #ifndef _DOUBLE_IS_32BITS |
---|
722 | d1 = z; |
---|
723 | #endif |
---|
724 | } |
---|
725 | #else |
---|
726 | if (k < Ebits + 16) |
---|
727 | { |
---|
728 | z = xa > xa0 ? *--xa : 0; |
---|
729 | d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; |
---|
730 | w = xa > xa0 ? *--xa : 0; |
---|
731 | y = xa > xa0 ? *--xa : 0; |
---|
732 | d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; |
---|
733 | goto ret_d; |
---|
734 | } |
---|
735 | z = xa > xa0 ? *--xa : 0; |
---|
736 | w = xa > xa0 ? *--xa : 0; |
---|
737 | k -= Ebits + 16; |
---|
738 | d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; |
---|
739 | y = xa > xa0 ? *--xa : 0; |
---|
740 | d1 = w << k + 16 | y << k; |
---|
741 | #endif |
---|
742 | ret_d: |
---|
743 | #ifdef VAX |
---|
744 | word0 (d) = d0 >> 16 | d0 << 16; |
---|
745 | word1 (d) = d1 >> 16 | d1 << 16; |
---|
746 | #else |
---|
747 | #undef d0 |
---|
748 | #undef d1 |
---|
749 | #endif |
---|
750 | return d.d; |
---|
751 | } |
---|
752 | |
---|
753 | _Bigint * |
---|
754 | d2b (struct _reent * ptr, |
---|
755 | double _d, |
---|
756 | int *e, |
---|
757 | int *bits) |
---|
758 | |
---|
759 | { |
---|
760 | union double_union d; |
---|
761 | _Bigint *b; |
---|
762 | int de, i, k; |
---|
763 | __ULong *x, y, z; |
---|
764 | #ifdef VAX |
---|
765 | __ULong d0, d1; |
---|
766 | #endif |
---|
767 | d.d = _d; |
---|
768 | #ifdef VAX |
---|
769 | d0 = word0 (d) >> 16 | word0 (d) << 16; |
---|
770 | d1 = word1 (d) >> 16 | word1 (d) << 16; |
---|
771 | #else |
---|
772 | #define d0 word0(d) |
---|
773 | #define d1 word1(d) |
---|
774 | d.d = _d; |
---|
775 | #endif |
---|
776 | |
---|
777 | #ifdef Pack_32 |
---|
778 | b = Balloc (ptr, 1); |
---|
779 | #else |
---|
780 | b = Balloc (ptr, 2); |
---|
781 | #endif |
---|
782 | x = b->_x; |
---|
783 | |
---|
784 | z = d0 & Frac_mask; |
---|
785 | d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ |
---|
786 | #ifdef Sudden_Underflow |
---|
787 | de = (int) (d0 >> Exp_shift); |
---|
788 | #ifndef IBM |
---|
789 | z |= Exp_msk11; |
---|
790 | #endif |
---|
791 | #else |
---|
792 | if ((de = (int) (d0 >> Exp_shift)) != 0) |
---|
793 | z |= Exp_msk1; |
---|
794 | #endif |
---|
795 | #ifdef Pack_32 |
---|
796 | #ifndef _DOUBLE_IS_32BITS |
---|
797 | if (d1) |
---|
798 | { |
---|
799 | y = d1; |
---|
800 | k = lo0bits (&y); |
---|
801 | if (k) |
---|
802 | { |
---|
803 | x[0] = y | z << (32 - k); |
---|
804 | z >>= k; |
---|
805 | } |
---|
806 | else |
---|
807 | x[0] = y; |
---|
808 | i = b->_wds = (x[1] = z) ? 2 : 1; |
---|
809 | } |
---|
810 | else |
---|
811 | #endif |
---|
812 | { |
---|
813 | #ifdef DEBUG |
---|
814 | if (!z) |
---|
815 | Bug ("Zero passed to d2b"); |
---|
816 | #endif |
---|
817 | k = lo0bits (&z); |
---|
818 | x[0] = z; |
---|
819 | i = b->_wds = 1; |
---|
820 | #ifndef _DOUBLE_IS_32BITS |
---|
821 | k += 32; |
---|
822 | #endif |
---|
823 | } |
---|
824 | #else |
---|
825 | if (d1) |
---|
826 | { |
---|
827 | y = d1; |
---|
828 | k = lo0bits (&y); |
---|
829 | if (k) |
---|
830 | if (k >= 16) |
---|
831 | { |
---|
832 | x[0] = y | z << 32 - k & 0xffff; |
---|
833 | x[1] = z >> k - 16 & 0xffff; |
---|
834 | x[2] = z >> k; |
---|
835 | i = 2; |
---|
836 | } |
---|
837 | else |
---|
838 | { |
---|
839 | x[0] = y & 0xffff; |
---|
840 | x[1] = y >> 16 | z << 16 - k & 0xffff; |
---|
841 | x[2] = z >> k & 0xffff; |
---|
842 | x[3] = z >> k + 16; |
---|
843 | i = 3; |
---|
844 | } |
---|
845 | else |
---|
846 | { |
---|
847 | x[0] = y & 0xffff; |
---|
848 | x[1] = y >> 16; |
---|
849 | x[2] = z & 0xffff; |
---|
850 | x[3] = z >> 16; |
---|
851 | i = 3; |
---|
852 | } |
---|
853 | } |
---|
854 | else |
---|
855 | { |
---|
856 | #ifdef DEBUG |
---|
857 | if (!z) |
---|
858 | Bug ("Zero passed to d2b"); |
---|
859 | #endif |
---|
860 | k = lo0bits (&z); |
---|
861 | if (k >= 16) |
---|
862 | { |
---|
863 | x[0] = z; |
---|
864 | i = 0; |
---|
865 | } |
---|
866 | else |
---|
867 | { |
---|
868 | x[0] = z & 0xffff; |
---|
869 | x[1] = z >> 16; |
---|
870 | i = 1; |
---|
871 | } |
---|
872 | k += 32; |
---|
873 | } |
---|
874 | while (!x[i]) |
---|
875 | --i; |
---|
876 | b->_wds = i + 1; |
---|
877 | #endif |
---|
878 | #ifndef Sudden_Underflow |
---|
879 | if (de) |
---|
880 | { |
---|
881 | #endif |
---|
882 | #ifdef IBM |
---|
883 | *e = (de - Bias - (P - 1) << 2) + k; |
---|
884 | *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask); |
---|
885 | #else |
---|
886 | *e = de - Bias - (P - 1) + k; |
---|
887 | *bits = P - k; |
---|
888 | #endif |
---|
889 | #ifndef Sudden_Underflow |
---|
890 | } |
---|
891 | else |
---|
892 | { |
---|
893 | *e = de - Bias - (P - 1) + 1 + k; |
---|
894 | #ifdef Pack_32 |
---|
895 | *bits = 32 * i - hi0bits (x[i - 1]); |
---|
896 | #else |
---|
897 | *bits = (i + 2) * 16 - hi0bits (x[i]); |
---|
898 | #endif |
---|
899 | } |
---|
900 | #endif |
---|
901 | return b; |
---|
902 | } |
---|
903 | #undef d0 |
---|
904 | #undef d1 |
---|
905 | |
---|
906 | double |
---|
907 | ratio (_Bigint * a, _Bigint * b) |
---|
908 | |
---|
909 | { |
---|
910 | union double_union da, db; |
---|
911 | int k, ka, kb; |
---|
912 | |
---|
913 | da.d = b2d (a, &ka); |
---|
914 | db.d = b2d (b, &kb); |
---|
915 | #ifdef Pack_32 |
---|
916 | k = ka - kb + 32 * (a->_wds - b->_wds); |
---|
917 | #else |
---|
918 | k = ka - kb + 16 * (a->_wds - b->_wds); |
---|
919 | #endif |
---|
920 | #ifdef IBM |
---|
921 | if (k > 0) |
---|
922 | { |
---|
923 | word0 (da) += (k >> 2) * Exp_msk1; |
---|
924 | if (k &= 3) |
---|
925 | da.d *= 1 << k; |
---|
926 | } |
---|
927 | else |
---|
928 | { |
---|
929 | k = -k; |
---|
930 | word0 (db) += (k >> 2) * Exp_msk1; |
---|
931 | if (k &= 3) |
---|
932 | db.d *= 1 << k; |
---|
933 | } |
---|
934 | #else |
---|
935 | if (k > 0) |
---|
936 | word0 (da) += k * Exp_msk1; |
---|
937 | else |
---|
938 | { |
---|
939 | k = -k; |
---|
940 | word0 (db) += k * Exp_msk1; |
---|
941 | } |
---|
942 | #endif |
---|
943 | return da.d / db.d; |
---|
944 | } |
---|
945 | |
---|
946 | |
---|
947 | const double |
---|
948 | tens[] = |
---|
949 | { |
---|
950 | 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, |
---|
951 | 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, |
---|
952 | 1e20, 1e21, 1e22, 1e23, 1e24 |
---|
953 | |
---|
954 | }; |
---|
955 | |
---|
956 | #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800) |
---|
957 | const double bigtens[] = |
---|
958 | {1e16, 1e32, 1e64, 1e128, 1e256}; |
---|
959 | |
---|
960 | const double tinytens[] = |
---|
961 | {1e-16, 1e-32, 1e-64, 1e-128, 1e-256}; |
---|
962 | #else |
---|
963 | const double bigtens[] = |
---|
964 | {1e16, 1e32}; |
---|
965 | |
---|
966 | const double tinytens[] = |
---|
967 | {1e-16, 1e-32}; |
---|
968 | #endif |
---|
969 | |
---|
970 | |
---|
971 | double |
---|
972 | _mprec_log10 (int dig) |
---|
973 | { |
---|
974 | double v = 1.0; |
---|
975 | if (dig < 24) |
---|
976 | return tens[dig]; |
---|
977 | while (dig > 0) |
---|
978 | { |
---|
979 | v *= 10; |
---|
980 | dig--; |
---|
981 | } |
---|
982 | return v; |
---|
983 | } |
---|
984 | |
---|
985 | void |
---|
986 | copybits (__ULong *c, |
---|
987 | int n, |
---|
988 | _Bigint *b) |
---|
989 | { |
---|
990 | __ULong *ce, *x, *xe; |
---|
991 | #ifdef Pack_16 |
---|
992 | int nw, nw1; |
---|
993 | #endif |
---|
994 | |
---|
995 | ce = c + ((n-1) >> kshift) + 1; |
---|
996 | x = b->_x; |
---|
997 | #ifdef Pack_32 |
---|
998 | xe = x + b->_wds; |
---|
999 | while(x < xe) |
---|
1000 | *c++ = *x++; |
---|
1001 | #else |
---|
1002 | nw = b->_wds; |
---|
1003 | nw1 = nw & 1; |
---|
1004 | for(xe = x + (nw - nw1); x < xe; x += 2) |
---|
1005 | Storeinc(c, x[1], x[0]); |
---|
1006 | if (nw1) |
---|
1007 | *c++ = *x; |
---|
1008 | #endif |
---|
1009 | while(c < ce) |
---|
1010 | *c++ = 0; |
---|
1011 | } |
---|
1012 | |
---|
1013 | __ULong |
---|
1014 | any_on (_Bigint *b, |
---|
1015 | int k) |
---|
1016 | { |
---|
1017 | int n, nwds; |
---|
1018 | __ULong *x, *x0, x1, x2; |
---|
1019 | |
---|
1020 | x = b->_x; |
---|
1021 | nwds = b->_wds; |
---|
1022 | n = k >> kshift; |
---|
1023 | if (n > nwds) |
---|
1024 | n = nwds; |
---|
1025 | else if (n < nwds && (k &= kmask)) { |
---|
1026 | x1 = x2 = x[n]; |
---|
1027 | x1 >>= k; |
---|
1028 | x1 <<= k; |
---|
1029 | if (x1 != x2) |
---|
1030 | return 1; |
---|
1031 | } |
---|
1032 | x0 = x; |
---|
1033 | x += n; |
---|
1034 | while(x > x0) |
---|
1035 | if (*--x) |
---|
1036 | return 1; |
---|
1037 | return 0; |
---|
1038 | } |
---|
1039 | |
---|