home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 20 / AACD20.BIN / AACD / Programming / Jikes / Source / src / double.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  2001-02-24  |  33.5 KB  |  1,357 lines
  1. // $Id: double.cpp,v 1.17 2001/02/14 21:25:11 mdejong 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://www.ibm.com/research/jikes.
  6. // Copyright (C) 1996, 1998, International Business Machines Corporation
  7. // and others.  All Rights Reserved.
  8. // You must accept the terms of that agreement to use this software.
  9. //
  10. //
  11. // NOTE: The IEEE 754 emulation code in double.h and double.cpp within
  12. // Jikes are adapted from code written by Alan M. Webb of IBM's Hursley
  13. // lab in porting the Sun JDK to System/390.
  14. //
  15. // In addition, the code for emulating the remainder operator, %, is
  16. // adapted from e_fmod.c, part of fdlibm, the Freely Distributable Math
  17. // Library mentioned in the documentation of java.lang.StrictMath.  The
  18. // original library is available at http://netlib2.cs.utk.edu/fdlibm.
  19. //
  20. // The code from fdlibm is copyrighted, as follows:
  21. // ====================================================
  22. // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  23. //
  24. // Developed at SunSoft, a Sun Microsystems, Inc. business.
  25. // Permission to use, copy, modify, and distribute this
  26. // software is freely granted, provided that this notice 
  27. // is preserved.
  28. // ====================================================
  29. //
  30. //
  31.  
  32. #include "double.h"
  33. #include "long.h"
  34.  
  35. #ifdef    HAVE_JIKES_NAMESPACE
  36. namespace Jikes {    // Open namespace Jikes block
  37. #endif
  38.  
  39. #ifndef HAVE_MEMBER_CONSTANTS
  40. // VC++ can't cope with constant class members
  41. u4 IEEEfloat::MAX_FRACT   = 0x01000000;
  42. u4 IEEEfloat::MAX_FRACT2  = 0x00800000;
  43. u4 IEEEfloat::MIN_INT_F   = 0xCF000000;
  44. i4 IEEEfloat::MIN_INT     = 0x80000000;
  45. i4 IEEEfloat::MAX_INT     = 0x7FFFFFFF;
  46.  
  47. u4 IEEEdouble::MAX_FRACT  = 0x00200000;
  48. u4 IEEEdouble::MAX_FRACT2 = 0x00100000;
  49. i4 IEEEdouble::MIN_INT    = 0x80000000;
  50. i4 IEEEdouble::MAX_INT    = 0x7FFFFFFF;
  51. #else
  52. // gcc bug 1877 can cause linker errors if
  53. // the following external decls are not used.
  54. const u4 IEEEfloat::MAX_FRACT;
  55. const u4 IEEEfloat::MAX_FRACT2;
  56. const u4 IEEEfloat::MIN_INT_F;
  57. const i4 IEEEfloat::MIN_INT;
  58. const i4 IEEEfloat::MAX_INT;
  59.  
  60. const u4 IEEEdouble::MAX_FRACT;
  61. const u4 IEEEdouble::MAX_FRACT2;
  62. const i4 IEEEdouble::MIN_INT;
  63. const i4 IEEEdouble::MAX_INT;
  64. #endif
  65.  
  66.  
  67. IEEEfloat::IEEEfloat(i4 a)
  68. {
  69. #ifdef HAVE_IEEE754
  70.     value.float_value = (float) a;
  71. #else
  72.     int    sign = 0;
  73.  
  74.     if (a < 0)
  75.     {
  76.         a  = -a; // even works for MIN_INT!
  77.         sign = 1;
  78.     }
  79.     if (a == 0)
  80.         *this = POSITIVE_ZERO();
  81.     else
  82.         *this = Normalize(sign, FRACT_BITS, (u4) a);
  83. #endif // HAVE_IEEE754
  84. }
  85.  
  86. IEEEfloat::IEEEfloat(LongInt a)
  87. {
  88. #ifdef HAVE_IEEE754
  89. # ifdef HAVE_UNSIGNED_LONG_LONG
  90.     value.float_value = (float) a.Words();
  91. # else
  92.     value.float_value = ((float) a.HighWord() * (float) 0x40000000 * 4.0f) +
  93.         (float) a.LowWord();
  94. # endif // HAVE_UNSIGNED_LONG_LONG
  95. #else
  96.     //
  97.     // Unfortunately, we cannot recycle the LongInt.DoubleValue() method, since
  98.     // in rare cases the double rounding puts us off by one bit.
  99.     //
  100.     int    sign = 0;
  101.  
  102.     if (a < 0)
  103.     {
  104.         a  = -a;
  105.         sign = 1;
  106.         if (a < 0) // special case MIN_LONG
  107.         {
  108.             value.word = 0xDF000000;
  109.             return;
  110.         }
  111.     }
  112.     if (a == 0)
  113.         *this = POSITIVE_ZERO();
  114.     else if (a.HighWord() == 0)
  115.         *this = Normalize(sign, FRACT_BITS, a.LowWord());
  116.     else
  117.     {
  118.         int exponent = FRACT_BITS, round = 0;
  119.         while (a.HighWord())
  120.         {
  121.             round |= (a.LowWord() & 0x000000ff) ? 1 : 0;
  122.             a >>= 8;
  123.             exponent += 8;
  124.         }
  125.         *this = Normalize(sign, exponent, a.LowWord() | round);
  126.     }
  127. #endif // HAVE_IEEE754
  128. }
  129.  
  130. IEEEfloat::IEEEfloat(IEEEdouble d)
  131. {
  132. #ifdef HAVE_IEEE754
  133.     value.float_value = (float) d.DoubleView();
  134. #else
  135.     // Either true zero, denormalized, or too small
  136.     if (d.Exponent() < -BIAS - 30)
  137.         *this = (d.IsPositive() ? POSITIVE_ZERO() : NEGATIVE_ZERO());
  138.     else
  139.     {
  140.         if (d.IsPositiveInfinity())
  141.             *this = POSITIVE_INFINITY();
  142.         else if (d.IsNegativeInfinity())
  143.             *this = NEGATIVE_INFINITY();
  144.         else if (d.IsNaN())
  145.             *this = NaN();
  146.         else 
  147.     {
  148.             //
  149.             // A regular, normalized number - do work on the parts
  150.             // Shift to 26th position, add implicit msb, rounding bits
  151.             //
  152.             LongInt fract = d.Fraction() << 5;
  153.         u4 fraction = fract.HighWord() | ((fract.LowWord()) ? 0x02000001
  154.                                                                 : 0x02000000);
  155.  
  156.         *this = Normalize(d.Sign(), d.Exponent() - 2, fraction);
  157.     }
  158.     }
  159. #endif // HAVE_IEEE754
  160. }
  161.  
  162. IEEEfloat::IEEEfloat(char *str, bool check_invalid)
  163. {
  164.     //
  165.     // This assumes that name already meets the format of JLS 3.10.2 (ie. no
  166.     // extra whitespace or invalid characters)
  167.     //
  168.     // TODO: This conversion is a temporary patch. Need volunteer to implement
  169.     //       Clinger algorithm from PLDI 1990.
  170.     //
  171.     value.float_value = (float) atof(str);
  172.  
  173.     //
  174.     // When parsing a literal in Java, a number that rounds to infinity or
  175.     // zero is invalid.  A true check_invalid sets any invalid number to NaN,
  176.     // to make the upstream processing easier.  Leave check_invalid false to
  177.     // allow parsing infinity or rounded zero from a string.
  178.     //
  179.     if (check_invalid)
  180.     {
  181.         if (IsInfinite())
  182.             *this = NaN();
  183.         if (IsZero())
  184.         {
  185.             for (char *p = str; *p; p++) {
  186.                 switch (*p) {
  187.                 case '-': case '+': case '0': case '.':
  188.                     break; // keep checking
  189.                 case 'e': case 'E': case 'd': case 'D': case 'f': case 'F':
  190.                     return; // got this far, str is 0
  191.                 default:
  192.                     *this = NaN(); // encountered non-zero digit. oops.
  193.                     return;
  194.                 }
  195.             }
  196.         }
  197.     }
  198. }
  199.  
  200. i4 IEEEfloat::IntValue()
  201. {
  202.     if (IsNaN())
  203.         return 0;
  204.             
  205.     int    sign = Sign(),
  206.         exponent = Exponent();
  207.  
  208.     if (IsInfinite())
  209.         return sign ? MIN_INT : MAX_INT;
  210.  
  211.     // This covers true zero and denorms.
  212.     if (exponent < 0)
  213.         return 0;
  214.  
  215.     i4 result = Fraction();
  216.  
  217.     if (exponent > FRACT_BITS)
  218.     result <<= (exponent - FRACT_BITS);
  219.     else if (exponent < FRACT_BITS)
  220.     result >>= (FRACT_BITS - exponent);
  221.  
  222.     return sign ? -result : result;
  223. }
  224.  
  225. LongInt IEEEfloat::LongValue()
  226. {
  227.     if (IsNaN())
  228.         return LongInt(0);
  229.             
  230.     int    sign = Sign(),
  231.         exponent = Exponent();
  232.  
  233.     if (IsInfinite())
  234.         return sign ? LongInt::MIN_LONG() : LongInt::MAX_LONG();
  235.  
  236.     // This covers true zero and denorms.
  237.     if (exponent < 0)
  238.         return LongInt(0);
  239.  
  240.     LongInt result(Fraction());
  241.  
  242.     if (exponent > FRACT_BITS)
  243.     result <<= (exponent - FRACT_BITS);
  244.     else if (exponent < FRACT_BITS)
  245.     result >>= (FRACT_BITS - exponent);
  246.  
  247.     return sign ? (LongInt) -result : result;
  248. }
  249.  
  250. IEEEfloat IEEEfloat::Normalize(int sign, int exponent, u4 fraction)
  251. {
  252.     bool round = false, sticky = false;
  253.  
  254.     assert(fraction);
  255.  
  256.     //
  257.     // Normalize right. MAX_FRACT is (FLT_MAX<<1)-FLT_MAX.
  258.     // So we need to shift on equal.
  259.     //
  260.     if (fraction >= MAX_FRACT)
  261.     {
  262.     while (fraction >= MAX_FRACT)
  263.     {
  264.         sticky |= round;
  265.         round = (fraction & 0x00000001) != 0;
  266.         fraction >>= 1;
  267.         exponent++;
  268.     }
  269.     if (round && (sticky || (fraction & 0x00000001)) && exponent > -BIAS)
  270.             //
  271.             // Capture any overflow caused by rounding. No other checks are
  272.             // required because if overflow occurred, the the low order bit
  273.             // was guaranteed to be zero. Do not round denorms yet.
  274.             //
  275.         if (++fraction >= MAX_FRACT)
  276.         {
  277.         fraction >>= 1;
  278.         exponent++;
  279.         }
  280.     }
  281.  
  282.     // Normalize left. MAX_FRACT2 is MAX_FRACT >> 1.
  283.     else
  284.     while (fraction < MAX_FRACT2)
  285.     {
  286.         fraction <<= 1;
  287.         exponent--;
  288.     }
  289.  
  290.     //
  291.     // Check and respond to overflow
  292.     //
  293.     if (exponent > BIAS)
  294.         return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
  295.  
  296.     //
  297.     // Check and respond to underflow 
  298.     //
  299.     if (exponent <= -BIAS)
  300.     {
  301.     while (exponent <= -BIAS)
  302.     {
  303.         sticky |= round;
  304.         round = (fraction & 0x00000001) != 0;
  305.         fraction >>= 1;
  306.         exponent++;
  307.     }
  308.  
  309.     if (round && (sticky || (fraction & 0x00000001)))
  310.         fraction++;
  311.  
  312.     exponent = -BIAS;
  313.  
  314.     if (fraction == 0)
  315.             return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  316.     }
  317.  
  318.     fraction &= 0x007FFFFF;
  319.     fraction |= ((exponent + BIAS) << FRACT_BITS);
  320.     if (sign)
  321.         fraction |= 0x80000000;
  322.     return IEEEfloat(fraction);
  323. }
  324.  
  325. int IEEEfloat::SplitInto(u4 &fraction)
  326. {
  327.     int exponent = Exponent();
  328.     fraction = Fraction();
  329.  
  330.     if (exponent == -BIAS)
  331.     {
  332.     exponent++;
  333.     while (fraction < MAX_FRACT2)
  334.     {
  335.         fraction <<= 1;
  336.         exponent--;
  337.     }
  338.     }
  339.  
  340.     return exponent;
  341. }
  342.  
  343. bool IEEEfloat::operator== (IEEEfloat op)
  344. {
  345.     // FIXME: This could be sped up by inlining
  346. #ifdef HAVE_IEEE754
  347.     return value.float_value == op.value.float_value;
  348. #else
  349.     return IsNaN() || op.IsNaN() ? false
  350.                                  : IsZero() && op.IsZero() ? true
  351.                                                            : value.word == op.value.word;
  352. #endif
  353. }
  354.  
  355. bool IEEEfloat::operator!= (IEEEfloat op)
  356. {
  357.     // FIXME: This could be sped up by inlining
  358. #ifdef HAVE_IEEE754
  359.     return value.float_value != op.value.float_value;
  360. #else
  361.     return !(*this == op);
  362. #endif
  363. }
  364.  
  365. bool IEEEfloat::operator< (IEEEfloat op)
  366. {
  367.     // FIXME: This could be sped up by inlining
  368. #ifdef HAVE_IEEE754
  369.     return (value.float_value < op.value.float_value);
  370. #else
  371.     if (IsNaN() || op.IsNaN())
  372.         return false; // NaNs are unordered
  373.     if (IsZero())
  374.         return op.IsZero() ? false : op.IsPositive();
  375.     if (op.IsZero())
  376.         return IsNegative();
  377.     // Exploit fact that remaining IEEE floating point sort as signed ints
  378.     return value.iword < op.value.iword;
  379. #endif
  380. }
  381.  
  382. bool IEEEfloat::operator<= (IEEEfloat op)
  383. {
  384.     // FIXME: This could be sped up by inlining
  385. #ifdef HAVE_IEEE754
  386.     return (value.float_value <= op.value.float_value);
  387. #else
  388.     return *this < op || *this == op;
  389. #endif
  390. }
  391.  
  392. bool IEEEfloat::operator> (IEEEfloat op)
  393. {
  394.     // FIXME: This could be sped up by inlining
  395. #ifdef HAVE_IEEE754
  396.     return (value.float_value > op.value.float_value);
  397. #else
  398.     if (IsNaN() || op.IsNaN())
  399.         return false; // NaNs are unordered.
  400.     if (IsZero())
  401.         return op.IsZero() ? false : op.IsNegative();
  402.     if (op.IsZero())
  403.         return IsPositive();
  404.     // Exploit fact that remaining IEEE floating point sort as signed ints
  405.     return value.iword > op.value.iword;
  406. #endif
  407. }
  408.  
  409. bool IEEEfloat::operator>= (IEEEfloat op)
  410. {
  411.     // FIXME: This could be sped up by inlining
  412. #ifdef HAVE_IEEE754
  413.     return (value.float_value >= op.value.float_value);
  414. #else
  415.     return *this > op || *this == op;
  416. #endif
  417. }
  418.  
  419. IEEEfloat IEEEfloat::operator+ (IEEEfloat op)
  420. {
  421. #ifdef HAVE_IEEE754
  422.     // FIXME: This could be sped up by inlining
  423.     return IEEEfloat(value.float_value + op.value.float_value);
  424. #else
  425.     if (IsNaN() || op.IsNaN())
  426.         return NaN(); // arithmetic on NaNs not allowed
  427.  
  428.     //
  429.     // Adding unlike infinities not allowed
  430.     // Adding infinities of same sign is infinity of that sign
  431.     // Adding finite and infinity produces infinity
  432.     //
  433.     if (IsInfinite())
  434.         return op.IsInfinite() && (Sign() != op.Sign()) ? NaN() : *this;
  435.     if (op.IsInfinite())
  436.         return op;
  437.  
  438.     //
  439.     // Adding zero is easy
  440.     //
  441.     if (IsZero())
  442.         return (op.IsZero()) ? (Sign() != op.Sign()) ? POSITIVE_ZERO()
  443.                                                      : *this
  444.                              : op;
  445.     if (op.IsZero())
  446.         return *this;
  447.  
  448.     //
  449.     // Now for the real work - do manipulations on copies
  450.     //
  451.     i4 x, y, round = 0;
  452.     int expx, expy, signx, signy;
  453.     
  454.     expx = SplitInto((u4 &) x);
  455.     expy = op.SplitInto((u4 &) y);
  456.     
  457.     signx = Sign();
  458.     signy = op.Sign();
  459.  
  460.     // If the exponents are far enough apart, the answer is easy
  461.     if (expx > expy + 25)
  462.         return *this;
  463.     if (expy > expx + 25)
  464.         return op;
  465.  
  466.     //
  467.     // Denormalize the fractions, so that the exponents are
  468.     // the same and then set the exponent for the result.
  469.     // Leave enough space for overflow and INT_MIN avoidance!
  470.     // 
  471.     if (signx)
  472.         x = -x;
  473.     if (signy)
  474.         y = -y;
  475.  
  476.     x <<= 6;
  477.     y <<= 6;
  478.  
  479.     if (expx > expy) {
  480.         round = y << (32 + expy - expx);
  481.         y >>= expx - expy;
  482.     }
  483.     else if (expy > expx)
  484.     {
  485.         round = x << (32 + expx - expy);
  486.         x >>= expy - expx;
  487.         expx = expy;
  488.     }
  489.  
  490.     //
  491.     // Do the arithmetic. The excess magnitude of 32-bit arithmetic means
  492.     // overflow is impossible (we only need 1 spare bit!). We ensure that
  493.     // pre-alignment avoids any question of INT_MIN negation problems.
  494.     //
  495.     x += y;
  496.     if (round)
  497.         x |= 1;
  498.  
  499.     //
  500.     // If the result is negative, then make the fraction positive again
  501.     // and remember the sign.
  502.     //
  503.     if (x < 0)
  504.     {
  505.         x = -x;
  506.         signx = 1;
  507.     }
  508.     else
  509.         signx = 0;
  510.  
  511.     if (x == 0)
  512.         return signx ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  513.  
  514.     //
  515.     // Time to normalize again! If we need to shift left (the addition was
  516.     // effectively a subtraction), then there cannot be any reason to round.
  517.     // If the number fits exactly we don't have anything to do either.
  518.     //
  519.     return Normalize(signx, expx - 6, (u4) x);
  520. #endif
  521. }
  522.  
  523. IEEEfloat IEEEfloat::operator- ()
  524. {
  525.     // FIXME: This could be sped up by inlining
  526. #ifdef HAVE_IEEE754
  527.     return IEEEfloat(-value.float_value);
  528. #else
  529.     if (IsNaN())
  530.         return *this;
  531.     return IEEEfloat(value.word ^ 0x80000000);
  532. #endif
  533. }
  534.  
  535. IEEEfloat IEEEfloat::operator* (IEEEfloat op)
  536. {
  537. #ifdef HAVE_IEEE754
  538.     return IEEEfloat(value.float_value * op.value.float_value);
  539. #else
  540.     if (IsNaN() || op.IsNaN())
  541.         return NaN(); // arithmetic on NaNs not allowed
  542.  
  543.     int sign = Sign() ^ op.Sign();
  544.  
  545.     //
  546.     // If either operand is zero or infinite, then the answer is easy.
  547.     //
  548.     if (IsZero())
  549.         return op.IsInfinite() ? NaN()
  550.                                : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  551.  
  552.     if (op.IsZero())
  553.         return IsInfinite() ? NaN()
  554.                             : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  555.  
  556.     if (IsInfinite() || op.IsInfinite())
  557.         return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
  558.  
  559.     //
  560.     // Now for the real work - do manipulations on copies
  561.     //
  562.     u4 x, y;
  563.     int exponent;
  564.  
  565.     exponent = SplitInto(x) + op.SplitInto(y);
  566.  
  567.     //
  568.     // The numbers to be multiplied are 24 bits in length (stored in 32 bit
  569.     // integers). Using ULongInt to perform the multiplication, the result
  570.     // will be 46-48 bits (unsigned); shift it back to 28 bits for Normalize,
  571.     // while folding the low 20 bits into the lsb for rounding purposes.
  572.     //
  573.  
  574.     ULongInt a = x,
  575.              b = y;
  576.     a *= b;
  577.     b = a & 0x000FFFFF;
  578.     a >>= 20;
  579.     x = a.LowWord() | ((b > 0) ? 1 : 0);
  580.     
  581.     return Normalize(sign, exponent - 3, x);
  582. #endif // HAVE_IEEE754
  583. }
  584.  
  585. IEEEfloat IEEEfloat::operator/ (IEEEfloat op)
  586. {
  587. #ifdef HAVE_IEEE754
  588.     return op.IsZero() ? ((IsNaN() || IsZero()) ? NaN()
  589.                                : (IsPositive() ^ op.IsPositive())
  590.                                        ? NEGATIVE_INFINITY()
  591.                                        : POSITIVE_INFINITY())
  592.                        : IEEEfloat(value.float_value / op.value.float_value);
  593. #else // HAVE_IEEE754
  594.     if (IsNaN() || op.IsNaN())
  595.         return NaN(); // arithmetic on NaNs not allowed
  596.  
  597.     int sign = Sign() ^ op.Sign();
  598.  
  599.     //
  600.     // Infinities and zeroes are special.
  601.     //
  602.  
  603.     if (IsInfinite())
  604.         return op.IsInfinite() ? NaN()
  605.                                : sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
  606.  
  607.     if (op.IsInfinite())
  608.         return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  609.  
  610.     if (IsZero())
  611.         return op.IsZero() ? NaN()
  612.                            : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  613.  
  614.     if (op.IsZero())
  615.         return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
  616.  
  617.     //
  618.     // Now for the real work - do manipulations on copies
  619.     //
  620.     u4 x, y;
  621.     int exponent;
  622.  
  623.     exponent = SplitInto(x) - op.SplitInto(y);
  624.  
  625.     u4  mask   = 0x80000000,
  626.         result = 0;
  627.  
  628.     // Because both values are normalised, a single shift guarantees results.
  629.  
  630.     if (x < y)
  631.     {
  632.     x <<= 1;
  633.     exponent--;
  634.     }
  635.  
  636.     //
  637.     // If the numerator is larger, then it is divisible.
  638.     // Reflect this in the result, and do the subtraction.
  639.     // Magnify the numerator again and reduce the mask.
  640.     //
  641.     while (mask)
  642.     {
  643.     if (x >= y)
  644.     {
  645.         result += mask;
  646.         x -= y;
  647.         if (x == 0)
  648.         break;
  649.     }
  650.         x <<= 1;
  651.         mask >>= 1;
  652.     }
  653.  
  654.     return Normalize(sign, exponent - 8, result);
  655. #endif // HAVE_IEEE754
  656. }
  657.  
  658. IEEEfloat IEEEfloat::operator% (IEEEfloat op)
  659. {
  660. #ifdef HAVE_IEEE754
  661.     return IEEEfloat((op.IsZero() ? NaN().value.float_value
  662.                                   : (float) fmod((double) value.float_value, (double) op.value.float_value)));
  663. #else // HAVE_IEEE754
  664.     if (IsNaN() || op.IsNaN())
  665.         return NaN(); // arithmetic on NaNs not allowed
  666.  
  667.     //
  668.     // Infinities and zeroes are special.
  669.     //
  670.  
  671.     if (IsInfinite() || op.IsZero())
  672.         return NaN();
  673.  
  674.     if (IsZero() || op.IsInfinite())
  675.         return *this;
  676.  
  677.     //
  678.     // Now for the real work - do manipulations on copies
  679.     // This algorithm is from fdlibm.c - see above notice
  680.     //
  681.  
  682.     int sign = Sign();
  683.     u4 x, y;
  684.     int expy, exponent;
  685.     i4 z;
  686.  
  687.     expy = op.SplitInto(y);
  688.     exponent = SplitInto(x) - expy;
  689.  
  690.     if (exponent < 0 || (exponent == 0 && x < y))
  691.         return *this;
  692.  
  693.     while (exponent--)
  694.     {
  695.         z = x - y;
  696.         if (z < 0)
  697.             x <<= 1;
  698.         else if (z == 0)
  699.             return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  700.         else
  701.             x = z + z;
  702.     }
  703.     z = x - y;
  704.     if (z >= 0)
  705.         x = z;
  706.  
  707.     if (x == 0)
  708.         return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  709.     return Normalize(sign, expy, x);
  710. #endif // HAVE_IEEE754
  711. }
  712.  
  713.  
  714.  
  715. IEEEdouble::IEEEdouble(IEEEfloat f)
  716. {
  717. #ifdef HAVE_IEEE754
  718.     value.double_value = (double) f.FloatView();
  719. #else
  720.     int sign = f.Sign();
  721.     int exponent = f.Exponent();
  722.  
  723.     if (exponent == -127)
  724.     {
  725.     if (f.IsZero())
  726.             *this = sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  727.     else
  728.         {
  729.             //
  730.             // This is a denormalized number, with exponent -126
  731.             // (1 - IEEEfloat::BIAS); shift it to fit double format
  732.             //
  733.         *this = Normalize(sign, -126, ULongInt(f.Fraction()) << 29);
  734.     }
  735.     }
  736.     else
  737.     {
  738.         if (f.IsInfinite())
  739.             *this = sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
  740.         else if (f.IsNaN())
  741.             *this = NaN();
  742.         else
  743.     {
  744.             // Regular, noramlized number.  Shift it to fit double format
  745.             *this = Normalize(sign, exponent, ULongInt(f.Fraction()) << 29);
  746.     }
  747.     }
  748. #endif // HAVE_IEEE754
  749. }
  750.  
  751. IEEEdouble::IEEEdouble(i4 a)
  752. {
  753. #ifdef HAVE_IEEE754
  754.     value.double_value = (double) a;
  755. #else
  756.     int    sign = 0;
  757.  
  758.     if (a < 0)
  759.     {
  760.         a  = -a; // even works for MIN_INT!
  761.         sign = 1;
  762.     }
  763.     if (a == 0)
  764.         *this = POSITIVE_ZERO();
  765.     else
  766.         *this = Normalize(sign, FRACT_BITS, ULongInt((u4) a));
  767. #endif // HAVE_IEEE754
  768. }
  769.  
  770. IEEEdouble::IEEEdouble(LongInt a)
  771. {
  772. #ifdef HAVE_IEEE754
  773. # ifdef HAVE_UNSIGNED_LONG_LONG
  774.     value.double_value = (double) a.Words();
  775. # else
  776.     value.double_value = ((double) a.HighWord() * (double) 0x40000000 * 4.0) +
  777.         (double) a.LowWord();
  778. # endif // HAVE_UNSIGNED_LONG_LONG
  779. #else
  780.     int    sign = 0;
  781.  
  782.     if (a < 0)
  783.     {
  784.         a  = -a; // even works for MIN_LONG!
  785.         sign = 1;
  786.     }
  787.     if (a == 0)
  788.         *this = POSITIVE_ZERO();
  789.     else
  790.         *this = Normalize(sign, FRACT_BITS, ULongInt(a));
  791. #endif // HAVE_IEEE754
  792. }
  793.  
  794. IEEEdouble::IEEEdouble(char *str, bool check_invalid)
  795. {
  796.     //
  797.     // TODO: This conversion is a temporary patch. Need volunteer to implement
  798.     //       Clinger algorithm from PLDI 1990.
  799.     //
  800.     value.double_value = atof(str);
  801.  
  802.     //
  803.     // When parsing a literal in Java, a number that rounds to infinity or
  804.     // zero is invalid.  A true check_invalid sets any invalid number to NaN,
  805.     // to make the upstream processing easier.  Leave check_invalid false to
  806.     // allow parsing infinity or rounded zero from a string.
  807.     //
  808.     if (check_invalid)
  809.     {
  810.         if (IsInfinite())
  811.             *this = NaN();
  812.         if (IsZero())
  813.         {
  814.             for (char *p = str; *p; p++) {
  815.                 switch (*p) {
  816.                 case '-': case '+': case '0': case '.':
  817.                     break; // keep checking
  818.                 case 'e': case 'E': case 'd': case 'D': case 'f': case 'F':
  819.                     return; // got this far, str is 0
  820.                 default:
  821.                     *this = NaN(); // encountered non-zero digit. oops.
  822.                     return;
  823.                 }
  824.             }
  825.         }
  826.     }
  827. }
  828.  
  829. i4 IEEEdouble::IntValue()
  830. {
  831.     if (IsNaN())
  832.         return 0;                                                             
  833.             
  834. #ifdef HAVE_IEEE754
  835.     if (value.double_value < MIN_INT)
  836.         return MIN_INT;
  837.     else if (value.double_value > MAX_INT)
  838.         return MAX_INT;
  839.     return (i4) value.double_value;
  840. #else
  841.     int    sign = Sign(),
  842.         exponent = Exponent();
  843.  
  844.     if (IsInfinite() || exponent > 30)
  845.         return sign ? MIN_INT : MAX_INT;
  846.  
  847.     // This includes true zero and denorms.
  848.     if (exponent < 0)
  849.         return 0;
  850.  
  851.     i4 result = (i4) (Fraction() >> (FRACT_BITS - exponent)).LowWord();
  852.  
  853.     return sign ? -result : result;
  854. #endif // HAVE_IEEE754
  855. }
  856.  
  857. LongInt IEEEdouble::LongValue()
  858. {
  859.     if (IsNaN())
  860.         return LongInt(0);
  861.             
  862.     int    sign = Sign(),
  863.         exponent = Exponent();
  864.  
  865.     if (IsInfinite() || exponent > 62)
  866.         return sign ? LongInt::MIN_LONG() : LongInt::MAX_LONG();
  867.  
  868.     // This covers true zero and denorms.
  869.     if (exponent < 0)
  870.         return LongInt(0);
  871.  
  872.     LongInt result = Fraction();
  873.  
  874.     if (exponent > FRACT_BITS)
  875.     result <<= (exponent - FRACT_BITS);
  876.     else if (exponent < FRACT_BITS)
  877.     result >>= (FRACT_BITS - exponent);
  878.  
  879.     return sign ? (LongInt) -result : result;
  880. }
  881.  
  882. IEEEdouble IEEEdouble::Normalize(int sign, int exponent, ULongInt fraction)
  883. {
  884.     bool round = false, sticky = false;
  885.  
  886.     assert(fraction != 0);
  887.  
  888.     //
  889.     // Normalize right. Shift until value < MAX_FRACT.
  890.     //
  891.     if (fraction.HighWord() >= MAX_FRACT)
  892.     {
  893.     while (fraction.HighWord() >= MAX_FRACT)
  894.     {
  895.         sticky |= round;
  896.         round = (fraction.LowWord() & 0x00000001) != 0;
  897.         fraction >>= 1;
  898.         exponent++;
  899.     }
  900.     if (round && (sticky || (fraction.LowWord() & 0x00000001)) && exponent > -BIAS)
  901.             //
  902.             // Capture any overflow caused by rounding. No other checks are
  903.             // required because if overflow occurred, the the low order bit
  904.             // was guaranteed to be zero. Do not round denorms yet.
  905.             //
  906.         if ((++fraction).HighWord() >= MAX_FRACT)
  907.         {
  908.         fraction >>= 1;
  909.         exponent++;
  910.         }
  911.     }
  912.  
  913.     // Normalize left. MAX_FRACT2 is MAX_FRACT >> 1.
  914.     else
  915.     while (fraction.HighWord() < MAX_FRACT2)
  916.     {
  917.         fraction <<= 1;
  918.         exponent--;
  919.     }
  920.  
  921.     //
  922.     // Check and respond to overflow
  923.     //
  924.     if (exponent > BIAS)
  925.         return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
  926.  
  927.     //
  928.     // Check and respond to underflow 
  929.     //
  930.     if (exponent <= -BIAS)
  931.     {
  932.     while (exponent <= -BIAS)
  933.     {
  934.         sticky |= round;
  935.         round = (fraction.LowWord() & 0x00000001) != 0;
  936.         fraction >>= 1;
  937.         exponent++;
  938.     }
  939.  
  940.     if (round && (sticky || (fraction.LowWord() & 0x00000001)))
  941.         fraction++;
  942.  
  943.     exponent = -BIAS;
  944.  
  945.     if (fraction == 0)
  946.             return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  947.     }
  948.  
  949.     return IEEEdouble((sign << 31) | ((exponent + BIAS) << 20)
  950.                                    | (fraction.HighWord() & 0x000FFFFF),
  951.                       fraction.LowWord());
  952. }
  953.  
  954. int IEEEdouble::SplitInto(BaseLong &fraction)
  955. {
  956.     int exponent = Exponent();
  957.     fraction = Fraction();
  958.  
  959.     if (exponent == -BIAS)
  960.     {
  961.     exponent++;
  962.     while (fraction.HighWord() < MAX_FRACT2)
  963.     {
  964.         fraction <<= 1;
  965.         exponent--;
  966.     }
  967.     }
  968.  
  969.     return exponent;
  970. }
  971.  
  972. bool IEEEdouble::operator== (IEEEdouble op)
  973. {
  974.     // FIXME: This could be sped up by inlining
  975. #ifdef HAVE_IEEE754
  976.     // TODO: Microsoft VC++ botches this, mixing 12.0 and NaN
  977.     return value.double_value == op.value.double_value;
  978. #else
  979.     return IsNaN() || op.IsNaN() ? false
  980.                                  : IsZero() && op.IsZero() ? true
  981.                                                            : (BaseLong) *this == (BaseLong) op;
  982. #endif
  983. }
  984.  
  985. bool IEEEdouble::operator!= (IEEEdouble op)
  986. {
  987.     // FIXME: This could be sped up by inlining
  988. #ifdef HAVE_IEEE754
  989.     return value.double_value != op.value.double_value;
  990. #else
  991.     return !(*this == op);
  992. #endif
  993. }
  994.  
  995. IEEEdouble IEEEdouble::operator+ (IEEEdouble op)
  996. {
  997. #ifdef HAVE_IEEE754
  998.     // FIXME: This could be sped up by inlining
  999.     return IEEEdouble(value.double_value + op.value.double_value);
  1000. #else
  1001.     if (IsNaN() || op.IsNaN())
  1002.         return NaN(); // arithmetic on NaNs not allowed
  1003.  
  1004.     //
  1005.     // Adding unlike infinities not allowed
  1006.     // Adding infinities of same sign is infinity of that sign
  1007.     // Adding finite and infinity produces infinity
  1008.     //
  1009.     if (IsInfinite())
  1010.         return op.IsInfinite() && (Sign() != op.Sign()) ? NaN() : *this;
  1011.     if (op.IsInfinite())
  1012.         return op;
  1013.  
  1014.     //
  1015.     // Adding zero is easy
  1016.     //
  1017.     if (IsZero())
  1018.         return (op.IsZero()) ? (Sign() != op.Sign()) ? POSITIVE_ZERO()
  1019.                                                      : *this
  1020.                              : op;
  1021.     if (op.IsZero())
  1022.         return *this;
  1023.  
  1024.     //
  1025.     // Now for the real work - do manipulations on copies
  1026.     //
  1027.     LongInt x, y, round = 0;
  1028.     int expx, expy, signx, signy;
  1029.     
  1030.     expx = SplitInto(x);
  1031.     expy = op.SplitInto(y);
  1032.     
  1033.     signx = Sign();
  1034.     signy = op.Sign();
  1035.  
  1036.     // If the exponents are far enough apart, the answer is easy
  1037.     if (expx > expy + 54)
  1038.         return *this;
  1039.     if (expy > expx + 54)
  1040.         return op;
  1041.  
  1042.     //
  1043.     // Denormalize the fractions, so that the exponents are
  1044.     // the same and then set the exponent for the result.
  1045.     // Leave enough space for overflow and LONG_MIN avoidance!
  1046.     // 
  1047.     if (signx)
  1048.         x = -x;
  1049.     if (signy)
  1050.         y = -y;
  1051.  
  1052.     x <<= 8;
  1053.     y <<= 8;
  1054.  
  1055.     if (expx > expy)
  1056.     {
  1057.         round = y << (64 + expy - expx);
  1058.         y >>= expx - expy;
  1059.     }
  1060.     else if (expy > expx)
  1061.     {
  1062.         round = x << (64 + expx - expy);
  1063.         x >>= expy - expx;
  1064.         expx = expy;
  1065.     }
  1066.  
  1067.     //
  1068.     // Do the arithmetic. The excess magnitude of 64-bit arithmetic means
  1069.     // overflow is impossible (we only need 1 spare bit!). We ensure that
  1070.     // pre-alignment avoids any question of LONG_MIN negation problems.
  1071.     //
  1072.     x += y;
  1073.     if (round != 0)
  1074.         x |= 1;
  1075.  
  1076.     //
  1077.     // If the result is negative, then make the fraction positive again
  1078.     // and remember the sign.
  1079.     //
  1080.     if (x < 0)
  1081.     {
  1082.         x = -x;
  1083.         signx = 1;
  1084.     }
  1085.     else
  1086.         signx = 0;
  1087.  
  1088.     if (x == 0)
  1089.         return signx ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  1090.  
  1091.     //
  1092.     // Time to normalize again! If we need to shift left (the addition was
  1093.     // effectively a subtraction), then there cannot be any reason to round.
  1094.     // If the number fits exactly we don't have anything to do either.
  1095.     //
  1096.     return Normalize(signx, expx - 8, (ULongInt) x);
  1097. #endif
  1098. }
  1099.  
  1100. IEEEdouble IEEEdouble::operator- ()
  1101. {
  1102.     // FIXME: This could be sped up by inlining
  1103. #ifdef HAVE_IEEE754
  1104.     return IEEEdouble(-value.double_value);
  1105. #else
  1106.     if (IsNaN())
  1107.         return *this;
  1108.     return IEEEdouble(HighWord() ^ 0x80000000, LowWord());
  1109. #endif // HAVE_IEEE754
  1110. }
  1111.  
  1112. IEEEdouble IEEEdouble::operator* (IEEEdouble op)
  1113. {
  1114. #ifdef HAVE_IEEE754
  1115.     return IEEEdouble(value.double_value * op.value.double_value);
  1116. #else
  1117.     if (IsNaN() || op.IsNaN())
  1118.         return NaN(); // arithmetic on NaNs not allowed
  1119.  
  1120.     int sign = Sign() ^ op.Sign();
  1121.  
  1122.     //
  1123.     // If either operand is zero or infinite, then the answer is easy.
  1124.     //
  1125.     if (IsZero())
  1126.         return op.IsInfinite() ? NaN()
  1127.                                : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  1128.  
  1129.     if (op.IsZero())
  1130.         return IsInfinite() ? NaN()
  1131.                             : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  1132.  
  1133.     if (IsInfinite() || op.IsInfinite())
  1134.         return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
  1135.  
  1136.     //
  1137.     // Now for the real work - do manipulations on copies
  1138.     //
  1139.     ULongInt x, y, z, w, pr1, pr2;
  1140.     u4 round;
  1141.     int exponent;
  1142.  
  1143.     exponent = SplitInto(x) + op.SplitInto(y);
  1144.  
  1145.     //
  1146.     // The numbers to be multiplied are 53 bits in length (stored in 64 bit
  1147.     // integers). Split them in 32 bit parts, then shift result into place.
  1148.     // Fold the low order bits into the lsb for rounding purposes.
  1149.     //
  1150.     z = ULongInt(x.LowWord());
  1151.     w = ULongInt(y.LowWord());
  1152.     x >>= 32;
  1153.     y >>= 32;
  1154.  
  1155.     pr1 = z * w; // low*low
  1156.     pr2 = z * y + x * w + pr1.HighWord(); // low*high + high*low + adjusted pr1
  1157.     round = pr1.LowWord();
  1158.     x = (x * y) << 20; // high*high, adjusted
  1159.  
  1160.     round |= pr2.LowWord() & 0xFFF;
  1161.     pr2 >>= 12;
  1162.  
  1163.     x += pr2 | (round ? 1 : 0);
  1164.  
  1165.     return Normalize(sign, exponent - 8, x);
  1166. #endif // HAVE_IEEE754
  1167. }
  1168.  
  1169. IEEEdouble IEEEdouble::operator/ (IEEEdouble op)
  1170. {
  1171. #ifdef HAVE_IEEE754
  1172.     return op.IsZero() ? ((IsNaN() || IsZero()) ? NaN()
  1173.                                : (IsPositive() ^ op.IsPositive())
  1174.                                        ? NEGATIVE_INFINITY()
  1175.                                        : POSITIVE_INFINITY())
  1176.                        : IEEEdouble(value.double_value / op.value.double_value);
  1177. #else // HAVE_IEEE754
  1178.     if (IsNaN() || op.IsNaN())
  1179.         return NaN(); // arithmetic on NaNs not allowed
  1180.  
  1181.     int sign = Sign() ^ op.Sign();
  1182.  
  1183.     //
  1184.     // Infinities and zeroes are special.
  1185.     //
  1186.  
  1187.     if (IsInfinite())
  1188.         return op.IsInfinite() ? NaN()
  1189.                                : sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
  1190.  
  1191.     if (op.IsInfinite())
  1192.         return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  1193.  
  1194.     if (IsZero())
  1195.         return op.IsZero() ? NaN()
  1196.                            : sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  1197.  
  1198.     if (op.IsZero())
  1199.         return sign ? NEGATIVE_INFINITY() : POSITIVE_INFINITY();
  1200.  
  1201.     //
  1202.     // Now for the real work - do manipulations on copies
  1203.     //
  1204.     ULongInt x, y;
  1205.     int exponent;
  1206.  
  1207.     exponent = SplitInto(x) - op.SplitInto(y);
  1208.  
  1209.     ULongInt mask  (0x80000000, 0x00000000),
  1210.              result(0x00000000, 0x00000000);
  1211.  
  1212.     // Because both values are normalised, a single shift guarantees results.
  1213.  
  1214.     if (x < y)
  1215.     {
  1216.     x <<= 1;
  1217.     exponent--;
  1218.     }
  1219.  
  1220.     //
  1221.     // If the numerator is larger, then it is divisible.
  1222.     // Reflect this in the result, and do the subtraction.
  1223.     // Magnify the numerator again and reduce the mask.
  1224.     //
  1225.     while (mask != 0)
  1226.     {
  1227.     if (x >= y)
  1228.     {
  1229.         result += mask;
  1230.         x -= y;
  1231.         if (x == 0)
  1232.         break;
  1233.     }
  1234.         x <<= 1;
  1235.         mask >>= 1;
  1236.     }
  1237.  
  1238.     return Normalize(sign, exponent - 11, result);
  1239. #endif // HAVE_IEEE754
  1240. }
  1241.  
  1242. IEEEdouble IEEEdouble::operator% (IEEEdouble op)
  1243. {
  1244. #ifdef HAVE_IEEE754
  1245.     return IEEEdouble((op.IsZero() ? NaN().value.double_value
  1246.                                    : fmod(value.double_value, op.value.double_value)));
  1247. #else // HAVE_IEEE754
  1248.     if (IsNaN() || op.IsNaN())
  1249.         return NaN(); // arithmetic on NaNs not allowed
  1250.  
  1251.     //
  1252.     // Infinities and zeroes are special.
  1253.     //
  1254.  
  1255.     if (IsInfinite() || op.IsZero())
  1256.         return NaN();
  1257.  
  1258.     if (IsZero() || op.IsInfinite())
  1259.         return *this;
  1260.  
  1261.     //
  1262.     // Now for the real work - do manipulations on copies
  1263.     // This algorithm is from fdlibm.c - see above notice
  1264.     //
  1265.  
  1266.     int sign = Sign();
  1267.     ULongInt x, y;
  1268.     int expy, exponent;
  1269.     LongInt z;
  1270.  
  1271.     expy = op.SplitInto(y);
  1272.     exponent = SplitInto(x) - expy;
  1273.  
  1274.     if (exponent < 0 || (exponent == 0 && x < y))
  1275.         return *this;
  1276.  
  1277.     while (exponent--)
  1278.     {
  1279.         z = x - y;
  1280.         if (z < 0)
  1281.             x <<= 1;
  1282.         else if (z == 0)
  1283.             return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  1284.         else
  1285.             x = z + z;
  1286.     }
  1287.     z = x - y;
  1288.     if (z >= 0)
  1289.         x = z;
  1290.  
  1291.     if (x == 0)
  1292.         return sign ? NEGATIVE_ZERO() : POSITIVE_ZERO();
  1293.     return Normalize(sign, expy, x);
  1294. #endif // HAVE_IEEE754
  1295. }
  1296.  
  1297. bool IEEEdouble::operator< (IEEEdouble op)
  1298. {
  1299.     // FIXME: This could be sped up by inlining
  1300. #ifdef HAVE_IEEE754
  1301.     return (value.double_value < op.value.double_value);
  1302. #else
  1303.     if (IsNaN() || op.IsNaN())
  1304.         return false; // NaNs are unordered
  1305.     if (IsZero())
  1306.         return op.IsZero() ? false : op.IsPositive();
  1307.     if (op.IsZero())
  1308.         return IsNegative();
  1309.     // Exploit fact that remaining IEEE floating point sort as signed ints
  1310.     u4 a = HighWord(), b = op.HighWord();
  1311.     return (a < b) || ((a == b) && LowWord() < op.LowWord());
  1312. #endif
  1313. }
  1314.  
  1315. bool IEEEdouble::operator<= (IEEEdouble op)
  1316. {
  1317.     // FIXME: This could be sped up by inlining
  1318. #ifdef HAVE_IEEE754
  1319.     return (value.double_value <= op.value.double_value);
  1320. #else
  1321.     return *this < op || *this == op;
  1322. #endif
  1323. }
  1324.  
  1325. bool IEEEdouble::operator> (IEEEdouble op)
  1326. {
  1327.     // FIXME: This could be sped up by inlining
  1328. #ifdef HAVE_IEEE754
  1329.     return (value.double_value > op.value.double_value);
  1330. #else
  1331.     if (IsNaN() || op.IsNaN())
  1332.         return false; // NaNs are unordered.
  1333.     if (IsZero())
  1334.         return op.IsZero() ? false : op.IsNegative();
  1335.     if (op.IsZero())
  1336.         return IsPositive();
  1337.     // Exploit fact that remaining IEEE floating point sort as signed ints
  1338.     u4 a = HighWord(), b = op.HighWord();
  1339.     return (a > b) || ((a == b) && LowWord() > op.LowWord());
  1340. #endif
  1341. }
  1342.  
  1343. bool IEEEdouble::operator>= (IEEEdouble op)
  1344. {
  1345.     // FIXME: This could be sped up by inlining
  1346. #ifdef HAVE_IEEE754
  1347.     return (value.double_value >= op.value.double_value);
  1348. #else
  1349.     return *this > op || *this == op;
  1350. #endif
  1351. }
  1352.  
  1353. #ifdef    HAVE_JIKES_NAMESPACE
  1354. }            // Close namespace Jikes block
  1355. #endif
  1356.  
  1357.