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