]>
Commit | Line | Data |
---|---|---|
f911ba98 TT |
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 | ||
80f798b4 | 83 | /*#include <_ansi.h>*/ |
21276379 | 84 | #include <assert.h> |
f911ba98 TT |
85 | #include <stdlib.h> |
86 | #include <string.h> | |
80f798b4 | 87 | /* #include <reent.h> */ |
f911ba98 TT |
88 | #include "mprec.h" |
89 | ||
90 | /* reent.c knows this value */ | |
80f798b4 | 91 | /* #define _Kmax 15 */ |
f911ba98 | 92 | |
80f798b4 TT |
93 | #define _reent _Jv_reent |
94 | #define _Bigint _Jv_Bigint | |
f911ba98 | 95 | |
54fde020 | 96 | #undef _REENT_CHECK_MP |
80f798b4 | 97 | #define _REENT_CHECK_MP(x) |
54fde020 | 98 | #undef _REENT_MP_FREELIST |
80f798b4 | 99 | #define _REENT_MP_FREELIST(x) ((x)->_freelist) |
54fde020 | 100 | #undef _REENT_MP_P5S |
80f798b4 | 101 | #define _REENT_MP_P5S(x) ((x)->_p5s) |
f911ba98 | 102 | |
54fde020 BE |
103 | #undef __ULong |
104 | #define __ULong unsigned long | |
105 | #undef __Long | |
106 | #define __Long long | |
f911ba98 | 107 | |
80f798b4 | 108 | static void * |
6b61b957 | 109 | mprec_calloc (void *ignore, size_t x1, size_t x2) |
80f798b4 TT |
110 | { |
111 | char *result = (char *) malloc (x1 * x2); | |
112 | memset (result, 0, x1 * x2); | |
113 | return result; | |
114 | } | |
f911ba98 | 115 | |
80f798b4 TT |
116 | _Bigint * |
117 | _DEFUN (Balloc, (ptr, k), struct _reent *ptr _AND int k) | |
118 | { | |
119 | int x; | |
120 | _Bigint *rv ; | |
121 | int new_k = k + 1; | |
122 | ||
123 | _REENT_CHECK_MP(ptr); | |
124 | if (_REENT_MP_FREELIST(ptr) == NULL) | |
125 | { | |
126 | /* Allocate a list of pointers to the mprec objects */ | |
6b61b957 | 127 | _REENT_MP_FREELIST(ptr) = (struct _Bigint **) mprec_calloc (ptr, |
80f798b4 TT |
128 | sizeof (struct _Bigint *), |
129 | new_k); | |
130 | if (_REENT_MP_FREELIST(ptr) == NULL) | |
131 | { | |
132 | return NULL; | |
133 | } | |
134 | ptr->_max_k = new_k; | |
135 | } | |
136 | else if (new_k > ptr->_max_k) | |
137 | { | |
138 | struct _Bigint **new_list | |
139 | = (struct _Bigint **) realloc (ptr->_freelist, | |
140 | new_k * sizeof (struct _Bigint *)); | |
141 | memset (&new_list[ptr->_max_k], 0, | |
142 | (new_k - ptr->_max_k) * sizeof (struct _Bigint *)); | |
143 | ptr->_freelist = new_list; | |
144 | ptr->_max_k = new_k; | |
f911ba98 | 145 | |
80f798b4 | 146 | } |
f911ba98 | 147 | |
80f798b4 | 148 | assert (k <= ptr->_max_k); |
f911ba98 | 149 | |
80f798b4 TT |
150 | if ((rv = _REENT_MP_FREELIST(ptr)[k]) != 0) |
151 | { | |
152 | _REENT_MP_FREELIST(ptr)[k] = rv->_next; | |
153 | } | |
154 | else | |
155 | { | |
156 | x = 1 << k; | |
157 | /* Allocate an mprec Bigint and stick in in the freelist */ | |
6b61b957 | 158 | rv = (_Bigint *) mprec_calloc (ptr, |
80f798b4 TT |
159 | 1, |
160 | sizeof (_Bigint) + | |
161 | (x-1) * sizeof(rv->_x)); | |
162 | if (rv == NULL) return NULL; | |
163 | rv->_k = k; | |
164 | rv->_maxwds = x; | |
165 | } | |
166 | rv->_sign = rv->_wds = 0; | |
f911ba98 TT |
167 | return rv; |
168 | } | |
169 | ||
f911ba98 | 170 | void |
80f798b4 | 171 | _DEFUN (Bfree, (ptr, v), struct _reent *ptr _AND _Bigint * v) |
f911ba98 | 172 | { |
80f798b4 TT |
173 | _REENT_CHECK_MP(ptr); |
174 | if (v) | |
175 | { | |
176 | v->_next = _REENT_MP_FREELIST(ptr)[v->_k]; | |
177 | _REENT_MP_FREELIST(ptr)[v->_k] = v; | |
178 | } | |
f911ba98 TT |
179 | } |
180 | ||
80f798b4 | 181 | _Bigint * |
f911ba98 | 182 | _DEFUN (multadd, (ptr, b, m, a), |
80f798b4 TT |
183 | struct _reent *ptr _AND |
184 | _Bigint * b _AND | |
f911ba98 TT |
185 | int m _AND |
186 | int a) | |
187 | { | |
188 | int i, wds; | |
80f798b4 | 189 | __ULong *x, y; |
f911ba98 | 190 | #ifdef Pack_32 |
80f798b4 | 191 | __ULong xi, z; |
f911ba98 | 192 | #endif |
80f798b4 | 193 | _Bigint *b1; |
f911ba98 TT |
194 | |
195 | wds = b->_wds; | |
196 | x = b->_x; | |
197 | i = 0; | |
198 | do | |
199 | { | |
200 | #ifdef Pack_32 | |
201 | xi = *x; | |
202 | y = (xi & 0xffff) * m + a; | |
203 | z = (xi >> 16) * m + (y >> 16); | |
204 | a = (int) (z >> 16); | |
205 | *x++ = (z << 16) + (y & 0xffff); | |
206 | #else | |
207 | y = *x * m + a; | |
208 | a = (int) (y >> 16); | |
209 | *x++ = y & 0xffff; | |
210 | #endif | |
211 | } | |
212 | while (++i < wds); | |
213 | if (a) | |
214 | { | |
215 | if (wds >= b->_maxwds) | |
216 | { | |
217 | b1 = Balloc (ptr, b->_k + 1); | |
218 | Bcopy (b1, b); | |
219 | Bfree (ptr, b); | |
220 | b = b1; | |
221 | } | |
222 | b->_x[wds++] = a; | |
223 | b->_wds = wds; | |
224 | } | |
225 | return b; | |
226 | } | |
227 | ||
80f798b4 | 228 | _Bigint * |
f911ba98 | 229 | _DEFUN (s2b, (ptr, s, nd0, nd, y9), |
80f798b4 | 230 | struct _reent * ptr _AND |
f911ba98 TT |
231 | _CONST char *s _AND |
232 | int nd0 _AND | |
233 | int nd _AND | |
80f798b4 | 234 | __ULong y9) |
f911ba98 | 235 | { |
80f798b4 | 236 | _Bigint *b; |
f911ba98 | 237 | int i, k; |
80f798b4 | 238 | __Long x, y; |
f911ba98 TT |
239 | |
240 | x = (nd + 8) / 9; | |
241 | for (k = 0, y = 1; x > y; y <<= 1, k++); | |
242 | #ifdef Pack_32 | |
243 | b = Balloc (ptr, k); | |
244 | b->_x[0] = y9; | |
245 | b->_wds = 1; | |
246 | #else | |
247 | b = Balloc (ptr, k + 1); | |
248 | b->_x[0] = y9 & 0xffff; | |
249 | b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1; | |
250 | #endif | |
251 | ||
252 | i = 9; | |
253 | if (9 < nd0) | |
254 | { | |
255 | s += 9; | |
256 | do | |
257 | b = multadd (ptr, b, 10, *s++ - '0'); | |
258 | while (++i < nd0); | |
259 | s++; | |
260 | } | |
261 | else | |
262 | s += 10; | |
263 | for (; i < nd; i++) | |
264 | b = multadd (ptr, b, 10, *s++ - '0'); | |
265 | return b; | |
266 | } | |
267 | ||
268 | int | |
269 | _DEFUN (hi0bits, | |
80f798b4 | 270 | (x), register __ULong x) |
f911ba98 TT |
271 | { |
272 | register int k = 0; | |
273 | ||
274 | if (!(x & 0xffff0000)) | |
275 | { | |
276 | k = 16; | |
277 | x <<= 16; | |
278 | } | |
279 | if (!(x & 0xff000000)) | |
280 | { | |
281 | k += 8; | |
282 | x <<= 8; | |
283 | } | |
284 | if (!(x & 0xf0000000)) | |
285 | { | |
286 | k += 4; | |
287 | x <<= 4; | |
288 | } | |
289 | if (!(x & 0xc0000000)) | |
290 | { | |
291 | k += 2; | |
292 | x <<= 2; | |
293 | } | |
294 | if (!(x & 0x80000000)) | |
295 | { | |
296 | k++; | |
297 | if (!(x & 0x40000000)) | |
298 | return 32; | |
299 | } | |
300 | return k; | |
301 | } | |
302 | ||
303 | int | |
80f798b4 | 304 | _DEFUN (lo0bits, (y), __ULong *y) |
f911ba98 TT |
305 | { |
306 | register int k; | |
80f798b4 | 307 | register __ULong x = *y; |
f911ba98 TT |
308 | |
309 | if (x & 7) | |
310 | { | |
311 | if (x & 1) | |
312 | return 0; | |
313 | if (x & 2) | |
314 | { | |
315 | *y = x >> 1; | |
316 | return 1; | |
317 | } | |
318 | *y = x >> 2; | |
319 | return 2; | |
320 | } | |
321 | k = 0; | |
322 | if (!(x & 0xffff)) | |
323 | { | |
324 | k = 16; | |
325 | x >>= 16; | |
326 | } | |
327 | if (!(x & 0xff)) | |
328 | { | |
329 | k += 8; | |
330 | x >>= 8; | |
331 | } | |
332 | if (!(x & 0xf)) | |
333 | { | |
334 | k += 4; | |
335 | x >>= 4; | |
336 | } | |
337 | if (!(x & 0x3)) | |
338 | { | |
339 | k += 2; | |
340 | x >>= 2; | |
341 | } | |
342 | if (!(x & 1)) | |
343 | { | |
344 | k++; | |
345 | x >>= 1; | |
80f798b4 | 346 | if (!x & 1) |
f911ba98 TT |
347 | return 32; |
348 | } | |
349 | *y = x; | |
350 | return k; | |
351 | } | |
352 | ||
80f798b4 TT |
353 | _Bigint * |
354 | _DEFUN (i2b, (ptr, i), struct _reent * ptr _AND int i) | |
f911ba98 | 355 | { |
80f798b4 | 356 | _Bigint *b; |
f911ba98 TT |
357 | |
358 | b = Balloc (ptr, 1); | |
359 | b->_x[0] = i; | |
360 | b->_wds = 1; | |
361 | return b; | |
362 | } | |
363 | ||
80f798b4 TT |
364 | _Bigint * |
365 | _DEFUN (mult, (ptr, a, b), struct _reent * ptr _AND _Bigint * a _AND _Bigint * b) | |
f911ba98 | 366 | { |
80f798b4 | 367 | _Bigint *c; |
f911ba98 | 368 | int k, wa, wb, wc; |
80f798b4 TT |
369 | __ULong carry, y, z; |
370 | __ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; | |
f911ba98 | 371 | #ifdef Pack_32 |
80f798b4 | 372 | __ULong z2; |
f911ba98 TT |
373 | #endif |
374 | ||
375 | if (a->_wds < b->_wds) | |
376 | { | |
377 | c = a; | |
378 | a = b; | |
379 | b = c; | |
380 | } | |
381 | k = a->_k; | |
382 | wa = a->_wds; | |
383 | wb = b->_wds; | |
384 | wc = wa + wb; | |
385 | if (wc > a->_maxwds) | |
386 | k++; | |
387 | c = Balloc (ptr, k); | |
388 | for (x = c->_x, xa = x + wc; x < xa; x++) | |
389 | *x = 0; | |
390 | xa = a->_x; | |
391 | xae = xa + wa; | |
392 | xb = b->_x; | |
393 | xbe = xb + wb; | |
394 | xc0 = c->_x; | |
395 | #ifdef Pack_32 | |
396 | for (; xb < xbe; xb++, xc0++) | |
397 | { | |
80f798b4 | 398 | if ((y = *xb & 0xffff) != 0) |
f911ba98 TT |
399 | { |
400 | x = xa; | |
401 | xc = xc0; | |
402 | carry = 0; | |
403 | do | |
404 | { | |
405 | z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; | |
406 | carry = z >> 16; | |
407 | z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; | |
408 | carry = z2 >> 16; | |
409 | Storeinc (xc, z2, z); | |
410 | } | |
411 | while (x < xae); | |
412 | *xc = carry; | |
413 | } | |
80f798b4 | 414 | if ((y = *xb >> 16) != 0) |
f911ba98 TT |
415 | { |
416 | x = xa; | |
417 | xc = xc0; | |
418 | carry = 0; | |
419 | z2 = *xc; | |
420 | do | |
421 | { | |
422 | z = (*x & 0xffff) * y + (*xc >> 16) + carry; | |
423 | carry = z >> 16; | |
424 | Storeinc (xc, z, z2); | |
425 | z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; | |
426 | carry = z2 >> 16; | |
427 | } | |
428 | while (x < xae); | |
429 | *xc = z2; | |
430 | } | |
431 | } | |
432 | #else | |
433 | for (; xb < xbe; xc0++) | |
434 | { | |
80f798b4 | 435 | if (y = *xb++) |
f911ba98 TT |
436 | { |
437 | x = xa; | |
438 | xc = xc0; | |
439 | carry = 0; | |
440 | do | |
441 | { | |
442 | z = *x++ * y + *xc + carry; | |
443 | carry = z >> 16; | |
444 | *xc++ = z & 0xffff; | |
445 | } | |
446 | while (x < xae); | |
447 | *xc = carry; | |
448 | } | |
449 | } | |
450 | #endif | |
451 | for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc); | |
452 | c->_wds = wc; | |
453 | return c; | |
454 | } | |
455 | ||
80f798b4 | 456 | _Bigint * |
f911ba98 | 457 | _DEFUN (pow5mult, |
80f798b4 | 458 | (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k) |
f911ba98 | 459 | { |
80f798b4 | 460 | _Bigint *b1, *p5, *p51; |
f911ba98 TT |
461 | int i; |
462 | static _CONST int p05[3] = {5, 25, 125}; | |
463 | ||
80f798b4 | 464 | if ((i = k & 3) != 0) |
f911ba98 TT |
465 | b = multadd (ptr, b, p05[i - 1], 0); |
466 | ||
467 | if (!(k >>= 2)) | |
468 | return b; | |
80f798b4 TT |
469 | _REENT_CHECK_MP(ptr); |
470 | if (!(p5 = _REENT_MP_P5S(ptr))) | |
f911ba98 TT |
471 | { |
472 | /* first time */ | |
80f798b4 | 473 | p5 = _REENT_MP_P5S(ptr) = i2b (ptr, 625); |
f911ba98 TT |
474 | p5->_next = 0; |
475 | } | |
476 | for (;;) | |
477 | { | |
478 | if (k & 1) | |
479 | { | |
480 | b1 = mult (ptr, b, p5); | |
481 | Bfree (ptr, b); | |
482 | b = b1; | |
483 | } | |
484 | if (!(k >>= 1)) | |
485 | break; | |
486 | if (!(p51 = p5->_next)) | |
487 | { | |
488 | p51 = p5->_next = mult (ptr, p5, p5); | |
489 | p51->_next = 0; | |
490 | } | |
491 | p5 = p51; | |
492 | } | |
493 | return b; | |
494 | } | |
495 | ||
80f798b4 TT |
496 | _Bigint * |
497 | _DEFUN (lshift, (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k) | |
f911ba98 TT |
498 | { |
499 | int i, k1, n, n1; | |
80f798b4 TT |
500 | _Bigint *b1; |
501 | __ULong *x, *x1, *xe, z; | |
f911ba98 TT |
502 | |
503 | #ifdef Pack_32 | |
504 | n = k >> 5; | |
505 | #else | |
506 | n = k >> 4; | |
507 | #endif | |
508 | k1 = b->_k; | |
509 | n1 = n + b->_wds + 1; | |
510 | for (i = b->_maxwds; n1 > i; i <<= 1) | |
511 | k1++; | |
512 | b1 = Balloc (ptr, k1); | |
513 | x1 = b1->_x; | |
514 | for (i = 0; i < n; i++) | |
515 | *x1++ = 0; | |
516 | x = b->_x; | |
517 | xe = x + b->_wds; | |
518 | #ifdef Pack_32 | |
519 | if (k &= 0x1f) | |
520 | { | |
521 | k1 = 32 - k; | |
522 | z = 0; | |
523 | do | |
524 | { | |
525 | *x1++ = *x << k | z; | |
526 | z = *x++ >> k1; | |
527 | } | |
528 | while (x < xe); | |
80f798b4 | 529 | if ((*x1 = z) != 0) |
f911ba98 TT |
530 | ++n1; |
531 | } | |
532 | #else | |
533 | if (k &= 0xf) | |
534 | { | |
535 | k1 = 16 - k; | |
536 | z = 0; | |
537 | do | |
538 | { | |
80f798b4 | 539 | *x1++ = *x << k & 0xffff | z; |
f911ba98 TT |
540 | z = *x++ >> k1; |
541 | } | |
542 | while (x < xe); | |
80f798b4 | 543 | if (*x1 = z) |
f911ba98 TT |
544 | ++n1; |
545 | } | |
546 | #endif | |
547 | else | |
548 | do | |
549 | *x1++ = *x++; | |
550 | while (x < xe); | |
551 | b1->_wds = n1 - 1; | |
552 | Bfree (ptr, b); | |
553 | return b1; | |
554 | } | |
555 | ||
556 | int | |
80f798b4 | 557 | _DEFUN (cmp, (a, b), _Bigint * a _AND _Bigint * b) |
f911ba98 | 558 | { |
80f798b4 | 559 | __ULong *xa, *xa0, *xb, *xb0; |
f911ba98 TT |
560 | int i, j; |
561 | ||
562 | i = a->_wds; | |
563 | j = b->_wds; | |
564 | #ifdef DEBUG | |
565 | if (i > 1 && !a->_x[i - 1]) | |
566 | Bug ("cmp called with a->_x[a->_wds-1] == 0"); | |
567 | if (j > 1 && !b->_x[j - 1]) | |
568 | Bug ("cmp called with b->_x[b->_wds-1] == 0"); | |
569 | #endif | |
570 | if (i -= j) | |
571 | return i; | |
572 | xa0 = a->_x; | |
573 | xa = xa0 + j; | |
574 | xb0 = b->_x; | |
575 | xb = xb0 + j; | |
576 | for (;;) | |
577 | { | |
578 | if (*--xa != *--xb) | |
579 | return *xa < *xb ? -1 : 1; | |
580 | if (xa <= xa0) | |
581 | break; | |
582 | } | |
583 | return 0; | |
584 | } | |
585 | ||
80f798b4 TT |
586 | _Bigint * |
587 | _DEFUN (diff, (ptr, a, b), struct _reent * ptr _AND | |
588 | _Bigint * a _AND _Bigint * b) | |
f911ba98 | 589 | { |
80f798b4 | 590 | _Bigint *c; |
f911ba98 | 591 | int i, wa, wb; |
80f798b4 TT |
592 | __Long borrow, y; /* We need signed shifts here. */ |
593 | __ULong *xa, *xae, *xb, *xbe, *xc; | |
f911ba98 | 594 | #ifdef Pack_32 |
80f798b4 | 595 | __Long z; |
f911ba98 TT |
596 | #endif |
597 | ||
598 | i = cmp (a, b); | |
599 | if (!i) | |
600 | { | |
601 | c = Balloc (ptr, 0); | |
602 | c->_wds = 1; | |
603 | c->_x[0] = 0; | |
604 | return c; | |
605 | } | |
606 | if (i < 0) | |
607 | { | |
608 | c = a; | |
609 | a = b; | |
610 | b = c; | |
611 | i = 1; | |
612 | } | |
613 | else | |
614 | i = 0; | |
615 | c = Balloc (ptr, a->_k); | |
616 | c->_sign = i; | |
617 | wa = a->_wds; | |
618 | xa = a->_x; | |
619 | xae = xa + wa; | |
620 | wb = b->_wds; | |
621 | xb = b->_x; | |
622 | xbe = xb + wb; | |
623 | xc = c->_x; | |
624 | borrow = 0; | |
625 | #ifdef Pack_32 | |
626 | do | |
627 | { | |
628 | y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; | |
629 | borrow = y >> 16; | |
630 | Sign_Extend (borrow, y); | |
631 | z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; | |
632 | borrow = z >> 16; | |
633 | Sign_Extend (borrow, z); | |
634 | Storeinc (xc, z, y); | |
635 | } | |
636 | while (xb < xbe); | |
637 | while (xa < xae) | |
638 | { | |
639 | y = (*xa & 0xffff) + borrow; | |
640 | borrow = y >> 16; | |
641 | Sign_Extend (borrow, y); | |
642 | z = (*xa++ >> 16) + borrow; | |
643 | borrow = z >> 16; | |
644 | Sign_Extend (borrow, z); | |
645 | Storeinc (xc, z, y); | |
646 | } | |
647 | #else | |
648 | do | |
649 | { | |
650 | y = *xa++ - *xb++ + borrow; | |
651 | borrow = y >> 16; | |
652 | Sign_Extend (borrow, y); | |
653 | *xc++ = y & 0xffff; | |
654 | } | |
655 | while (xb < xbe); | |
656 | while (xa < xae) | |
657 | { | |
658 | y = *xa++ + borrow; | |
659 | borrow = y >> 16; | |
660 | Sign_Extend (borrow, y); | |
661 | *xc++ = y & 0xffff; | |
662 | } | |
663 | #endif | |
664 | while (!*--xc) | |
665 | wa--; | |
666 | c->_wds = wa; | |
667 | return c; | |
668 | } | |
669 | ||
670 | double | |
671 | _DEFUN (ulp, (_x), double _x) | |
672 | { | |
673 | union double_union x, a; | |
902f7d1a | 674 | register int32_t L; |
f911ba98 TT |
675 | |
676 | x.d = _x; | |
677 | ||
678 | L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1; | |
679 | #ifndef Sudden_Underflow | |
680 | if (L > 0) | |
681 | { | |
682 | #endif | |
683 | #ifdef IBM | |
684 | L |= Exp_msk1 >> 4; | |
685 | #endif | |
686 | word0 (a) = L; | |
687 | #ifndef _DOUBLE_IS_32BITS | |
688 | word1 (a) = 0; | |
689 | #endif | |
690 | ||
691 | #ifndef Sudden_Underflow | |
692 | } | |
693 | else | |
694 | { | |
695 | L = -L >> Exp_shift; | |
696 | if (L < Exp_shift) | |
697 | { | |
698 | word0 (a) = 0x80000 >> L; | |
699 | #ifndef _DOUBLE_IS_32BITS | |
700 | word1 (a) = 0; | |
701 | #endif | |
702 | } | |
703 | else | |
704 | { | |
705 | word0 (a) = 0; | |
706 | L -= Exp_shift; | |
707 | #ifndef _DOUBLE_IS_32BITS | |
80f798b4 | 708 | word1 (a) = L >= 31 ? 1 : 1 << (31 - L); |
f911ba98 TT |
709 | #endif |
710 | } | |
711 | } | |
712 | #endif | |
713 | return a.d; | |
714 | } | |
715 | ||
716 | double | |
717 | _DEFUN (b2d, (a, e), | |
80f798b4 | 718 | _Bigint * a _AND int *e) |
f911ba98 | 719 | { |
80f798b4 | 720 | __ULong *xa, *xa0, w, y, z; |
f911ba98 TT |
721 | int k; |
722 | union double_union d; | |
723 | #ifdef VAX | |
80f798b4 | 724 | __ULong d0, d1; |
f911ba98 TT |
725 | #else |
726 | #define d0 word0(d) | |
727 | #define d1 word1(d) | |
728 | #endif | |
729 | ||
730 | xa0 = a->_x; | |
731 | xa = xa0 + a->_wds; | |
732 | y = *--xa; | |
733 | #ifdef DEBUG | |
734 | if (!y) | |
735 | Bug ("zero y in b2d"); | |
736 | #endif | |
737 | k = hi0bits (y); | |
738 | *e = 32 - k; | |
739 | #ifdef Pack_32 | |
740 | if (k < Ebits) | |
741 | { | |
742 | d0 = Exp_1 | y >> (Ebits - k); | |
743 | w = xa > xa0 ? *--xa : 0; | |
744 | #ifndef _DOUBLE_IS_32BITS | |
80f798b4 | 745 | d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k); |
f911ba98 TT |
746 | #endif |
747 | goto ret_d; | |
748 | } | |
749 | z = xa > xa0 ? *--xa : 0; | |
750 | if (k -= Ebits) | |
751 | { | |
752 | d0 = Exp_1 | y << k | z >> (32 - k); | |
753 | y = xa > xa0 ? *--xa : 0; | |
754 | #ifndef _DOUBLE_IS_32BITS | |
755 | d1 = z << k | y >> (32 - k); | |
756 | #endif | |
757 | } | |
758 | else | |
759 | { | |
760 | d0 = Exp_1 | y; | |
761 | #ifndef _DOUBLE_IS_32BITS | |
762 | d1 = z; | |
763 | #endif | |
764 | } | |
765 | #else | |
766 | if (k < Ebits + 16) | |
767 | { | |
768 | z = xa > xa0 ? *--xa : 0; | |
80f798b4 | 769 | d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; |
f911ba98 TT |
770 | w = xa > xa0 ? *--xa : 0; |
771 | y = xa > xa0 ? *--xa : 0; | |
80f798b4 | 772 | d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; |
f911ba98 TT |
773 | goto ret_d; |
774 | } | |
775 | z = xa > xa0 ? *--xa : 0; | |
776 | w = xa > xa0 ? *--xa : 0; | |
777 | k -= Ebits + 16; | |
80f798b4 | 778 | d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; |
f911ba98 | 779 | y = xa > xa0 ? *--xa : 0; |
80f798b4 | 780 | d1 = w << k + 16 | y << k; |
f911ba98 TT |
781 | #endif |
782 | ret_d: | |
783 | #ifdef VAX | |
784 | word0 (d) = d0 >> 16 | d0 << 16; | |
785 | word1 (d) = d1 >> 16 | d1 << 16; | |
786 | #else | |
787 | #undef d0 | |
788 | #undef d1 | |
789 | #endif | |
790 | return d.d; | |
791 | } | |
792 | ||
80f798b4 | 793 | _Bigint * |
f911ba98 TT |
794 | _DEFUN (d2b, |
795 | (ptr, _d, e, bits), | |
80f798b4 | 796 | struct _reent * ptr _AND |
f911ba98 TT |
797 | double _d _AND |
798 | int *e _AND | |
799 | int *bits) | |
800 | ||
801 | { | |
802 | union double_union d; | |
80f798b4 | 803 | _Bigint *b; |
f911ba98 | 804 | int de, i, k; |
80f798b4 | 805 | __ULong *x, y, z; |
f911ba98 | 806 | #ifdef VAX |
80f798b4 TT |
807 | __ULong d0, d1; |
808 | #endif | |
f911ba98 | 809 | d.d = _d; |
80f798b4 | 810 | #ifdef VAX |
f911ba98 TT |
811 | d0 = word0 (d) >> 16 | word0 (d) << 16; |
812 | d1 = word1 (d) >> 16 | word1 (d) << 16; | |
813 | #else | |
814 | #define d0 word0(d) | |
815 | #define d1 word1(d) | |
816 | d.d = _d; | |
817 | #endif | |
818 | ||
819 | #ifdef Pack_32 | |
820 | b = Balloc (ptr, 1); | |
821 | #else | |
822 | b = Balloc (ptr, 2); | |
823 | #endif | |
824 | x = b->_x; | |
825 | ||
826 | z = d0 & Frac_mask; | |
827 | d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ | |
828 | #ifdef Sudden_Underflow | |
829 | de = (int) (d0 >> Exp_shift); | |
830 | #ifndef IBM | |
831 | z |= Exp_msk11; | |
832 | #endif | |
833 | #else | |
80f798b4 | 834 | if ((de = (int) (d0 >> Exp_shift)) != 0) |
f911ba98 TT |
835 | z |= Exp_msk1; |
836 | #endif | |
837 | #ifdef Pack_32 | |
838 | #ifndef _DOUBLE_IS_32BITS | |
80f798b4 | 839 | if (d1) |
f911ba98 | 840 | { |
80f798b4 TT |
841 | y = d1; |
842 | k = lo0bits (&y); | |
843 | if (k) | |
f911ba98 | 844 | { |
80f798b4 | 845 | x[0] = y | z << (32 - k); |
f911ba98 TT |
846 | z >>= k; |
847 | } | |
848 | else | |
849 | x[0] = y; | |
850 | i = b->_wds = (x[1] = z) ? 2 : 1; | |
851 | } | |
852 | else | |
853 | #endif | |
854 | { | |
855 | #ifdef DEBUG | |
856 | if (!z) | |
857 | Bug ("Zero passed to d2b"); | |
858 | #endif | |
859 | k = lo0bits (&z); | |
860 | x[0] = z; | |
861 | i = b->_wds = 1; | |
862 | #ifndef _DOUBLE_IS_32BITS | |
863 | k += 32; | |
864 | #endif | |
865 | } | |
866 | #else | |
80f798b4 | 867 | if (d1) |
f911ba98 | 868 | { |
80f798b4 TT |
869 | y = d1; |
870 | k = lo0bits (&y); | |
871 | if (k) | |
f911ba98 TT |
872 | if (k >= 16) |
873 | { | |
80f798b4 TT |
874 | x[0] = y | z << 32 - k & 0xffff; |
875 | x[1] = z >> k - 16 & 0xffff; | |
f911ba98 TT |
876 | x[2] = z >> k; |
877 | i = 2; | |
878 | } | |
879 | else | |
880 | { | |
881 | x[0] = y & 0xffff; | |
80f798b4 | 882 | x[1] = y >> 16 | z << 16 - k & 0xffff; |
f911ba98 | 883 | x[2] = z >> k & 0xffff; |
80f798b4 | 884 | x[3] = z >> k + 16; |
f911ba98 TT |
885 | i = 3; |
886 | } | |
887 | else | |
888 | { | |
889 | x[0] = y & 0xffff; | |
890 | x[1] = y >> 16; | |
891 | x[2] = z & 0xffff; | |
892 | x[3] = z >> 16; | |
893 | i = 3; | |
894 | } | |
895 | } | |
896 | else | |
897 | { | |
898 | #ifdef DEBUG | |
899 | if (!z) | |
900 | Bug ("Zero passed to d2b"); | |
901 | #endif | |
902 | k = lo0bits (&z); | |
903 | if (k >= 16) | |
904 | { | |
905 | x[0] = z; | |
906 | i = 0; | |
907 | } | |
908 | else | |
909 | { | |
910 | x[0] = z & 0xffff; | |
911 | x[1] = z >> 16; | |
912 | i = 1; | |
913 | } | |
914 | k += 32; | |
915 | } | |
916 | while (!x[i]) | |
917 | --i; | |
918 | b->_wds = i + 1; | |
919 | #endif | |
920 | #ifndef Sudden_Underflow | |
921 | if (de) | |
922 | { | |
923 | #endif | |
924 | #ifdef IBM | |
925 | *e = (de - Bias - (P - 1) << 2) + k; | |
926 | *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask); | |
927 | #else | |
928 | *e = de - Bias - (P - 1) + k; | |
929 | *bits = P - k; | |
930 | #endif | |
931 | #ifndef Sudden_Underflow | |
932 | } | |
933 | else | |
934 | { | |
935 | *e = de - Bias - (P - 1) + 1 + k; | |
936 | #ifdef Pack_32 | |
937 | *bits = 32 * i - hi0bits (x[i - 1]); | |
938 | #else | |
939 | *bits = (i + 2) * 16 - hi0bits (x[i]); | |
940 | #endif | |
941 | } | |
942 | #endif | |
943 | return b; | |
944 | } | |
945 | #undef d0 | |
946 | #undef d1 | |
947 | ||
948 | double | |
80f798b4 | 949 | _DEFUN (ratio, (a, b), _Bigint * a _AND _Bigint * b) |
f911ba98 TT |
950 | |
951 | { | |
952 | union double_union da, db; | |
953 | int k, ka, kb; | |
954 | ||
955 | da.d = b2d (a, &ka); | |
956 | db.d = b2d (b, &kb); | |
957 | #ifdef Pack_32 | |
958 | k = ka - kb + 32 * (a->_wds - b->_wds); | |
959 | #else | |
960 | k = ka - kb + 16 * (a->_wds - b->_wds); | |
961 | #endif | |
962 | #ifdef IBM | |
963 | if (k > 0) | |
964 | { | |
965 | word0 (da) += (k >> 2) * Exp_msk1; | |
966 | if (k &= 3) | |
967 | da.d *= 1 << k; | |
968 | } | |
969 | else | |
970 | { | |
971 | k = -k; | |
972 | word0 (db) += (k >> 2) * Exp_msk1; | |
973 | if (k &= 3) | |
974 | db.d *= 1 << k; | |
975 | } | |
976 | #else | |
977 | if (k > 0) | |
978 | word0 (da) += k * Exp_msk1; | |
979 | else | |
980 | { | |
981 | k = -k; | |
982 | word0 (db) += k * Exp_msk1; | |
983 | } | |
984 | #endif | |
985 | return da.d / db.d; | |
986 | } | |
987 | ||
988 | ||
989 | _CONST double | |
990 | tens[] = | |
991 | { | |
992 | 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, | |
993 | 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, | |
994 | 1e20, 1e21, 1e22, 1e23, 1e24 | |
995 | ||
996 | }; | |
997 | ||
998 | #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800) | |
999 | _CONST double bigtens[] = | |
1000 | {1e16, 1e32, 1e64, 1e128, 1e256}; | |
1001 | ||
1002 | _CONST double tinytens[] = | |
1003 | {1e-16, 1e-32, 1e-64, 1e-128, 1e-256}; | |
1004 | #else | |
1005 | _CONST double bigtens[] = | |
1006 | {1e16, 1e32}; | |
1007 | ||
1008 | _CONST double tinytens[] = | |
1009 | {1e-16, 1e-32}; | |
1010 | #endif | |
1011 | ||
1012 | ||
80f798b4 TT |
1013 | double |
1014 | _DEFUN (_mprec_log10, (dig), | |
1015 | int dig) | |
1016 | { | |
1017 | double v = 1.0; | |
1018 | if (dig < 24) | |
1019 | return tens[dig]; | |
1020 | while (dig > 0) | |
1021 | { | |
1022 | v *= 10; | |
1023 | dig--; | |
1024 | } | |
1025 | return v; | |
1026 | } |