source: trunk/gcc/libjava/java/lang/mprec.c

Last change on this file was 2, checked in by bird, 22 years ago

Initial revision

  • Property cvs2svn:cvs-rev set to 1.1
  • Property svn:eol-style set to native
  • Property svn:executable set to *
File size: 17.5 KB
Line 
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
120void
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
221int
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
256int
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
508int
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
622double
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
668double
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
734ret_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
894double
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
Note: See TracBrowser for help on using the repository browser.