1 // $Id: double.cpp,v 1.35 2004/06/02 11:26:21 elliott-oss Exp $
2 //
3 // This software is subject to the terms of the IBM Jikes Compiler
4 // License Agreement available at the following URL:
5 // http://ibm.com/developerworks/opensource/jikes.
6 // Copyright (C) 1996, 2004 IBM Corporation and others. All Rights Reserved.
7 // You must accept the terms of that agreement to use this software.
8 //
9 //
10 // NOTES: The IEEE 754 emulation code in double.h and double.cpp within
11 // Jikes are adapted from code written by Alan M. Webb of IBM's Hursley
12 // lab in porting the Sun JDK to System/390.
13 //
14 //
15 //
16 // In addition, the code for emulating the remainder operator, %, is
17 // adapted from e_fmod.c, part of fdlibm, the Freely Distributable Math
18 // Library mentioned in the documentation of java.lang.StrictMath. The
19 // original library is available at http://netlib2.cs.utk.edu/fdlibm.
20 //
21 // The code from fdlibm is copyrighted, as follows:
22 // ====================================================
23 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
24 //
25 // Developed at SunSoft, a Sun Microsystems, Inc. business.
26 // Permission to use, copy, modify, and distribute this
27 // software is freely granted, provided that this notice
28 // is preserved.
29 // ====================================================
30 //
31 //
32 //
33 // Likewise, the code for accurate conversions between floating point
34 // and decimal strings, in double.h, double.cpp, platform.h, and
35 // platform.cpp, is adapted from dtoa.c. The original code can be
36 // found at http://netlib2.cs.utk.edu/fp/dtoa.c.
37 //
38 // The code in dtoa.c is copyrighted as follows:
39 //****************************************************************
40 //*
41 //* The author of this software is David M. Gay.
42 //*
43 //* Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
44 //*
45 //* Permission to use, copy, modify, and distribute this software for any
46 //* purpose without fee is hereby granted, provided that this entire notice
47 //* is included in all copies of any software which is or includes a copy
48 //* or modification of this software and in all copies of the supporting
49 //* documentation for such software.
50 //*
51 //* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
52 //* WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
53 //* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
54 //* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
55 //*
56 //***************************************************************/
57 //
58 //
59
60 #include "double.h"
61 #include "long.h"
62 #include "code.h"
63
64 #ifdef HAVE_JIKES_NAMESPACE
65 namespace Jikes { // Open namespace Jikes block
66 #endif
67
68 #ifndef HAVE_MEMBER_CONSTANTS
69 // VC++ can't cope with constant class members
70 IEEEfloat IEEEfloat::tens[] = {
71 1e0f, 1e1f, 1e2f, 1e3f, 1e4f, 1e5f, 1e6f,
72 1e7f, 1e8f, 1e9f, 1e10f
73 };
74 IEEEfloat IEEEfloat::bigtens[] = {
75 1e8f, 1e16f, 1e32f
76 };
77
78 IEEEdouble IEEEdouble::tens[] = {
79 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7,
80 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15,
81 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22
82 };
83 IEEEdouble IEEEdouble::bigtens[] = {
84 1e16, 1e32, 1e64, 1e128, 1e256
85 };
86 #else
87 const IEEEfloat IEEEfloat::tens[] = {
88 1e0f, 1e1f, 1e2f, 1e3f, 1e4f, 1e5f, 1e6f,
89 1e7f, 1e8f, 1e9f, 1e10f
90 };
91 const IEEEfloat IEEEfloat::bigtens[] = {
92 1e8f, 1e16f, 1e32f
93 };
94
95 const IEEEdouble IEEEdouble::tens[] = {
96 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7,
97 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15,
98 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22
99 };
100 const IEEEdouble IEEEdouble::bigtens[] = {
101 1e16, 1e32, 1e64, 1e128, 1e256
102 };
103 #endif
104
105
IEEEfloat(i4 a)106 IEEEfloat::IEEEfloat(i4 a)
107 {
108 #ifdef HAVE_IEEE754
109 value.float_value = (float) a;
110 #else
111 int sign = 0;
112
113 if (a < 0)
114 {
115 a = -a; // even works for MIN_INT!
116 sign = 1;
117 }
118 *this = (a == 0) ? POSITIVE_ZERO()
119 : Normalize(sign, FRACT_SIZE, (u4) a);
120 #endif // HAVE_IEEE754
121 }
122
IEEEfloat(const LongInt & a)123 IEEEfloat::IEEEfloat(const LongInt &a)
124 {
125 #ifdef HAVE_IEEE754
126 # ifdef HAVE_64BIT_TYPES
127 value.float_value = (float)(i8) a.Words();
128 # else
129 value.float_value = ((float)(i4) a.HighWord() * (float) 0x40000000 * 4.0f)
130 + (float) a.LowWord();
131 # endif // HAVE_64BIT_TYPES
132 #else
133 //
134 // Unfortunately, we cannot recycle the LongInt.DoubleValue() method, since
135 // in rare cases the double rounding puts us off by one bit.
136 //
137 int sign = 0;
138 LongInt l = a;
139
140 if (a < 0)
141 {
142 l = -a;
143 sign = 1;
144 if (l < 0) // special case MIN_LONG
145 {
146 value.word = MIN_LONG_F;
147 return;
148 }
149 }
150 if (l == 0)
151 *this = POSITIVE_ZERO();
152 else if (l.HighWord() == 0)
153 *this = Normalize(sign, FRACT_SIZE, l.LowWord());
154 else
155 {
156 int exponent = FRACT_SIZE - 1, sticky = 0;
157 while (l.HighWord())
158 {
159 sticky |= (l.LowWord() & BYTE_MASK) ? 1 : 0;
160 l >>= 8;
161 exponent += 8;
162 }
163 u4 low = l.LowWord();
164 if ((i4) low < 0)
165 {
166 sticky |= (low & BYTE_MASK) ? 1 : 0;
167 low >>= 8;
168 exponent += 8;
169 }
170 *this = Normalize(sign, exponent, low + low + sticky);
171 }
172 #endif // HAVE_IEEE754
173 }
174
IEEEfloat(const IEEEdouble & d)175 IEEEfloat::IEEEfloat(const IEEEdouble &d)
176 {
177 #ifdef HAVE_IEEE754
178 value.float_value = (float) d.DoubleView();
179 #else
180 // Either true zero, denormalized, or too small
181 if (d.Exponent() < -BIAS - 30)
182 *this = (d.IsPositive() ? POSITIVE_ZERO() : NEGATIVE_ZERO());
183 else
184 {
185 if (d.IsPositiveInfinity())
186 *this = POSITIVE_INFINITY();
187 else if (d.IsNegativeInfinity())
188 *this = NEGATIVE_INFINITY();
189 else if (d.IsNaN())
190 *this = NaN();
191 else
192 {
193 //
194 // A regular, normalized number - do work on the parts
195 // Shift to 26th position, add implicit msb, rounding bits
196 //
197 LongInt fract = d.Fraction() << 5;
198 u4 fraction = fract.HighWord() | (fract.LowWord() ? 0x02000001
199 : 0x02000000);
200
201 *this = Normalize(d.Sign(), d.Exponent() - 2, fraction);
202 }
203 }
204 #endif // HAVE_IEEE754
205 }
206
Adjust(const BigInt & delta,const BigInt & bs,const bool dsign)207 bool IEEEfloat::Adjust(const BigInt &delta, const BigInt &bs, const bool dsign)
208 {
209 IEEEfloat aadj, aadj1;
210 i4 y;
211
212 aadj = Ratio(delta, bs);
213 if (aadj <= 2)
214 {
215 if (dsign)
216 aadj = aadj1 = 1;
217 else if (FractBits())
218 {
219 if (value.word == 1)
220 {
221 // underflow
222 *this = POSITIVE_ZERO();
223 return true;
224 }
225 aadj = 1;
226 aadj1 = -1;
227 }
228 else
229 {
230 //
231 // special case - mantissa is power of 2
232 //
233 if (aadj < 1.0f)
234 aadj = 0.5f;
235 else
236 aadj *= 0.5f;
237 aadj1 = -aadj;
238 }
239 }
240 else
241 {
242 aadj *= 0.5f;
243 aadj1 = dsign ? aadj : -aadj;
244 }
245 y = Exponent();
246 //
247 // Check for overflow
248 //
249 if (y == BIAS)
250 {
251 IEEEfloat tmp(*this);
252 value.word -= (FRACT_SIZE + 1) * MIN_FRACT;
253 *this += aadj1 * Ulp();
254 if (Exponent() >= BIAS - FRACT_SIZE)
255 {
256 if (tmp.value.word == POS_INF - 1)
257 {
258 // overflow
259 *this = POSITIVE_INFINITY();
260 return true;
261 }
262 value.word = POS_INF - 1;
263 return false;
264 }
265 else
266 value.word += (FRACT_SIZE + 1) * MIN_FRACT;
267 }
268 else
269 {
270 //
271 // Compute adj so that the IEEE rounding rules will
272 // correctly round *this + adj in some half-way cases.
273 // If *this * Ulp() is denormalized, we must adjust aadj
274 // to avoid trouble from bits lost to denormalization.
275 //
276 if (y <= FRACT_SIZE - BIAS || aadj > 1)
277 {
278 aadj1 = IEEEfloat((aadj + 0.5f).IntValue());
279 if (! dsign)
280 aadj1 = -aadj1;
281 }
282 *this += aadj1 * Ulp();
283 }
284 if (y == Exponent())
285 {
286 //
287 // Can we stop now?
288 // The tolerances below are conservative.
289 //
290 aadj -= aadj.IntValue();
291 if (dsign || FractBits())
292 {
293 if (aadj < .4999999f || aadj > .5000001f)
294 return true;
295 }
296 else if (aadj < .4999999f / 2)
297 return true;
298 }
299 return false;
300 }
301
IEEEfloat(const char * str,bool check_invalid)302 IEEEfloat::IEEEfloat(const char *str, bool check_invalid)
303 {
304 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2,
305 e, e1, i, j, k;
306 int nd, // number of digits in mantissa (except extra '0's)
307 nd0, // number of digits before '.' (except leading '0's)
308 nf, // number of digits after '.' (except trailing '0's)
309 nz; // number of consecutive '0' digits
310 bool nz0, // whether leading zero exists
311 roundup, // whether to round up long string
312 sign, // if the string represents a negative
313 dsign, // the sign of delta
314 esign; // the sign of the exponent
315 char c;
316 const char *s, *s0, *s1;
317 i4 L;
318 i4 y, z;
319
320 sign = nz0 = roundup = false;
321 nz = 0;
322
323 //
324 // consume whitespace
325 //
326 for (s = str; ; s++)
327 {
328 switch (*s)
329 {
330 case U_MINUS:
331 sign = true;
332 // fallthrough
333 case U_PLUS:
334 if (*++s)
335 goto break2;
336 // fallthrough
337 case U_NU:
338 *this = NaN();
339 return;
340 case U_SP:
341 case U_HT:
342 case U_FF:
343 case U_LF:
344 case U_CR:
345 continue;
346 default:
347 goto break2;
348 }
349 }
350 break2:
351 //
352 // Consume leading zeroes.
353 //
354 if (*s == U_0)
355 {
356 nz0 = true;
357 while (*++s == U_0);
358
359 if (! *s)
360 {
361 *this = sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
362 return;
363 }
364 //
365 // Parse hexadecimal floating point.
366 //
367 else if (*s == U_x || *s == U_X)
368 {
369 i4 fraction = 0;
370 // Exponent adjustment, based on '.' location.
371 int shift = FRACT_SIZE;
372 bool seen_dot = false;
373 while (*++s && *s == U_0); // Consume leading zeroes.
374 c = *s;
375 if (c == U_DOT)
376 {
377 seen_dot = true;
378 while (*++s && *s == U_0)
379 shift -= 4;
380 c = *s;
381 }
382 if (! c || c == U_p || c == U_P)
383 {
384 *this = sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
385 return;
386 }
387 // To avoid overflow, stop after enough bits have been read.
388 for (i = 0; i < (FRACT_SIZE >> 2) + 2 && c; i++, c = *++s)
389 {
390 if (c == U_DOT)
391 {
392 if (seen_dot)
393 break;
394 c = *++s;
395 seen_dot = true;
396 }
397 int value;
398 if (Code::IsHexDigit(c))
399 value = Code::Value(c);
400 else break;
401 if (seen_dot)
402 shift -= 4;
403 fraction = (fraction << 4) + value;
404 }
405 // Round any remaining bits.
406 bool sticky = false;
407 while (c == U_DOT || Code::IsHexDigit(c))
408 {
409 if (c == U_DOT)
410 {
411 if (seen_dot)
412 break;
413 seen_dot = true;
414 }
415 else
416 {
417 if (! seen_dot)
418 shift += 4;
419 if (c != U_0)
420 sticky = true;
421 }
422 c = *++s;
423 }
424 assert(fraction != 0);
425 if (sticky)
426 fraction++;
427 // On to the expononet.
428 int exponent = 0;
429 esign = false;
430 if (c == U_p || c == U_P)
431 {
432 if (*++s == U_MINUS)
433 {
434 esign = true;
435 s++;
436 }
437 else if (*s == U_PLUS)
438 s++;
439 while ((c = *s++))
440 {
441 if (! Code::IsDecimalDigit(c))
442 break;
443 exponent = exponent * 10 + c - U_0;
444 // Check for exponent overflow
445 if (exponent + shift > 19999)
446 {
447 if (check_invalid)
448 *this = NaN();
449 else
450 {
451 *this = esign ? POSITIVE_ZERO()
452 : POSITIVE_INFINITY();
453 if (sign)
454 *this = - *this;
455 }
456 return;
457 }
458 }
459 }
460 if (esign)
461 exponent = - exponent;
462 *this = Normalize(sign, exponent + shift, fraction);
463 if (check_invalid && (IsZero() || IsInfinite()))
464 *this = NaN();
465 return;
466 }
467 }
468
469 //
470 // parse before '.'
471 //
472 s0 = s;
473 y = z = 0;
474 for (nd = nf = 0; Code::IsDecimalDigit(c = *s); nd++, s++)
475 if (nd < 8)
476 y = 10 * y + c - U_0;
477 nd0 = nd;
478 if (c == U_DOT)
479 {
480 //
481 // parse after '.'
482 //
483 c = *++s;
484 if (!nd)
485 {
486 for ( ; c == U_0; c = *++s)
487 nz++;
488 if (c > U_0 && c <= U_9)
489 {
490 s0 = s;
491 nf += nz;
492 nz = 0;
493 }
494 }
495 for ( ; Code::IsDecimalDigit(c); c = *++s)
496 {
497 nz++;
498 if (c -= U_0)
499 {
500 nf += nz;
501 for (i = 1; i < nz; i++)
502 if (nd++ < 8)
503 y *= 10;
504 if (nd++ < 8)
505 y = 10 * y + c;
506 nz = 0;
507 }
508 }
509 }
510 //
511 // consume exponent
512 //
513 e = 0;
514 if (c == U_e || c == U_E) {
515 str = s;
516 esign = false;
517 switch (c = *++s)
518 {
519 case U_MINUS:
520 esign = true;
521 // fallthrough
522 case U_PLUS:
523 c = *++s;
524 }
525 if (Code::IsDecimalDigit(c))
526 {
527 if (!nd && !nz && !nz0)
528 {
529 *this = sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
530 return;
531 }
532 while (c == U_0)
533 c = *++s;
534 if (c > U_0 && c <= U_9)
535 {
536 L = c - U_0;
537 s1 = s;
538 while (Code::IsDecimalDigit(c = *++s))
539 L = 10 * L + c - U_0;
540 //
541 // Avoid confusion from exponents so large that e might
542 // overflow
543 if (s - s1 > 8 || L > 19999)
544 e = 19999;
545 else
546 e = L;
547 if (esign)
548 e = -e;
549 }
550 else
551 e = 0;
552 }
553 else
554 s = str;
555 }
556 if (!nd) {
557 *this = (!nz && !nz0) ? NaN()
558 : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
559 return;
560 }
561 //
562 // for long strings, round away all digits beyond maximum precise string
563 // for fraction, there are n decimal digits after '.' if lsb is 2^n;
564 // but the first m of these digits are '0', for d = x*10^m
565 // so, digits after MAX_DIGITS may be ignored
566 //
567 if (nd > MAX_DIGITS)
568 {
569 k = nd - MAX_DIGITS;
570 i = MAX_DIGITS - nd0;
571 if (i <= 0)
572 {
573 // decimal after last precise digit
574 nd0 = MAX_DIGITS;
575 nf = i;
576 j = 0;
577 }
578 else if (i == MAX_DIGITS)
579 {
580 // decimal before first precise digit
581 nf -= k;
582 j = 0;
583 }
584 else
585 {
586 // decimal inside precise digits
587 nf -= k;
588 j = 1;
589 }
590 roundup = s0[MAX_DIGITS - 1 + j] != U_4;
591 nd = MAX_DIGITS;
592 }
593 e1 = e -= nf;
594
595 //
596 // Now we have nd0 digits, starting at s0, followed by a
597 // decimal point, followed by nd-nd0 digits. The number we're
598 // after is the integer represented by those digits times 10**e
599 //
600 if (! nd0)
601 nd0 = nd;
602 k = nd < 8 ? nd : 8;
603 *this = IEEEfloat(y);
604 if (nd < 8)
605 {
606 if (!e)
607 {
608 if (sign)
609 *this = -*this;
610 return;
611 }
612 if (e > 0)
613 {
614 if (e <= 10)
615 {
616 *this *= tens[e];
617 if (sign)
618 *this = -*this;
619 return;
620 }
621 i = 7 - nd;
622 if (e <= 10 + i)
623 {
624 e -= i;
625 *this *= tens[i];
626 *this *= tens[e];
627 if (sign)
628 *this = -*this;
629 return;
630 }
631 }
632 else if (e >= -10 )
633 {
634 *this /= tens[-e];
635 if (sign)
636 *this = -*this;
637 return;
638 }
639 }
640 e1 += nd - k;
641
642 //
643 // Get starting approximation: *this * 10**e1
644 //
645 if (e1 > 0)
646 {
647 i = e1 & 7;
648 if (i)
649 *this *= tens[i];
650 if (e1 >>= 3)
651 {
652 if (e1 > MAX_DEC_EXP >> 3)
653 {
654 *this = check_invalid ? NaN()
655 : sign ? NEGATIVE_INFINITY()
656 : POSITIVE_INFINITY();
657 return;
658 }
659 for (j = 0; e1 > 1; j++, e1 >>= 1)
660 if (e1 & 1)
661 *this *= bigtens[j];
662 //
663 // The last multiplication could overflow.
664 //
665 value.word -= (FRACT_SIZE + 1) * MIN_FRACT;
666 *this *= bigtens[j];
667 z = Exponent();
668 if (z > BIAS - FRACT_SIZE)
669 {
670 *this = check_invalid ? NaN()
671 : sign ? NEGATIVE_INFINITY()
672 : POSITIVE_INFINITY();
673 return;
674 }
675 if (z > BIAS - FRACT_SIZE - 1)
676 value.word = POS_INF - 1;
677 else
678 value.word += (FRACT_SIZE + 1) * MIN_FRACT;
679 }
680 }
681 else if (e1 < 0)
682 {
683 e1 = -e1;
684 i = e1 & 7;
685 if (i)
686 *this /= tens[i];
687 if (e1 >>= 3)
688 {
689 if (e1 >= 1 << 3)
690 {
691 *this = check_invalid ? NaN()
692 : sign ? NEGATIVE_ZERO()
693 : POSITIVE_ZERO();
694 return;
695 }
696 for (j = 0; e1 > 1; j++, e1 >>= 1)
697 if (e1 & 1)
698 *this /= bigtens[j];
699 //
700 // The last multiplication could underflow.
701 //
702 IEEEfloat tmp(*this);
703 *this /= bigtens[j];
704 if (IsZero())
705 {
706 *this = tmp * 2;
707 *this /= bigtens[j];
708 if (IsZero())
709 {
710 *this = check_invalid ? NaN()
711 : sign ? NEGATIVE_ZERO()
712 : POSITIVE_ZERO();
713 return;
714 }
715 value.word = 1;
716 }
717 }
718 }
719
720 //
721 // Now the hard part -- adjusting *this to the correct value.
722 // Put digits into bd: true value = bd * 10^e
723 //
724 BigInt bd0(s0, nd0, nd, y, 8);
725 if (roundup)
726 ++bd0;
727 while (true) {
728 BigInt bd(bd0);
729 BigInt bb(*this, bbe, bbbits); // *this = bb * 2^bbe
730 BigInt bs(1);
731 if (e >= 0)
732 {
733 bb2 = bb5 = 0;
734 bd2 = bd5 = e;
735 }
736 else
737 {
738 bb2 = bb5 = -e;
739 bd2 = bd5 = 0;
740 }
741 if (bbe >= 0)
742 bb2 += bbe;
743 else
744 bd2 -= bbe;
745 bs2 = bb2;
746 j = bbe;
747 i = j + bbbits - 1; // logb(*this)
748 if (i < 1 - BIAS) // denormal
749 j += BIAS + FRACT_SIZE;
750 else
751 j = FRACT_SIZE + 2 - bbbits;
752 bb2 += j;
753 bd2 += j;
754 i = bb2 < bd2 ? bb2 : bd2;
755 if (i > bs2)
756 i = bs2;
757 if (i > 0)
758 {
759 bb2 -= i;
760 bd2 -= i;
761 bs2 -= i;
762 }
763 if (bb5 > 0) {
764 bs.pow5mult(bb5);
765 bb *= bs;
766 }
767 if (bb2 > 0)
768 bb <<= bb2;
769 if (bd5 > 0)
770 bd.pow5mult(bd5);
771 if (bd2 > 0)
772 bd <<= bd2;
773 if (bs2 > 0)
774 bs <<= bs2;
775 BigInt delta = bb - bd;
776 dsign = delta.IsNegative();
777 if (dsign)
778 delta.IsNegative(false);
779 i = delta.compareTo(bs);
780 //
781 // Error is less than half an ulp -- check for
782 // special case of mantissa a power of two.
783 //
784 if (i < 0)
785 {
786 if (dsign || FractBits() || Exponent() <= 1 - BIAS
787 || delta.IsZero())
788 break;
789 delta <<= 1;
790 if (delta.compareTo(bs) > 0)
791 // boundary case -- decrement exponent
792 value.word--;
793 break;
794 }
795 //
796 // exactly half-way between
797 //
798 else if (i == 0)
799 {
800 if (value.word & 1)
801 value.word += dsign ? 1 : -1;
802 break;
803 }
804 //
805 // more than 1/2 ulp off - try again
806 //
807 // This is broken into a separate method because mingw gcc 2.95.2
808 // has an ICE caused by register over-allocation if it is inline.
809 //
810 if (Adjust(delta, bs, dsign))
811 break;
812 }
813 if (check_invalid && (IsZero() || IsInfinite()))
814 *this = NaN();
815 else if (sign)
816 *this = -*this;
817 }
818
IntValue() const819 i4 IEEEfloat::IntValue() const
820 {
821 if (IsNaN())
822 return 0;
823
824 int sign = Sign(),
825 exponent = Exponent();
826
827 if (exponent > 30)
828 return sign ? Int::MIN_INT() : Int::MAX_INT();
829
830 // This covers true zero and denorms.
831 if (exponent < 0)
832 return 0;
833
834 i4 result = Fraction();
835
836 if (exponent > FRACT_SIZE)
837 result <<= (exponent - FRACT_SIZE);
838 else if (exponent < FRACT_SIZE)
839 result >>= (FRACT_SIZE - exponent);
840
841 return sign ? -result : result;
842 }
843
LongValue() const844 LongInt IEEEfloat::LongValue() const
845 {
846 if (IsNaN())
847 return LongInt(0);
848
849 int sign = Sign(),
850 exponent = Exponent();
851
852 if (exponent > 62)
853 return sign ? LongInt::MIN_LONG() : LongInt::MAX_LONG();
854
855 // This covers true zero and denorms.
856 if (exponent < 0)
857 return LongInt(0);
858
859 LongInt result(Fraction());
860
861 if (exponent > FRACT_SIZE)
862 result <<= (exponent - FRACT_SIZE);
863 else if (exponent < FRACT_SIZE)
864 result >>= (FRACT_SIZE - exponent);
865
866 return sign ? (LongInt) -result : result;
867 }
868
Normalize(int sign,int exponent,u4 fraction)869 IEEEfloat IEEEfloat::Normalize(int sign, int exponent, u4 fraction)
870 {
871 bool round = false, sticky = false;
872
873 assert(fraction != 0);
874
875 //
876 // Normalize right. Shift until value < MAX_FRACT.
877 //
878 if (fraction >= MAX_FRACT)
879 {
880 while (fraction >= MAX_FRACT)
881 {
882 sticky |= round;
883 round = (fraction & 1) != 0;
884 fraction >>= 1;
885 exponent++;
886 }
887 if (round && (sticky || (fraction & 1)) && exponent > -BIAS)
888 //
889 // Capture any overflow caused by rounding. No other checks are
890 // required because if overflow occurred, the the low order bit
891 // was guaranteed to be zero. Do not round denorms yet.
892 //
893 if (++fraction >= MAX_FRACT)
894 {
895 fraction >>= 1;
896 exponent++;
897 }
898 }
899
900 //
901 // Normalize left. Shift until value >= MIN_FRACT.
902 //
903 else
904 while (fraction < MIN_FRACT)
905 {
906 fraction <<= 1;
907 exponent--;
908 }
909
910 //
911 // Check and respond to overflow
912 //
913 if (exponent > BIAS)
914 return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
915
916 //
917 // Check and respond to underflow
918 //
919 if (exponent <= -BIAS)
920 {
921 if (exponent < -BIAS - FRACT_SIZE)
922 return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
923 while (exponent <= -BIAS)
924 {
925 sticky |= round;
926 round = (fraction & 1) != 0;
927 fraction >>= 1;
928 exponent++;
929 }
930 if (round && (sticky || (fraction & 1)))
931 fraction++;
932 exponent = -BIAS;
933 if (fraction == 0)
934 return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
935 }
936
937 fraction &= FRACT_BITS;
938 fraction |= ((exponent + BIAS) << FRACT_SIZE);
939 if (sign)
940 fraction |= SIGN_BIT;
941 return IEEEfloat(fraction);
942 }
943
SplitInto(u4 & fraction) const944 int IEEEfloat::SplitInto(u4 &fraction) const
945 {
946 int exponent = Exponent();
947 fraction = Fraction();
948
949 if (exponent == -BIAS)
950 {
951 exponent++;
952 while (fraction < MIN_FRACT)
953 {
954 fraction <<= 1;
955 exponent--;
956 }
957 }
958
959 return exponent;
960 }
961
Ulp() const962 IEEEfloat IEEEfloat::Ulp() const
963 {
964 i4 L;
965 IEEEfloat f;
966 f.value.float_value = value.float_value;
967
968 L = (i4) f.ExpBits() - FRACT_SIZE * MIN_FRACT;
969 if (L > 0)
970 f.value.iword = L;
971 else
972 {
973 L = -L >> FRACT_SIZE;
974 f.value.iword = L >= (i4) FRACT_SIZE ? 1 : 0x400000 >> L;
975 }
976 return f;
977 }
978
Ratio(const BigInt & a,const BigInt & b)979 IEEEfloat IEEEfloat::Ratio(const BigInt &a, const BigInt &b)
980 {
981 IEEEfloat fa, fb;
982 int k;
983 fa = a.FloatValue();
984 fb = b.FloatValue();
985 k = b.hi0bits() - a.hi0bits() + 32 * (a.wds - b.wds);
986 if (k > 0)
987 fa.value.word += k * MIN_FRACT;
988 else
989 fb.value.word -= k * MIN_FRACT;
990 return fa / fb;
991 }
992
operator ==(const IEEEfloat op) const993 bool IEEEfloat::operator== (const IEEEfloat op) const
994 {
995 // FIXME: This could be sped up by inlining
996 #ifdef HAVE_IEEE754
997 return value.float_value == op.value.float_value;
998 #else
999 return (IsNaN() || op.IsNaN() ? false
1000 : IsZero() && op.IsZero() ? true
1001 : value.word == op.value.word);
1002 #endif
1003 }
1004
operator !=(const IEEEfloat op) const1005 bool IEEEfloat::operator!= (const IEEEfloat op) const
1006 {
1007 // FIXME: This could be sped up by inlining
1008 #ifdef HAVE_IEEE754
1009 return value.float_value != op.value.float_value;
1010 #else
1011 return ! (*this == op);
1012 #endif
1013 }
1014
operator <(const IEEEfloat op) const1015 bool IEEEfloat::operator< (const IEEEfloat op) const
1016 {
1017 // FIXME: This could be sped up by inlining
1018 #ifdef HAVE_IEEE754
1019 return (value.float_value < op.value.float_value);
1020 #else
1021 if (IsNaN() || op.IsNaN())
1022 return false; // NaNs are unordered
1023 if (IsZero() && op.IsZero())
1024 return false;
1025 // Exploit fact that all other IEEE floating point numbers sort like
1026 // ints after worrying about sign.
1027 if (IsNegative())
1028 return op.IsPositive() ||
1029 (value.word & ABS_BITS) > (op.value.word & ABS_BITS);
1030 return op.IsPositive() &&
1031 (value.word & ABS_BITS) < (op.value.word & ABS_BITS);
1032 #endif
1033 }
1034
operator <=(const IEEEfloat op) const1035 bool IEEEfloat::operator<= (const IEEEfloat op) const
1036 {
1037 // FIXME: This could be sped up by inlining
1038 #ifdef HAVE_IEEE754
1039 return (value.float_value <= op.value.float_value);
1040 #else
1041 return *this < op || *this == op;
1042 #endif
1043 }
1044
operator >(const IEEEfloat op) const1045 bool IEEEfloat::operator> (const IEEEfloat op) const
1046 {
1047 // FIXME: This could be sped up by inlining
1048 #ifdef HAVE_IEEE754
1049 return (value.float_value > op.value.float_value);
1050 #else
1051 if (IsNaN() || op.IsNaN())
1052 return false; // NaNs are unordered.
1053 if (IsZero() && op.IsZero())
1054 return false;
1055 // Exploit fact that all other IEEE floating point numbers sort like
1056 // ints after worrying about sign.
1057 if (IsPositive())
1058 return op.IsNegative() ||
1059 (value.word & ABS_BITS) > (op.value.word & ABS_BITS);
1060 return op.IsNegative() &&
1061 (value.word & ABS_BITS) < (op.value.word & ABS_BITS);
1062 #endif
1063 }
1064
operator >=(const IEEEfloat op) const1065 bool IEEEfloat::operator>= (const IEEEfloat op) const
1066 {
1067 // FIXME: This could be sped up by inlining
1068 #ifdef HAVE_IEEE754
1069 return (value.float_value >= op.value.float_value);
1070 #else
1071 return *this > op || *this == op;
1072 #endif
1073 }
1074
operator +(const IEEEfloat op) const1075 IEEEfloat IEEEfloat::operator+ (const IEEEfloat op) const
1076 {
1077 #ifdef HAVE_IEEE754
1078 // FIXME: This could be sped up by inlining
1079 return IEEEfloat(value.float_value + op.value.float_value);
1080 #else
1081 if (IsNaN() || op.IsNaN())
1082 return NaN(); // arithmetic on NaNs not allowed
1083
1084 //
1085 // Adding unlike infinities not allowed
1086 // Adding infinities of same sign is infinity of that sign
1087 // Adding finite and infinity produces infinity
1088 //
1089 if (IsInfinite())
1090 return op.IsInfinite() && (Sign() != op.Sign()) ? NaN() : *this;
1091 if (op.IsInfinite())
1092 return op;
1093
1094 //
1095 // Adding zero is easy
1096 //
1097 if (IsZero())
1098 return (op.IsZero()) ? (Sign() != op.Sign()) ? POSITIVE_ZERO()
1099 : *this
1100 : op;
1101 if (op.IsZero())
1102 return *this;
1103
1104 //
1105 // Now for the real work - do manipulations on copies
1106 //
1107 i4 x, y, round = 0;
1108 int expx, expy, signx, signy;
1109
1110 expx = SplitInto((u4 &) x);
1111 expy = op.SplitInto((u4 &) y);
1112
1113 signx = Sign();
1114 signy = op.Sign();
1115
1116 // If the exponents are far enough apart, the answer is easy
1117 if (expx > expy + 25)
1118 return *this;
1119 if (expy > expx + 25)
1120 return op;
1121
1122 //
1123 // Denormalize the fractions, so that the exponents are
1124 // the same and then set the exponent for the result.
1125 // Leave enough space for overflow and INT_MIN avoidance!
1126 //
1127 if (signx)
1128 x = -x;
1129 if (signy)
1130 y = -y;
1131
1132 x <<= 6;
1133 y <<= 6;
1134
1135 if (expx > expy) {
1136 round = y << (32 + expy - expx);
1137 y >>= expx - expy;
1138 }
1139 else if (expy > expx)
1140 {
1141 round = x << (32 + expx - expy);
1142 x >>= expy - expx;
1143 expx = expy;
1144 }
1145
1146 //
1147 // Do the arithmetic. The excess magnitude of 32-bit arithmetic means
1148 // overflow is impossible (we only need 1 spare bit!). We ensure that
1149 // pre-alignment avoids any question of INT_MIN negation problems.
1150 //
1151 x += y;
1152 if (round)
1153 x |= 1;
1154
1155 //
1156 // If the result is negative, then make the fraction positive again
1157 // and remember the sign.
1158 //
1159 if (x < 0)
1160 {
1161 x = -x;
1162 signx = 1;
1163 }
1164 else
1165 signx = 0;
1166
1167 if (x == 0)
1168 return signx ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1169
1170 //
1171 // Time to normalize again! If we need to shift left (the addition was
1172 // effectively a subtraction), then there cannot be any reason to round.
1173 // If the number fits exactly we don't have anything to do either.
1174 //
1175 return Normalize(signx, expx - 6, (u4) x);
1176 #endif
1177 }
1178
operator -() const1179 IEEEfloat IEEEfloat::operator- () const
1180 {
1181 // FIXME: This could be sped up by inlining
1182 #ifdef HAVE_IEEE754
1183 return IEEEfloat(-value.float_value);
1184 #else
1185 if (IsNaN())
1186 return *this;
1187 return IEEEfloat(value.word ^ SIGN_BIT);
1188 #endif
1189 }
1190
operator *(const IEEEfloat op) const1191 IEEEfloat IEEEfloat::operator* (const IEEEfloat op) const
1192 {
1193 #ifdef HAVE_IEEE754
1194 return IEEEfloat(value.float_value * op.value.float_value);
1195 #else
1196 if (IsNaN() || op.IsNaN())
1197 return NaN(); // arithmetic on NaNs not allowed
1198
1199 int sign = Sign() ^ op.Sign();
1200
1201 //
1202 // If either operand is zero or infinite, then the answer is easy.
1203 //
1204 if (IsZero())
1205 return op.IsInfinite() ? NaN()
1206 : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1207
1208 if (op.IsZero())
1209 return IsInfinite() ? NaN()
1210 : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1211
1212 if (IsInfinite() || op.IsInfinite())
1213 return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
1214
1215 //
1216 // Now for the real work - do manipulations on copies
1217 //
1218 u4 x, y;
1219 int exponent;
1220
1221 exponent = SplitInto(x) + op.SplitInto(y);
1222
1223 //
1224 // The numbers to be multiplied are 24 bits in length (stored in 32 bit
1225 // integers). Using ULongInt to perform the multiplication, the result
1226 // will be 46-48 bits (unsigned); shift it back to 28 bits for Normalize,
1227 // while folding the low 20 bits into the lsb for rounding purposes.
1228 //
1229
1230 ULongInt a = x,
1231 b = y;
1232 a *= b;
1233 b = a & 0xfffff;
1234 a >>= 20;
1235 x = a.LowWord() | ((b > 0) ? 1 : 0);
1236
1237 return Normalize(sign, exponent - 3, x);
1238 #endif // HAVE_IEEE754
1239 }
1240
operator /(const IEEEfloat op) const1241 IEEEfloat IEEEfloat::operator/ (const IEEEfloat op) const
1242 {
1243 #ifdef HAVE_IEEE754
1244 return op.IsZero() ? ((IsNaN() || IsZero()) ? NaN()
1245 : (IsPositive() ^ op.IsPositive())
1246 ? NEGATIVE_INFINITY()
1247 : POSITIVE_INFINITY())
1248 : IEEEfloat(value.float_value / op.value.float_value);
1249 #else // HAVE_IEEE754
1250 if (IsNaN() || op.IsNaN())
1251 return NaN(); // arithmetic on NaNs not allowed
1252
1253 int sign = Sign() ^ op.Sign();
1254
1255 //
1256 // Infinities and zeroes are special.
1257 //
1258
1259 if (IsInfinite())
1260 return (op.IsInfinite() ? NaN()
1261 : sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY());
1262
1263 if (op.IsInfinite())
1264 return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1265
1266 if (IsZero())
1267 return op.IsZero() ? NaN()
1268 : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1269
1270 if (op.IsZero())
1271 return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
1272
1273 //
1274 // Now for the real work - do manipulations on copies
1275 //
1276 u4 x, y;
1277 int exponent;
1278
1279 exponent = SplitInto(x) - op.SplitInto(y);
1280
1281 u4 mask = 0x80000000,
1282 result = 0;
1283
1284 // Because both values are normalised, a single shift guarantees results.
1285
1286 if (x < y)
1287 {
1288 x <<= 1;
1289 exponent--;
1290 }
1291
1292 //
1293 // If the numerator is larger, then it is divisible.
1294 // Reflect this in the result, and do the subtraction.
1295 // Magnify the numerator again and reduce the mask.
1296 //
1297 while (mask)
1298 {
1299 if (x >= y)
1300 {
1301 result += mask;
1302 x -= y;
1303 if (x == 0)
1304 break;
1305 }
1306 x <<= 1;
1307 mask >>= 1;
1308 }
1309
1310 return Normalize(sign, exponent - 8, result);
1311 #endif // HAVE_IEEE754
1312 }
1313
operator %(const IEEEfloat op) const1314 IEEEfloat IEEEfloat::operator% (const IEEEfloat op) const
1315 {
1316 #ifdef HAVE_IEEE754
1317 return IEEEfloat((op.IsZero() ? NaN().value.float_value
1318 : (float) fmod((double) value.float_value,
1319 (double) op.value.float_value)));
1320 #else // HAVE_IEEE754
1321 if (IsNaN() || op.IsNaN())
1322 return NaN(); // arithmetic on NaNs not allowed
1323
1324 //
1325 // Infinities and zeroes are special.
1326 //
1327
1328 if (IsInfinite() || op.IsZero())
1329 return NaN();
1330
1331 if (IsZero() || op.IsInfinite())
1332 return *this;
1333
1334 //
1335 // Now for the real work - do manipulations on copies
1336 // This algorithm is from fdlibm.c - see above notice
1337 //
1338
1339 int sign = Sign();
1340 u4 x, y;
1341 int expy, exponent;
1342 i4 z;
1343
1344 expy = op.SplitInto(y);
1345 exponent = SplitInto(x) - expy;
1346
1347 if (exponent < 0 || (exponent == 0 && x < y))
1348 return *this;
1349
1350 while (exponent--)
1351 {
1352 z = x - y;
1353 if (z < 0)
1354 x <<= 1;
1355 else if (z == 0)
1356 return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1357 else
1358 x = z + z;
1359 }
1360 z = x - y;
1361 if (z >= 0)
1362 x = z;
1363
1364 if (x == 0)
1365 return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1366 return Normalize(sign, expy, x);
1367 #endif // HAVE_IEEE754
1368 }
1369
1370
1371
IEEEdouble(const IEEEfloat & f)1372 IEEEdouble::IEEEdouble(const IEEEfloat &f)
1373 {
1374 #ifdef HAVE_IEEE754
1375 value.double_value = (double) f.FloatView();
1376 #else
1377 int sign = f.Sign();
1378 int exponent = f.Exponent();
1379
1380 if (exponent == -IEEEfloat::Bias())
1381 {
1382 if (f.IsZero())
1383 *this = sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1384 else
1385 {
1386 //
1387 // This is a denormalized number, shift it to fit double format
1388 //
1389 *this = Normalize(sign, 1 - IEEEfloat::Bias(),
1390 ULongInt(f.Fraction()) << 29);
1391 }
1392 }
1393 else
1394 {
1395 if (f.IsInfinite())
1396 *this = sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
1397 else if (f.IsNaN())
1398 *this = NaN();
1399 else
1400 {
1401 // Regular, normalized number. Shift it to fit double format
1402 *this = Normalize(sign, exponent, ULongInt(f.Fraction()) << 29);
1403 }
1404 }
1405 #endif // HAVE_IEEE754
1406 }
1407
IEEEdouble(i4 a)1408 IEEEdouble::IEEEdouble(i4 a)
1409 {
1410 #ifdef HAVE_IEEE754
1411 value.double_value = (double) a;
1412 #else
1413 int sign = 0;
1414
1415 if (a < 0)
1416 {
1417 a = -a; // even works for MIN_INT!
1418 sign = 1;
1419 }
1420 *this = (a == 0) ? POSITIVE_ZERO()
1421 : Normalize(sign, FRACT_SIZE, ULongInt((u4) a));
1422 #endif // HAVE_IEEE754
1423 }
1424
IEEEdouble(const LongInt & a)1425 IEEEdouble::IEEEdouble(const LongInt &a)
1426 {
1427 #ifdef HAVE_IEEE754
1428 # ifdef HAVE_64BIT_TYPES
1429 value.double_value = (double)(i8) a.Words();
1430 # else
1431 value.double_value = ((double)(i4) a.HighWord() * (double) 0x40000000 *
1432 4.0) + (double) a.LowWord();
1433 # endif // HAVE_64BIT_TYPES
1434 #else
1435 int sign = 0;
1436 LongInt l = a;
1437
1438 if (a < 0)
1439 {
1440 l = -a; // even works for MIN_LONG!
1441 sign = 1;
1442 }
1443 *this = (l == 0) ? POSITIVE_ZERO()
1444 : Normalize(sign, FRACT_SIZE, ULongInt(l));
1445 #endif // HAVE_IEEE754
1446 }
1447
IEEEdouble(const char * str,bool check_invalid)1448 IEEEdouble::IEEEdouble(const char *str, bool check_invalid)
1449 {
1450 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2,
1451 e, e1, i, j, k;
1452 int nd, // number of digits in mantissa (except extra '0's)
1453 nd0, // number of digits before '.' (except leading '0's)
1454 nf, // number of digits after '.' (except trailing '0's)
1455 nz; // number of consecutive '0' digits
1456 bool nz0, // whether leading zero exists
1457 roundup, // whether to round up long string
1458 sign, // if the string represents a negative
1459 dsign, // the sign of delta
1460 esign; // the sign of the exponent
1461 char c;
1462 const char *s, *s0, *s1;
1463 IEEEdouble aadj, aadj1;
1464 i4 L;
1465 i4 y, z;
1466
1467 sign = nz0 = roundup = false;
1468 nz = 0;
1469
1470 //
1471 // consume whitespace
1472 //
1473 for (s = str; ; s++)
1474 {
1475 switch (*s)
1476 {
1477 case U_MINUS:
1478 sign = true;
1479 // fallthrough
1480 case U_PLUS:
1481 if (*++s)
1482 goto break2;
1483 // fallthrough
1484 case U_NU:
1485 *this = NaN();
1486 return;
1487 case U_SPACE:
1488 case U_HT:
1489 case U_FF:
1490 case U_LF:
1491 case U_CR:
1492 continue;
1493 default:
1494 goto break2;
1495 }
1496 }
1497 break2:
1498 //
1499 // Consume leading zeroes.
1500 //
1501 if (*s == U_0)
1502 {
1503 nz0 = true;
1504 while (*++s == U_0);
1505
1506 if (! *s)
1507 {
1508 *this = sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1509 return;
1510 }
1511 //
1512 // Parse hexadecimal floating point.
1513 //
1514 else if (*s == U_x || *s == U_X)
1515 {
1516 LongInt fraction = 0;
1517 // Exponent adjustment, based on '.' location.
1518 int shift = FRACT_SIZE;
1519 bool seen_dot = false;
1520 while (*++s && *s == U_0); // Consume leading zeroes.
1521 c = *s;
1522 if (c == U_DOT)
1523 {
1524 seen_dot = true;
1525 while (*++s && *s == U_0)
1526 shift -= 4;
1527 c = *s;
1528 }
1529 if (! c || c == U_p || c == U_P)
1530 {
1531 *this = sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1532 return;
1533 }
1534 // To avoid overflow, stop after enough bits have been read.
1535 for (i = 0; i < (FRACT_SIZE >> 2) + 2 && c; i++, c = *++s)
1536 {
1537 if (c == U_DOT)
1538 {
1539 if (seen_dot)
1540 break;
1541 c = *++s;
1542 seen_dot = true;
1543 }
1544 int value;
1545 if (Code::IsHexDigit(c))
1546 value = Code::Value(c);
1547 else break;
1548 if (seen_dot)
1549 shift -= 4;
1550 fraction = (fraction << 4) + value;
1551 }
1552 // Round any remaining bits.
1553 bool sticky = false;
1554 while (c == U_DOT || Code::IsHexDigit(c))
1555 {
1556 if (c == U_DOT)
1557 {
1558 if (seen_dot)
1559 break;
1560 seen_dot = true;
1561 }
1562 else
1563 {
1564 if (! seen_dot)
1565 shift += 4;
1566 if (c != U_0)
1567 sticky = true;
1568 }
1569 c = *++s;
1570 }
1571 assert(fraction != 0);
1572 if (sticky)
1573 fraction++;
1574 // On to the expononet.
1575 int exponent = 0;
1576 esign = false;
1577 if (c == U_p || c == U_P)
1578 {
1579 if (*++s == U_MINUS)
1580 {
1581 esign = true;
1582 s++;
1583 }
1584 else if (*s == U_PLUS)
1585 s++;
1586 while ((c = *s++))
1587 {
1588 if (! Code::IsDecimalDigit(c))
1589 break;
1590 exponent = exponent * 10 + c - U_0;
1591 // Check for exponent overflow
1592 if (exponent + shift > 19999)
1593 {
1594 if (check_invalid)
1595 *this = NaN();
1596 else
1597 {
1598 *this = esign ? POSITIVE_ZERO()
1599 : POSITIVE_INFINITY();
1600 if (sign)
1601 *this = - *this;
1602 }
1603 return;
1604 }
1605 }
1606 }
1607 if (esign)
1608 exponent = - exponent;
1609 *this = Normalize(sign, exponent + shift, fraction);
1610 if (check_invalid && (IsZero() || IsInfinite()))
1611 *this = NaN();
1612 return;
1613 }
1614 }
1615
1616 //
1617 // parse before '.'
1618 //
1619 s0 = s;
1620 y = z = 0;
1621 for (nd = nf = 0; Code::IsDecimalDigit(c = *s); nd++, s++)
1622 if (nd < 9)
1623 y = 10 * y + c - U_0;
1624 else if (nd < 16)
1625 z = 10 * z + c - U_0;
1626 nd0 = nd;
1627 if (c == U_DOT)
1628 {
1629 //
1630 // parse after '.'
1631 //
1632 c = *++s;
1633 if (!nd)
1634 {
1635 for ( ; c == U_0; c = *++s)
1636 nz++;
1637 if (c > U_0 && c <= U_9)
1638 {
1639 s0 = s;
1640 nf += nz;
1641 nz = 0;
1642 }
1643 }
1644 for ( ; Code::IsDecimalDigit(c); c = *++s)
1645 {
1646 nz++;
1647 if (c -= U_0)
1648 {
1649 nf += nz;
1650 for (i = 1; i < nz; i++)
1651 if (nd++ < 9)
1652 y *= 10;
1653 else if (nd <= 16)
1654 z *= 10;
1655 if (nd++ < 9)
1656 y = 10 * y + c;
1657 else if (nd <= 16)
1658 z = 10 * z + c;
1659 nz = 0;
1660 }
1661 }
1662 }
1663 //
1664 // consume exponent
1665 //
1666 e = 0;
1667 if (c == U_e || c == U_E) {
1668 str = s;
1669 esign = false;
1670 switch (c = *++s)
1671 {
1672 case U_MINUS:
1673 esign = true;
1674 // fallthrough
1675 case U_PLUS:
1676 c = *++s;
1677 }
1678 if (Code::IsDecimalDigit(c))
1679 {
1680 if (!nd && !nz && !nz0)
1681 {
1682 *this = sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1683 return;
1684 }
1685 while (c == U_0)
1686 c = *++s;
1687 if (c > U_0 && c <= U_9)
1688 {
1689 L = c - U_0;
1690 s1 = s;
1691 while (Code::IsDecimalDigit(c = *++s))
1692 L = 10 * L + c - U_0;
1693 //
1694 // Avoid confusion from exponents so large that e might
1695 // overflow
1696 if (s - s1 > 8 || L > 19999)
1697 e = 19999;
1698 else
1699 e = L;
1700 if (esign)
1701 e = -e;
1702 }
1703 else
1704 e = 0;
1705 }
1706 else
1707 s = str;
1708 }
1709 if (!nd) {
1710 *this = (!nz && !nz0) ? NaN()
1711 : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
1712 return;
1713 }
1714 //
1715 // for long strings, round away all digits beyond maximum precise string
1716 // for fraction, there are n decimal digits after '.' if lsb is 2^n;
1717 // but the first m of these digits are '0', for d = x*10^m
1718 // so, digits after MAX_DIGITS may be ignored
1719 //
1720 if (nd > MAX_DIGITS)
1721 {
1722 k = nd - MAX_DIGITS;
1723 i = MAX_DIGITS - nd0;
1724 if (i <= 0)
1725 {
1726 // decimal after last precise digit
1727 nd0 = MAX_DIGITS;
1728 nf = i;
1729 j = 0;
1730 }
1731 else if (i == MAX_DIGITS)
1732 {
1733 // decimal before first precise digit
1734 nf -= k;
1735 j = 0;
1736 }
1737 else
1738 {
1739 // decimal inside precise digits
1740 nf -= k;
1741 j = 1;
1742 }
1743 roundup = s0[MAX_DIGITS - 1 + j] != U_4;
1744 nd = MAX_DIGITS;
1745 }
1746 e1 = e -= nf;
1747
1748 //
1749 // Now we have nd0 digits, starting at s0, followed by a
1750 // decimal point, followed by nd-nd0 digits. The number we're
1751 // after is the integer represented by those digits times 10**e
1752 //
1753 if (! nd0)
1754 nd0 = nd;
1755 k = nd < 16 ? nd : 16;
1756 *this = IEEEdouble(y);
1757 if (k > 9)
1758 *this = *this * tens[k - 9] + z;
1759 if (nd < 16)
1760 {
1761 if (!e)
1762 {
1763 if (sign)
1764 *this = -*this;
1765 return;
1766 }
1767 if (e > 0)
1768 {
1769 if (e <= 22)
1770 {
1771 *this *= tens[e];
1772 if (sign)
1773 *this = -*this;
1774 return;
1775 }
1776 i = 15 - nd;
1777 if (e <= 22 + i)
1778 {
1779 e -= i;
1780 *this *= tens[i];
1781 *this *= tens[e];
1782 if (sign)
1783 *this = -*this;
1784 return;
1785 }
1786 }
1787 else if (e >= -22 )
1788 {
1789 *this /= tens[-e];
1790 if (sign)
1791 *this = -*this;
1792 return;
1793 }
1794 }
1795 e1 += nd - k;
1796
1797 //
1798 // Get starting approximation: *this * 10**e1
1799 //
1800 if (e1 > 0)
1801 {
1802 i = e1 & 0xf;
1803 if (i)
1804 *this *= tens[i];
1805 if (e1 >>= 4)
1806 {
1807 if (e1 > MAX_DEC_EXP >> 4)
1808 {
1809 *this = check_invalid ? NaN()
1810 : sign ? NEGATIVE_INFINITY()
1811 : POSITIVE_INFINITY();
1812 return;
1813 }
1814 for (j = 0; e1 > 1; j++, e1 >>= 1)
1815 if (e1 & 1)
1816 *this *= bigtens[j];
1817 //
1818 // The last multiplication could overflow.
1819 //
1820 setHighWord(HighWord() - (FRACT_SIZE + 1) * MIN_FRACT);
1821 *this *= bigtens[j];
1822 z = Exponent();
1823 if (z > BIAS - FRACT_SIZE)
1824 {
1825 *this = check_invalid ? NaN()
1826 : sign ? NEGATIVE_INFINITY()
1827 : POSITIVE_INFINITY();
1828 return;
1829 }
1830 if (z > BIAS - FRACT_SIZE - 1)
1831 setHighAndLowWords(POS_INF_HI - 1, ZERO_LO - 1);
1832 else
1833 setHighWord(HighWord() + (FRACT_SIZE + 1) * MIN_FRACT);
1834 }
1835 }
1836 else if (e1 < 0)
1837 {
1838 e1 = -e1;
1839 i = e1 & 0xf;
1840 if (i)
1841 *this /= tens[i];
1842 if (e1 >>= 4)
1843 {
1844 if (e1 >= 1 << 5)
1845 {
1846 *this = check_invalid ? NaN()
1847 : sign ? NEGATIVE_ZERO()
1848 : POSITIVE_ZERO();
1849 return;
1850 }
1851 for (j = 0; e1 > 1; j++, e1 >>= 1)
1852 if (e1 & 1)
1853 *this /= bigtens[j];
1854 //
1855 // The last multiplication could underflow.
1856 //
1857 IEEEdouble tmp(*this);
1858 *this /= bigtens[j];
1859 if (IsZero())
1860 {
1861 *this = tmp * 2;
1862 *this /= bigtens[j];
1863 if (IsZero())
1864 {
1865 *this = check_invalid ? NaN()
1866 : sign ? NEGATIVE_ZERO()
1867 : POSITIVE_ZERO();
1868 return;
1869 }
1870 setHighAndLowWords(0, 1);
1871 }
1872 }
1873 }
1874
1875 //
1876 // Now the hard part -- adjusting *this to the correct value.
1877 // Put digits into bd: true value = bd * 10^e
1878 //
1879 BigInt bd0(s0, nd0, nd, y, 9);
1880 if (roundup)
1881 ++bd0;
1882 while (true) {
1883 BigInt bd(bd0);
1884 BigInt bb(*this, bbe, bbbits); // *this = bb * 2^bbe
1885 BigInt bs(1);
1886 if (e >= 0)
1887 {
1888 bb2 = bb5 = 0;
1889 bd2 = bd5 = e;
1890 }
1891 else
1892 {
1893 bb2 = bb5 = -e;
1894 bd2 = bd5 = 0;
1895 }
1896 if (bbe >= 0)
1897 bb2 += bbe;
1898 else
1899 bd2 -= bbe;
1900 bs2 = bb2;
1901 j = bbe;
1902 i = j + bbbits - 1; // logb(*this)
1903 if (i < 1 - BIAS) // denormal
1904 j += BIAS + FRACT_SIZE;
1905 else
1906 j = FRACT_SIZE + 2 - bbbits;
1907 bb2 += j;
1908 bd2 += j;
1909 i = bb2 < bd2 ? bb2 : bd2;
1910 if (i > bs2)
1911 i = bs2;
1912 if (i > 0)
1913 {
1914 bb2 -= i;
1915 bd2 -= i;
1916 bs2 -= i;
1917 }
1918 if (bb5 > 0) {
1919 bs.pow5mult(bb5);
1920 bb *= bs;
1921 }
1922 if (bb2 > 0)
1923 bb <<= bb2;
1924 if (bd5 > 0)
1925 bd.pow5mult(bd5);
1926 if (bd2 > 0)
1927 bd <<= bd2;
1928 if (bs2 > 0)
1929 bs <<= bs2;
1930 BigInt delta = bb - bd;
1931 dsign = delta.IsNegative();
1932 if (dsign)
1933 delta.IsNegative(false);
1934 i = delta.compareTo(bs);
1935 //
1936 // Error is less than half an ulp -- check for
1937 // special case of mantissa a power of two.
1938 //
1939 if (i < 0)
1940 {
1941 if (dsign || LowWord() || FractBits()
1942 || Exponent() <= 1 - BIAS || delta.IsZero())
1943 break;
1944 delta <<= 1;
1945 if (delta.compareTo(bs) > 0)
1946 // boundary case -- decrement exponent
1947 BaseLong::operator --();
1948 break;
1949 }
1950 //
1951 // exactly half-way between
1952 //
1953 else if (i == 0)
1954 {
1955 if (LowWord() & 1)
1956 BaseLong::operator +=(dsign ? 1 : -1);
1957 break;
1958 }
1959 //
1960 // more than 1/2 ulp off - try again
1961 //
1962 aadj = Ratio(delta, bs);
1963 if (aadj <= 2)
1964 {
1965 if (dsign)
1966 aadj = aadj1 = 1;
1967 else if (FractBits() || LowWord())
1968 {
1969 if (!HighWord() && LowWord() == 1)
1970 {
1971 // underflow
1972 *this = POSITIVE_ZERO();
1973 break;
1974 }
1975 aadj = 1;
1976 aadj1 = -1;
1977 }
1978 else
1979 {
1980 //
1981 // special case - mantissa is power of 2
1982 //
1983 if (aadj < 1.0)
1984 aadj = 0.5;
1985 else
1986 aadj *= 0.5;
1987 aadj1 = -aadj;
1988 }
1989 }
1990 else
1991 {
1992 aadj *= 0.5;
1993 //
1994 // gcc 2.95.2 has an ICE from too many registers with this line:
1995 // aadj1 = dsign ? aadj : -aadj;
1996 //
1997 if (dsign)
1998 aadj1 = aadj;
1999 else
2000 aadj1 = -aadj;
2001 }
2002 y = Exponent();
2003 //
2004 // Check for overflow
2005 //
2006 if (y == BIAS)
2007 {
2008 IEEEdouble tmp(*this);
2009 setHighWord(HighWord() - (FRACT_SIZE + 1) * MIN_FRACT);
2010 *this += aadj1 * Ulp();
2011 if (Exponent() >= BIAS - FRACT_SIZE)
2012 {
2013 if (tmp.HighWord() == POS_INF_HI - 1 &&
2014 tmp.LowWord() == ZERO_LO - 1)
2015 {
2016 // overflow
2017 *this = POSITIVE_INFINITY();
2018 break;
2019 }
2020 setHighAndLowWords(POS_INF_HI - 1, ZERO_LO - 1);
2021 continue;
2022 }
2023 else
2024 setHighWord(HighWord() + (FRACT_SIZE + 1) * MIN_FRACT);
2025 }
2026 else
2027 {
2028 //
2029 // Compute adj so that the IEEE rounding rules will
2030 // correctly round *this + adj in some half-way cases.
2031 // If *this * Ulp() is denormalized, we must adjust aadj
2032 // to avoid trouble from bits lost to denormalization.
2033 //
2034 if (y <= FRACT_SIZE - BIAS || aadj > 1)
2035 {
2036 aadj1 = IEEEdouble((aadj + 0.5).IntValue());
2037 if (! dsign)
2038 aadj1 = -aadj1;
2039 }
2040 *this += aadj1 * Ulp();
2041 }
2042 if (y == Exponent())
2043 {
2044 //
2045 // Can we stop now?
2046 // The tolerances below are conservative.
2047 //
2048 aadj -= aadj.IntValue();
2049 if (dsign || FractBits() || LowWord())
2050 {
2051 if (aadj < .4999999 || aadj > .5000001)
2052 break;
2053 }
2054 else if (aadj < .4999999 / 2)
2055 break;
2056 }
2057 }
2058 if (check_invalid && (IsZero() || IsInfinite()))
2059 *this = NaN();
2060 else if (sign)
2061 *this = -*this;
2062 }
2063
IntValue() const2064 i4 IEEEdouble::IntValue() const
2065 {
2066 if (IsNaN())
2067 return 0;
2068
2069 #ifdef HAVE_IEEE754
2070 if (value.double_value < (double)(i4) Int::MIN_INT())
2071 return Int::MIN_INT();
2072 else if (value.double_value > (double) Int::MAX_INT())
2073 return Int::MAX_INT();
2074 return (i4) value.double_value;
2075 #else
2076 int sign = Sign(),
2077 exponent = Exponent();
2078
2079 if (exponent > 30)
2080 return sign ? Int::MIN_INT() : Int::MAX_INT();
2081
2082 // This includes true zero and denorms.
2083 if (exponent < 0)
2084 return 0;
2085
2086 i4 result = (i4) (Fraction() >> (FRACT_SIZE - exponent)).LowWord();
2087
2088 return sign ? -result : result;
2089 #endif // ! HAVE_IEEE754
2090 }
2091
LongValue() const2092 LongInt IEEEdouble::LongValue() const
2093 {
2094 if (IsNaN())
2095 return LongInt(0);
2096
2097 int sign = Sign(),
2098 exponent = Exponent();
2099
2100 if (exponent > 62)
2101 return sign ? LongInt::MIN_LONG() : LongInt::MAX_LONG();
2102
2103 // This covers true zero and denorms.
2104 if (exponent < 0)
2105 return LongInt(0);
2106
2107 LongInt result = Fraction();
2108
2109 if (exponent > (int) FRACT_SIZE)
2110 result <<= (exponent - FRACT_SIZE);
2111 else if (exponent < (int) FRACT_SIZE)
2112 result >>= (FRACT_SIZE - exponent);
2113
2114 return sign ? (LongInt) -result : result;
2115 }
2116
Normalize(int sign,int exponent,ULongInt fraction)2117 IEEEdouble IEEEdouble::Normalize(int sign, int exponent, ULongInt fraction)
2118 {
2119 bool round = false, sticky = false;
2120
2121 assert(fraction != 0);
2122
2123 //
2124 // Normalize right. Shift until value < MAX_FRACT.
2125 //
2126 if (fraction.HighWord() >= MAX_FRACT)
2127 {
2128 while (fraction.HighWord() >= MAX_FRACT)
2129 {
2130 sticky |= round;
2131 round = (fraction.LowWord() & 1) != 0;
2132 fraction >>= 1;
2133 exponent++;
2134 }
2135 if (round && (sticky || (fraction.LowWord() & 1)) &&
2136 exponent > -(int) BIAS)
2137 {
2138 //
2139 // Capture any overflow caused by rounding. No other checks are
2140 // required because if overflow occurred, the the low order bit
2141 // was guaranteed to be zero. Do not round denorms yet.
2142 //
2143 if ((++fraction).HighWord() >= MAX_FRACT)
2144 {
2145 fraction >>= 1;
2146 exponent++;
2147 }
2148 }
2149 }
2150
2151 //
2152 // Normalize left. Shift until value >= MIN_FRACT.
2153 //
2154 else
2155 while (fraction.HighWord() < MIN_FRACT)
2156 {
2157 fraction <<= 1;
2158 exponent--;
2159 }
2160
2161 //
2162 // Check and respond to overflow
2163 //
2164 if (exponent > BIAS)
2165 return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
2166
2167 //
2168 // Check and respond to underflow
2169 //
2170 if (exponent <= -BIAS)
2171 {
2172 if (exponent < -BIAS - FRACT_SIZE)
2173 return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
2174 while (exponent <= -BIAS)
2175 {
2176 sticky |= round;
2177 round = (fraction.LowWord() & 1) != 0;
2178 fraction >>= 1;
2179 exponent++;
2180 }
2181 if (round && (sticky || (fraction.LowWord() & 1)))
2182 fraction++;
2183 exponent = -BIAS;
2184 if (fraction == 0)
2185 return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
2186 }
2187
2188 return IEEEdouble((sign << 31) | ((exponent + BIAS) << FRACT_SIZE_HI)
2189 | (fraction.HighWord() & FRACT_BITS),
2190 fraction.LowWord());
2191 }
2192
SplitInto(BaseLong & fraction) const2193 int IEEEdouble::SplitInto(BaseLong &fraction) const
2194 {
2195 int exponent = Exponent();
2196 fraction = Fraction();
2197
2198 if (exponent == -(int) BIAS)
2199 {
2200 exponent++;
2201 while (fraction.HighWord() < MIN_FRACT)
2202 {
2203 fraction <<= 1;
2204 exponent--;
2205 }
2206 }
2207
2208 return exponent;
2209 }
2210
Ulp() const2211 IEEEdouble IEEEdouble::Ulp() const
2212 {
2213 i4 L;
2214 IEEEdouble d;
2215 d.value.double_value = value.double_value;
2216
2217 L = (i4) d.ExpBits() - FRACT_SIZE * MIN_FRACT;
2218 if (L > 0)
2219 d.setHighAndLowWords((u4) L, 0);
2220 else
2221 {
2222 L = -L >> FRACT_SIZE_HI;
2223 if (L < (i4) FRACT_SIZE_HI)
2224 d.setHighAndLowWords(MIN_FRACT >> (L + 1), 0);
2225 else
2226 {
2227 L -= FRACT_SIZE_HI;
2228 d.setHighAndLowWords(0, L >= 31 ? 1 : 1 << (31 - L));
2229 }
2230 }
2231 return d;
2232 }
2233
Ratio(const BigInt & a,const BigInt & b)2234 IEEEdouble IEEEdouble::Ratio(const BigInt &a, const BigInt &b)
2235 {
2236 IEEEdouble da, db;
2237 int k;
2238 da = a.DoubleValue();
2239 db = b.DoubleValue();
2240 k = b.hi0bits() - a.hi0bits() + 32 * (a.wds - b.wds);
2241 if (k > 0)
2242 da.setHighWord(da.HighWord() + k * MIN_FRACT);
2243 else
2244 db.setHighWord(db.HighWord() - k * MIN_FRACT);
2245 return da / db;
2246 }
2247
operator ==(const IEEEdouble op) const2248 bool IEEEdouble::operator== (const IEEEdouble op) const
2249 {
2250 // FIXME: This could be sped up by inlining
2251 #ifdef HAVE_IEEE754
2252 // TODO: Microsoft VC++ botches this, mixing 12.0 and NaN
2253 return value.double_value == op.value.double_value;
2254 #else
2255 return (IsNaN() || op.IsNaN() ? false
2256 : IsZero() && op.IsZero() ? true
2257 : (BaseLong) *this == (BaseLong) op);
2258 #endif
2259 }
2260
operator !=(const IEEEdouble op) const2261 bool IEEEdouble::operator!= (const IEEEdouble op) const
2262 {
2263 // FIXME: This could be sped up by inlining
2264 #ifdef HAVE_IEEE754
2265 return value.double_value != op.value.double_value;
2266 #else
2267 return !(*this == op);
2268 #endif
2269 }
2270
operator +(const IEEEdouble op) const2271 IEEEdouble IEEEdouble::operator+ (const IEEEdouble op) const
2272 {
2273 #ifdef HAVE_IEEE754
2274 // FIXME: This could be sped up by inlining
2275 return IEEEdouble(value.double_value + op.value.double_value);
2276 #else
2277 if (IsNaN() || op.IsNaN())
2278 return NaN(); // arithmetic on NaNs not allowed
2279
2280 //
2281 // Adding unlike infinities not allowed
2282 // Adding infinities of same sign is infinity of that sign
2283 // Adding finite and infinity produces infinity
2284 //
2285 if (IsInfinite())
2286 return op.IsInfinite() && (Sign() != op.Sign()) ? NaN() : *this;
2287 if (op.IsInfinite())
2288 return op;
2289
2290 //
2291 // Adding zero is easy
2292 //
2293 if (IsZero())
2294 return (op.IsZero()) ? (Sign() != op.Sign()) ? POSITIVE_ZERO()
2295 : *this
2296 : op;
2297 if (op.IsZero())
2298 return *this;
2299
2300 //
2301 // Now for the real work - do manipulations on copies
2302 //
2303 LongInt x, y, round = 0;
2304 int expx, expy, signx, signy;
2305
2306 expx = SplitInto(x);
2307 expy = op.SplitInto(y);
2308
2309 signx = Sign();
2310 signy = op.Sign();
2311
2312 // If the exponents are far enough apart, the answer is easy
2313 if (expx > expy + 54)
2314 return *this;
2315 if (expy > expx + 54)
2316 return op;
2317
2318 //
2319 // Denormalize the fractions, so that the exponents are
2320 // the same and then set the exponent for the result.
2321 // Leave enough space for overflow and LONG_MIN avoidance!
2322 //
2323 if (signx)
2324 x = -x;
2325 if (signy)
2326 y = -y;
2327
2328 x <<= 8;
2329 y <<= 8;
2330
2331 if (expx > expy)
2332 {
2333 round = y << (64 + expy - expx);
2334 y >>= expx - expy;
2335 }
2336 else if (expy > expx)
2337 {
2338 round = x << (64 + expx - expy);
2339 x >>= expy - expx;
2340 expx = expy;
2341 }
2342
2343 //
2344 // Do the arithmetic. The excess magnitude of 64-bit arithmetic means
2345 // overflow is impossible (we only need 1 spare bit!). We ensure that
2346 // pre-alignment avoids any question of LONG_MIN negation problems.
2347 //
2348 x += y;
2349 if (round != 0)
2350 x |= 1;
2351
2352 //
2353 // If the result is negative, then make the fraction positive again
2354 // and remember the sign.
2355 //
2356 if (x < 0)
2357 {
2358 x = -x;
2359 signx = 1;
2360 }
2361 else
2362 signx = 0;
2363
2364 if (x == 0)
2365 return signx ? NEGATIVE_ZERO() : POSITIVE_ZERO();
2366
2367 //
2368 // Time to normalize again! If we need to shift left (the addition was
2369 // effectively a subtraction), then there cannot be any reason to round.
2370 // If the number fits exactly we don't have anything to do either.
2371 //
2372 return Normalize(signx, expx - 8, (ULongInt) x);
2373 #endif
2374 }
2375
operator -() const2376 IEEEdouble IEEEdouble::operator- () const
2377 {
2378 // FIXME: This could be sped up by inlining
2379 #ifdef HAVE_IEEE754
2380 return IEEEdouble(-value.double_value);
2381 #else
2382 if (IsNaN())
2383 return *this;
2384 return IEEEdouble(HighWord() ^ SIGN_BIT, LowWord());
2385 #endif // HAVE_IEEE754
2386 }
2387
operator *(const IEEEdouble op) const2388 IEEEdouble IEEEdouble::operator* (const IEEEdouble op) const
2389 {
2390 #ifdef HAVE_IEEE754
2391 return IEEEdouble(value.double_value * op.value.double_value);
2392 #else
2393 if (IsNaN() || op.IsNaN())
2394 return NaN(); // arithmetic on NaNs not allowed
2395
2396 int sign = Sign() ^ op.Sign();
2397
2398 //
2399 // If either operand is zero or infinite, then the answer is easy.
2400 //
2401 if (IsZero())
2402 return op.IsInfinite() ? NaN()
2403 : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
2404
2405 if (op.IsZero())
2406 return IsInfinite() ? NaN()
2407 : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
2408
2409 if (IsInfinite() || op.IsInfinite())
2410 return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
2411
2412 //
2413 // Now for the real work - do manipulations on copies
2414 //
2415 ULongInt x, y, z, w, pr1, pr2;
2416 u4 round;
2417 int exponent;
2418
2419 exponent = SplitInto(x) + op.SplitInto(y);
2420
2421 //
2422 // The numbers to be multiplied are 53 bits in length (stored in 64 bit
2423 // integers). Split them in 32 bit parts, then shift result into place.
2424 // Fold the low order bits into the lsb for rounding purposes.
2425 //
2426 z = ULongInt(x.LowWord());
2427 w = ULongInt(y.LowWord());
2428 x >>= 32;
2429 y >>= 32;
2430
2431 pr1 = z * w; // low*low
2432 pr2 = z * y + x * w + pr1.HighWord(); // low*high + high*low + adjusted pr1
2433 round = pr1.LowWord();
2434 x = (x * y) << 20; // high*high, adjusted
2435
2436 round |= pr2.LowWord() & 0xFFF;
2437 pr2 >>= 12;
2438
2439 x += pr2 | (round ? 1 : 0);
2440
2441 return Normalize(sign, exponent - 8, x);
2442 #endif // HAVE_IEEE754
2443 }
2444
operator /(const IEEEdouble op) const2445 IEEEdouble IEEEdouble::operator/ (const IEEEdouble op) const
2446 {
2447 #ifdef HAVE_IEEE754
2448 return op.IsZero() ? ((IsNaN() || IsZero()) ? NaN()
2449 : ((IsPositive() ^ op.IsPositive())
2450 ? NEGATIVE_INFINITY()
2451 : POSITIVE_INFINITY()))
2452 : IEEEdouble(value.double_value / op.value.double_value);
2453 #else // HAVE_IEEE754
2454 if (IsNaN() || op.IsNaN())
2455 return NaN(); // arithmetic on NaNs not allowed
2456
2457 int sign = Sign() ^ op.Sign();
2458
2459 //
2460 // Infinities and zeroes are special.
2461 //
2462
2463 if (IsInfinite())
2464 return (op.IsInfinite() ? NaN()
2465 : sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY());
2466
2467 if (op.IsInfinite())
2468 return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
2469
2470 if (IsZero())
2471 return op.IsZero() ? NaN()
2472 : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
2473
2474 if (op.IsZero())
2475 return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
2476
2477 //
2478 // Now for the real work - do manipulations on copies
2479 //
2480 ULongInt x, y;
2481 int exponent;
2482
2483 exponent = SplitInto(x) - op.SplitInto(y);
2484
2485 ULongInt mask (0x80000000, 0x00000000),
2486 result(0x00000000, 0x00000000);
2487
2488 // Because both values are normalised, a single shift guarantees results.
2489
2490 if (x < y)
2491 {
2492 x <<= 1;
2493 exponent--;
2494 }
2495
2496 //
2497 // If the numerator is larger, then it is divisible.
2498 // Reflect this in the result, and do the subtraction.
2499 // Magnify the numerator again and reduce the mask.
2500 //
2501 while (mask != 0)
2502 {
2503 if (x >= y)
2504 {
2505 result += mask;
2506 x -= y;
2507 if (x == 0)
2508 break;
2509 }
2510 x <<= 1;
2511 mask >>= 1;
2512 }
2513
2514 return Normalize(sign, exponent - 11, result);
2515 #endif // HAVE_IEEE754
2516 }
2517
operator %(const IEEEdouble op) const2518 IEEEdouble IEEEdouble::operator% (const IEEEdouble op) const
2519 {
2520 #ifdef HAVE_IEEE754
2521 return IEEEdouble((op.IsZero() ? NaN().value.double_value
2522 : fmod(value.double_value, op.value.double_value)));
2523 #else // HAVE_IEEE754
2524 if (IsNaN() || op.IsNaN())
2525 return NaN(); // arithmetic on NaNs not allowed
2526
2527 //
2528 // Infinities and zeroes are special.
2529 //
2530
2531 if (IsInfinite() || op.IsZero())
2532 return NaN();
2533
2534 if (IsZero() || op.IsInfinite())
2535 return *this;
2536
2537 //
2538 // Now for the real work - do manipulations on copies
2539 // This algorithm is from fdlibm.c - see above notice
2540 //
2541
2542 int sign = Sign();
2543 ULongInt x, y;
2544 int expy, exponent;
2545 LongInt z;
2546
2547 expy = op.SplitInto(y);
2548 exponent = SplitInto(x) - expy;
2549
2550 if (exponent < 0 || (exponent == 0 && x < y))
2551 return *this;
2552
2553 while (exponent--)
2554 {
2555 z = x - y;
2556 if (z < 0)
2557 x <<= 1;
2558 else if (z == 0)
2559 return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
2560 else
2561 x = z + z;
2562 }
2563 z = x - y;
2564 if (z >= 0)
2565 x = z;
2566
2567 if (x == 0)
2568 return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
2569 return Normalize(sign, expy, x);
2570 #endif // HAVE_IEEE754
2571 }
2572
operator <(const IEEEdouble op) const2573 bool IEEEdouble::operator< (const IEEEdouble op) const
2574 {
2575 // FIXME: This could be sped up by inlining
2576 #ifdef HAVE_IEEE754
2577 return (value.double_value < op.value.double_value);
2578 #else
2579 if (IsNaN() || op.IsNaN())
2580 return false; // NaNs are unordered
2581 if (IsZero() && op.IsZero())
2582 return false;
2583 // Exploit fact that all other IEEE floating point numbers sort like
2584 // ints after worrying about sign.
2585 u4 a = HighWord() & ABS_BITS, b = op.HighWord() & ABS_BITS;
2586 if (IsNegative())
2587 return op.IsPositive() ||
2588 (a > b || (a == b && LowWord() > op.LowWord()));
2589 return op.IsPositive() &&
2590 (a < b || (a == b && LowWord() < op.LowWord()));
2591 #endif
2592 }
2593
operator <=(const IEEEdouble op) const2594 bool IEEEdouble::operator<= (const IEEEdouble op) const
2595 {
2596 // FIXME: This could be sped up by inlining
2597 #ifdef HAVE_IEEE754
2598 return (value.double_value <= op.value.double_value);
2599 #else
2600 return *this < op || *this == op;
2601 #endif
2602 }
2603
operator >(const IEEEdouble op) const2604 bool IEEEdouble::operator> (const IEEEdouble op) const
2605 {
2606 // FIXME: This could be sped up by inlining
2607 #ifdef HAVE_IEEE754
2608 return (value.double_value > op.value.double_value);
2609 #else
2610 if (IsNaN() || op.IsNaN())
2611 return false; // NaNs are unordered.
2612 if (IsZero() && op.IsZero())
2613 return false;
2614 // Exploit fact that all other IEEE floating point numbers sort like
2615 // ints after worrying about sign.
2616 u4 a = HighWord() & ABS_BITS, b = op.HighWord() & ABS_BITS;
2617 if (IsPositive())
2618 return op.IsNegative() ||
2619 (a > b || (a == b && LowWord() > op.LowWord()));
2620 return op.IsNegative() &&
2621 (a < b || (a == b && LowWord() < op.LowWord()));
2622 #endif
2623 }
2624
operator >=(const IEEEdouble op) const2625 bool IEEEdouble::operator>= (const IEEEdouble op) const
2626 {
2627 // FIXME: This could be sped up by inlining
2628 #ifdef HAVE_IEEE754
2629 return (value.double_value >= op.value.double_value);
2630 #else
2631 return *this > op || *this == op;
2632 #endif
2633 }
2634
2635
2636 #ifndef HAVE_MEMBER_CONSTANTS
2637 u4 BigInt::fives[] = {
2638 1, 5, 25, 125, 625, 3125, 15625, 78125
2639 };
2640 BigInt *BigInt::bigfives[] = {
2641 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL
2642 };
2643 #else // HAVE_MEMBER_CONSTANTS
2644 const u4 BigInt::fives[] = {
2645 1, 5, 25, 125, 625, 3125, 15625, 78125
2646 };
2647 const BigInt *BigInt::bigfives[] = {
2648 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL
2649 };
2650 #endif // HAVE_MEMBER_CONSTANTS
2651
BigInt(const IEEEfloat & f,int & e,int & bits)2652 BigInt::BigInt(const IEEEfloat &f, int &e, int &bits) : data(NULL)
2653 {
2654 int fe, k;
2655 u4 z;
2656 resize(0);
2657 z = f.Fraction();
2658 fe = f.Exponent();
2659 k = lo0bits(z);
2660 data[0] = z;
2661 wds = 1;
2662 if (fe + IEEEfloat::Bias())
2663 {
2664 e = k + fe - IEEEfloat::FractSize();
2665 bits = IEEEfloat::FractSize() - k + 1;
2666 }
2667 else
2668 {
2669 e = k + fe - IEEEfloat::FractSize() + 1;
2670 bits = 32 - hi0bits(z);
2671 }
2672 }
2673
BigInt(const IEEEdouble & d,int & e,int & bits)2674 BigInt::BigInt(const IEEEdouble &d, int &e, int &bits) : data(NULL)
2675 {
2676 int de, k;
2677 LongInt x;
2678 u4 y, z;
2679
2680 resize(1);
2681 x = d.Fraction();
2682 z = x.HighWord();
2683 de = d.Exponent();
2684 if ((y = x.LowWord()) != 0)
2685 {
2686 if ((k = lo0bits(y)) != 0)
2687 {
2688 data[0] = y | z << (32 - k);
2689 z >>= k;
2690 }
2691 else
2692 data[0] = y;
2693 wds = (data[1] = z) ? 2 : 1;
2694 }
2695 else
2696 {
2697 k = lo0bits(z) + 32;
2698 data[0] = z;
2699 wds = 1;
2700 }
2701 if (de + IEEEdouble::Bias())
2702 {
2703 e = k + de - IEEEdouble::FractSize();
2704 bits = IEEEdouble::FractSize() - k + 1;
2705 }
2706 else
2707 {
2708 e = k + de - IEEEdouble::FractSize() + 1;
2709 bits = 32 * wds - hi0bits(data[wds - 1]);
2710 }
2711 }
2712
BigInt(const char * s,int nd0,int nd,u4 start,int startsize)2713 BigInt::BigInt(const char *s, int nd0, int nd, u4 start, int startsize) :
2714 data(NULL)
2715 {
2716 int i, k;
2717 u4 x, y;
2718 x = (nd + 8) / 9;
2719 for (k = 0, y = 1; x > y; y <<= 1, k++);
2720
2721 resize(k);
2722 data[0] = start;
2723 wds = 1;
2724 i = startsize;
2725 if (startsize < nd0)
2726 {
2727 s += startsize;
2728 do
2729 multadd(10, *s++ - U_0);
2730 while (++i < nd0);
2731 s++;
2732 }
2733 else
2734 s += startsize + 1;
2735
2736 for ( ; i < nd; i++)
2737 multadd(10, *s++ - U_0);
2738 }
2739
operator =(const BigInt & b)2740 BigInt &BigInt::operator =(const BigInt &b)
2741 {
2742 if (this != &b)
2743 {
2744 k = b.k;
2745 maxwds = b.maxwds;
2746 neg = b.neg;
2747 wds = b.wds;
2748 delete data;
2749 data = new u4[maxwds];
2750 memcpy(data, b.data, wds * sizeof(u4));
2751 }
2752 return *this;
2753 }
2754
hi0bits(u4 x)2755 int BigInt::hi0bits(u4 x)
2756 {
2757 int k = 0;
2758 if (! (x & 0xffff0000))
2759 {
2760 k = 16;
2761 x <<= 16;
2762 }
2763 if (! (x & 0xff000000))
2764 {
2765 k += 8;
2766 x <<= 8;
2767 }
2768 if (! (x & 0xf0000000))
2769 {
2770 k += 4;
2771 x <<= 4;
2772 }
2773 if (! (x & 0xc0000000))
2774 {
2775 k += 2;
2776 x <<= 2;
2777 }
2778 if (! (x & 0x80000000))
2779 {
2780 k++;
2781 if (! (x & 0x40000000))
2782 return 32;
2783 }
2784 return k;
2785 }
2786
lo0bits(u4 & y)2787 int BigInt::lo0bits(u4 &y)
2788 {
2789 int k;
2790 if (y & 7)
2791 {
2792 if (y & 1)
2793 return 0;
2794 if (y & 2)
2795 {
2796 y >>= 1;
2797 return 1;
2798 }
2799 y >>= 2;
2800 return 2;
2801 }
2802 k = 0;
2803 if (! (y & 0xffff))
2804 {
2805 k = 16;
2806 y >>= 16;
2807 }
2808 if (! (y & 0xff))
2809 {
2810 k += 8;
2811 y >>= 8;
2812 }
2813 if (! (y & 0xf))
2814 {
2815 k += 4;
2816 y >>= 4;
2817 }
2818 if (! (y & 0x3))
2819 {
2820 k += 2;
2821 y >>= 2;
2822 }
2823 if (! (y & 1))
2824 {
2825 k++;
2826 y >>= 1;
2827 if (! y)
2828 return 32;
2829 }
2830 return k;
2831 }
2832
operator +(const unsigned op) const2833 BigInt BigInt::operator +(const unsigned op) const
2834 {
2835 int i = 0; // counter
2836 u4 carry = op; // carry between words
2837 ULongInt sum; // sum
2838 BigInt result(*this);
2839 u4 *x = result.data; // access to data
2840
2841 do
2842 {
2843 sum = ULongInt(*x) + carry;
2844 carry = sum.HighWord();
2845 *x++ = sum.LowWord();
2846 } while (carry && ++i < wds);
2847 if (carry && i == wds)
2848 {
2849 if (wds == maxwds)
2850 {
2851 result.maxwds = 1 << (++result.k);
2852 x = new u4[result.maxwds];
2853 memcpy(x, result.data, wds * sizeof(u4));
2854 delete result.data;
2855 result.data = x;
2856 }
2857 result.data[result.wds++] = carry;
2858 }
2859 return result;
2860 }
2861
operator -(const BigInt & op) const2862 BigInt BigInt::operator -(const BigInt &op) const
2863 {
2864 const BigInt *a = this, *b = &op;
2865 BigInt zero(0);
2866 int i, wa, wb;
2867 u4 *xa, *xae, *xb, *xbe, *xc;
2868 u4 borrow;
2869 ULongInt y;
2870
2871 i = a -> compareTo(op);
2872 if (! i)
2873 return zero;
2874 if (i < 0)
2875 {
2876 const BigInt *tmp = a;
2877 a = b;
2878 b = tmp;
2879 i = 1;
2880 }
2881 else
2882 i = 0;
2883
2884 BigInt c(0);
2885 c.resize(a -> k);
2886 c.neg = i != 0;
2887 wa = a -> wds;
2888 xa = a -> data;
2889 xae = xa + wa;
2890 wb = b -> wds;
2891 xb = b -> data;
2892 xbe = xb + wb;
2893 xc = c.data;
2894 borrow = 0;
2895 do
2896 {
2897 y = ULongInt(*xa++) - *xb++ - borrow;
2898 borrow = y.HighWord() & 1;
2899 *xc++ = y.LowWord();
2900 } while (xb < xbe);
2901 while (xa < xae)
2902 {
2903 y = ULongInt(*xa++) - borrow;
2904 borrow = y.HighWord() & 1;
2905 *xc++ = y.LowWord();
2906 }
2907 while (!*--xc)
2908 wa--;
2909 c.wds = wa;
2910 return c;
2911 }
2912
operator *(unsigned op) const2913 BigInt BigInt::operator *(unsigned op) const
2914 {
2915 int i = 0; // counter
2916 u4 carry = 0; // carry between words
2917 ULongInt product; // product
2918 BigInt result(*this);
2919 u4 *x = result.data; // access to data
2920 ULongInt factor = (u4) op; // avoid creating object multiple times
2921
2922 do
2923 {
2924 product = ULongInt(*x) * factor + carry;
2925 carry = product.HighWord();
2926 *x++ = product.LowWord();
2927 } while (++i < wds);
2928 if (carry)
2929 {
2930 if (wds == maxwds)
2931 {
2932 result.maxwds = 1 << (++result.k);
2933 x = new u4[result.maxwds];
2934 memcpy(x, result.data, wds * sizeof(u4));
2935 delete result.data;
2936 result.data = x;
2937 }
2938 result.data[result.wds++] = carry;
2939 }
2940 return result;
2941 }
2942
operator *(const BigInt & op) const2943 BigInt BigInt::operator *(const BigInt &op) const
2944 {
2945 const BigInt *a = this, *b = &op;
2946 int k; // c -> k
2947 int wa, wb, wc; // wds in each of a, b, c
2948 u4 *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
2949 u4 y, carry;
2950 ULongInt z;
2951
2952 if (a -> wds < b -> wds)
2953 {
2954 const BigInt *tmp = a;
2955 a = b;
2956 b = tmp;
2957 }
2958 k = a -> k;
2959 wa = a -> wds;
2960 wb = b -> wds;
2961 wc = wa + wb;
2962 if (wc > a -> maxwds)
2963 k++;
2964 BigInt c(0);
2965 c.resize(k);
2966 c.neg = a -> neg ^ b -> neg;
2967 for (x = c.data, xa = x + wc; x < xa; x++)
2968 *x = 0;
2969 xa = a -> data;
2970 xae = xa + wa;
2971 xb = b -> data;
2972 xbe = xb + wb;
2973 xc0 = c.data;
2974 for ( ; xb < xbe; xc0++)
2975 {
2976 if ((y = *xb++) != 0)
2977 {
2978 x = xa;
2979 xc = xc0;
2980 carry = 0;
2981 do
2982 {
2983 z = ULongInt(*x++) * y + *xc + carry;
2984 carry = z.HighWord();
2985 *xc++ = z.LowWord();
2986 } while (x < xae);
2987 *xc = carry;
2988 }
2989 }
2990 for (xc0 = c.data, xc = xc0 + wc; wc > 0 && !*--xc; --wc);
2991
2992 c.wds = wc;
2993 return c;
2994 }
2995
operator <<(unsigned op) const2996 BigInt BigInt::operator <<(unsigned op) const
2997 {
2998 int i, k1, n, n1;
2999 u4 *x, *x1, *xe, z;
3000 n = op >> 5;
3001 k1 = k;
3002 n1 = n + wds + 1;
3003 for (i = maxwds; n1 > i; i <<= 1)
3004 k1++;
3005 BigInt result(*this);
3006 result.maxwds = 1 << k1;
3007 result.k = k1;
3008 delete result.data;
3009 result.data = x1 = new u4[result.maxwds];
3010 for (i = 0; i < n; i++)
3011 *x1++ = 0;
3012 x = data;
3013 xe = x + wds;
3014 if (op &= 0x1f)
3015 {
3016 k1 = 32 - op;
3017 z = 0;
3018 do
3019 {
3020 *x1++ = *x << op | z;
3021 z = *x++ >> k1;
3022 } while (x < xe);
3023 if ((*x1 = z) != 0)
3024 ++n1;
3025 }
3026 else
3027 do
3028 *x1++ = *x++;
3029 while (x < xe);
3030 result.wds = n1 - 1;
3031 return result;
3032 }
3033
multadd(unsigned m,unsigned a)3034 BigInt &BigInt::multadd(unsigned m, unsigned a)
3035 {
3036 int i = 0; // counter
3037 u4 *x = data; // access to data
3038 u4 carry = a; // carry between words
3039 ULongInt product; // product
3040 ULongInt factor = (u4) m; // avoid creating object multiple times
3041
3042 do
3043 {
3044 product = ULongInt(*x) * factor + carry;
3045 carry = product.HighWord();
3046 *x++ = product.LowWord();
3047 } while (++i < wds);
3048 if (carry)
3049 {
3050 if (wds == maxwds)
3051 {
3052 maxwds = 1 << (++k);
3053 x = new u4[maxwds];
3054 memcpy(x, data, wds * sizeof(u4));
3055 delete data;
3056 data = x;
3057 }
3058 data[wds++] = carry;
3059 }
3060 return *this;
3061 }
3062
pow5mult(unsigned k)3063 BigInt &BigInt::pow5mult(unsigned k)
3064 {
3065 const BigInt *p5;
3066 int i;
3067 assert(k < 0x800);
3068 if ((i = k & 0x7) != 0)
3069 *this *= fives[i];
3070 if (! (k >>= 3))
3071 return *this;
3072 if (! (p5 = bigfives[i = 0]))
3073 p5 = bigfives[i] = new BigInt(390625);
3074 while (true)
3075 {
3076 if (k & 1)
3077 *this *= *p5;
3078 if (! (k >>= 1))
3079 break;
3080 if (! (p5 = bigfives[++i]))
3081 p5 = bigfives[i] = new BigInt(*bigfives[i-1] * *bigfives[i-1]);
3082 }
3083 return *this;
3084 }
3085
compareTo(const BigInt & b) const3086 int BigInt::compareTo(const BigInt &b) const
3087 {
3088 u4 *xa, *xa0, *xb, *xb0;
3089 int i, j;
3090 i = wds;
3091 j = b.wds;
3092 if (i -= j)
3093 return i;
3094 xa0 = data;
3095 xa = xa0 + j;
3096 xb0 = b.data;
3097 xb = xb0 + j;
3098 while (xa > xa0)
3099 if (*--xa != *--xb)
3100 return *xa < *xb ? -1 : 1;
3101 return 0;
3102 }
3103
quorem(const BigInt & S)3104 int BigInt::quorem(const BigInt &S)
3105 {
3106 int n;
3107 u4 *bx, *bxe, q, *sx, *sxe;
3108 u4 borrow, carry;
3109 ULongInt y, ys;
3110 n = S.wds;
3111 if (wds < n)
3112 return 0;
3113 sx = S.data;
3114 sxe = sx + --n;
3115 bx = data;
3116 bxe = bx + n;
3117 q = *bxe / (*sxe + 1);
3118 if (q)
3119 {
3120 borrow = 0;
3121 carry = 0;
3122 do
3123 {
3124 ys = ULongInt(*sx++) * q + carry;
3125 carry = ys.HighWord();
3126 y = ULongInt(*bx) - ys.LowWord() - borrow;
3127 borrow = y.HighWord() & 1;
3128 *bx++ = y.LowWord();
3129 } while (sx <= sxe);
3130 if (!*bxe)
3131 {
3132 bx = data;
3133 while (--bxe > bx && !*bxe)
3134 --n;
3135 wds = n;
3136 }
3137 }
3138 if (compareTo(S) >= 0)
3139 {
3140 q++;
3141 borrow = 0;
3142 carry = 0;
3143 bx = data;
3144 sx = S.data;
3145 do
3146 {
3147 ys = ULongInt(*sx++) + carry;
3148 carry = ys.HighWord();
3149 y = ULongInt(*bx) - ys.LowWord() - borrow;
3150 borrow = y.HighWord() & 1;
3151 *bx++ = y.LowWord();
3152 } while (sx <= sxe);
3153 bx = data;
3154 bxe = bx + n;
3155 if (!*bxe)
3156 {
3157 while (--bxe > bx && !*bxe)
3158 --n;
3159 wds = n;
3160 }
3161 }
3162 return q;
3163 }
3164
FloatValue() const3165 IEEEfloat BigInt::FloatValue() const
3166 {
3167 u4 *xa, y, z;
3168 int k;
3169 xa = data + wds;
3170 y = *--xa;
3171 k = hi0bits(y);
3172 if (k < 8)
3173 return IEEEfloat(0x3f800000 | y >> (8 - k));
3174 z = xa > data ? *--xa : 0;
3175 if (k -= 8)
3176 return IEEEfloat(0x3f800000 | y << k | z >> (32 - k));
3177 else
3178 return IEEEfloat(0x3f800000 | y);
3179 }
3180
DoubleValue() const3181 IEEEdouble BigInt::DoubleValue() const
3182 {
3183 u4 *xa, w, y, z, hi, lo;
3184 int k;
3185 xa = data + wds;
3186 y = *--xa;
3187 k = hi0bits(y);
3188 if (k < 11)
3189 {
3190 hi = 0x3ff00000 | y >> (11 - k);
3191 w = xa > data ? *--xa : 0;
3192 lo = y << (32 - 11 + k) | w >> (11 - k);
3193 return IEEEdouble(hi, lo);
3194 }
3195 z = xa > data ? *--xa : 0;
3196 if (k -= 11)
3197 {
3198 hi = 0x3ff00000 | y << k | z >> (32 - k);
3199 y = xa > data ? *--xa : 0;
3200 lo = z << k | y >> (32 - k);
3201 }
3202 else
3203 {
3204 hi = 0x3ff00000 | y;
3205 lo = z;
3206 }
3207 return IEEEdouble(hi, lo);
3208 }
3209
3210
3211
3212 #ifdef HAVE_JIKES_NAMESPACE
3213 } // Close namespace Jikes block
3214 #endif
3215
3216