New Long Emulation -> Interest in Math/Double missing methods in Double/Math?

[Ray Cromwell]

Now that the GWT compiler supports long types in 1.5 with 64 bit accuracy, is there any interest to have implementations of the following methods?

Double.doubleToLongBits()
Double.longBitsToDouble()
Math.getExponent()
Math.IEEEremainder()
Math.cbrt()
Math.ulp()
Math.nextAfter()
Math.nextUp()

I have reasonably efficient implementations of the first 4 functions (meaning not brute force loops) which are probably not bit-for-bit ulp-for-ulp correct vs x86 or fdlibm, but I haven't found a failure case yet. The first 2 enable the last 3 to be implemented easily. Of course, doubleToLongBits won't be as efficient as C (not are the emulated longs), but if you need 2-way reversible longbits->double functions, they might suffice.

Here is a sample implementation of doubleToLongBits(), it contains loops, but the initial guess gets the approximation within a factor of 10, so the loop will usually only a few iterations, usually 2. The implementations of
getExponent() and IEEEremainder() are efficient. Whether they are 100% correct remains to be seen, but they seem to be very good approximations.

-Ray


public static long doubleToLongBits(final double v) {
    if (Double.isNaN(v)) {
      // IEEE754, NaN exponent bits all 1s, and mantissa is non-zero
      return 0x0FFFl << 51;
    }

    long sign = (v < 0 ? 0x1l << 63 : 0);
    long exponent = 0;

    double absV = Math.abs(v);
    // IEEE754 infinite numbers, exponent all 1s, mantissa is 0
    if (Double.isInfinite(v)) {
      exponent = 0x07FFl << 52;
    } else {
      if (absV == 0.0) {
        // IEEE754, exponent is 0, mantissa is zero
        // we don't handle negative zero at the moment, it is treated as
        // positive zero
        exponent = 0l;
      } else {
        // get an approximation to the exponent
        int guess = (int) Math.floor(Math.log(absV) / Math.log(2));
        // force it to -1023, 1023 interval (<= -1023 = denorm/zero)
        guess = Math.max(-1023, Math.min(guess, 1023));
     
        // divide away exponent guess
        double exp = Math.pow(2, guess);
        absV = absV / exp;
       
        // while the number is still bigger than a normalized number
        // increment exponent guess
        while (absV > 2.0) {
          guess++;
          absV /= 2.0;
        }
        // if the number is smaller than a normalized number
        // decrement exponent
        while (absV < 1 && guess > 1024) {
          guess--;
          absV *= 2;
        }
        exponent = (guess + 1023l) << 52;
      }
    }
    // if denorm
    if (exponent <= 0) {
      absV /= 2;
    }

    // the input value has now been stripped of its exponent
    // and is in the range [0,2), we strip off the leading decimal
    // and use the remainer as a percentage of the significand value (2^52)
    long mantissa = (long)((absV % 1) * Math.pow(2, 52));
    return sign | exponent | (mantissa & 0xfffffffffffffl);
  }

[Bruce Johnson]

That would rule. Time may be running short for GWT 1.5, but we'd be grateful for patches whether it makes it in or not.

[Ian Petersen]

Double.doubleToLongBits() would be nice to have because then Eclipse's
"Generate hashCode and equals" wizard could be used in client-side
code without having to modify the resulting code. Using it that way
would probably require caching results in Double, though (similar to
the HashCache class in GWT's String implementation).

Math.nextAfter and .nextUp confused me at first because I'm still on
Java 1.5 and they're part of 1.6, but they'd be useful in implementing
a java.text.ChoiceFormat work-alike.

I don't have a need for the others.

Ian

--
Tired of pop-ups, security holes, and spyware?
Try Firefox: http://www.getfirefox.com

[Ray Cromwell]

Ask and ye shall receive... :)

Attached is the implementation of Double.doubleToLongBits and Double.longBitsToDouble, plus patch to DoubleTest.java to include tests against actual values returned from JDK1.6 execution of those methods.

-Ray

[Ray Cromwell]

Opps, the test cases are misnamed and wrong (assertEquals on double comparisons), here's the fix.

-Ray

[Scott Blum]

Freeland or Lex, could one of you review?

[Ray Cromwell]

One comment is that those constant Math.pow()'s could be saved as a
static final constant (e.g. Math.pow(2, 52) ), likewise, the bitshift
masks can be constants, unless the GWT compiler evaluates them at
compile time.

-Ray

On Mar 19, 5:08 pm, "Scott Blum" <sco...@google.com> wrote:
> Freeland or Lex, could one of you review?
>

> On Wed, Mar 19, 2008 at 6:57 PM, Ray Cromwell <cromwell...@gmail.com> wrote:
>
> > Opps, the test cases are misnamed and wrong (assertEquals on double
> > comparisons), here's the fix.
>
> > -Ray
>

> > On Wed, Mar 19, 2008 at 3:41 PM, Ray Cromwell <cromwell...@gmail.com>

> > wrote:
>
> > > Ask and ye shall receive... :)
>
> > > Attached is the implementation of Double.doubleToLongBits and
> > > Double.longBitsToDouble, plus patch to DoubleTest.java to include tests
> > > against actual values returned from JDK1.6 execution of those methods.
>
> > > -Ray
>
> > > On Wed, Mar 19, 2008 at 2:15 PM, Bruce Johnson <br...@google.com> wrote:
>
> > > > That would rule. Time may be running short for GWT 1.5, but we'd be
> > > > grateful for patches whether it makes it in or not.
>

> > > > On Wed, Mar 19, 2008 at 4:55 PM, Ray Cromwell <cromwell...@gmail.com>

[Ray Cromwell]

During my implementation of methods in the Math class I discovered that the implementation of Math.hypot() has a bug. It won't correctly handle values of a and b which might overflow or underflow. For example, if either a or b is of order sqrt(MAX_VALUE) or sqrt(MIN_VALUE), then hypot() will incorrect return Infinity.

Here is a fixed implementation (which can probably be optimized a little bit, removing some of the special cases)

private static final double HALF_MAX_POWER = Math.pow(2, 500);
private static final double HALF_MIN_POWER = Math.pow(2, -500);

public static double hypot(double a, double b) {
        if(Double.isNaN(a-b)) {
            return Double.NaN;
        }
        else if(Math.abs(a-b) == a || Math.abs(a-b) == b) {
            return a+b;
        }
        // if overflow/underflow not possible
        else if(a > HALF_MIN_POWER && a < HALF_MAX_POWER
           && b > HALF_MIN_POWER && b < HALF_MAX_POWER) {
            return Math.sqrt(a*a + b*b);
        }
        // else, rewrite a*a + b*b as
        // hypot = a*a * (1 + b*b/a*a) where a > b
        // then sqrt(hypot) as a * sqrt(1 + b*b/a*a)
        double denom = a > b ? a : b;
        double ratio = a > b ? b / a : a / b;
        double scale = 1;
        // however if the ratio of b/a would squared
        // overflow/underflow, we need to scale
        // so we multiply by 2^-500 or 2^500 if the
        // exponent exceeds 2^500 or 2^-500 respectively
        if(ratio < HALF_MIN_POWER) {
            scale = HALF_MAX_POWER;
        }
        else if(ratio > HALF_MAX_POWER) {
            scale = HALF_MIN_POWER;
        }
        ratio *= scale;
        // hypot = sqrt(1/scale * a*a * scale * (1 + b*b/a*a))
        //       = sqrt(a*a/scale * (scale + b*b/a*a*scale)
        //       = a * sqrt(1/scale) * sqrt(scale + b*b/a*a*scale)
        return denom*Math.sqrt(1/scale)*Math.sqrt(scale + ratio*ratio);
    }

-Ray

[Ray Cromwell]

BTW, since the code doesn't document exactly why it is supposed to work, here's my 'theory' on how this algorithm works to help review. To tell the truth, I'm not sure how correct it is around the margins, since I coded it according to the theory, and then 'fixed' it when some tests didn't work as expected.

Let d = a double value, then it may be written as d = 1.m * 2 ^ e. That is, the leading mantissa bit is always a 1 (except in denorms), multiplied by a 2^e shift.

We first want to retrieve e, so we take log2 both sides (log2(x) = log(x)/log(2))

log2(d) = log2(1.m * 2^e) = log2(1.m) + log2(2^e) = log2(1.m) + e

Note: log2(1.m) is always less than 1

therefore, floor(log2(d)) = floor(log2(1.m) + e) = floor(e) = e

Once we have e, we now want to compute the mantissa, so going back to d = 1.m * 2^e

d/2^e = 1.m * 2^e / 2^e = 1.m = mantissa

Taking d/2^e % 1 strips off the leading 1 and gives us the fractional mantissa.

Now, if we forget about denorms, and perhaps some imprecision in the log2(x) operation, this seems clear that it will work to recover the mantissa and exponent.

Next, multiply the fractional mantissa * 2^52 to convert it to a binary fraction. Add 1023 to the exponent to bias it. Combine with sign, and shift into appropriate IEEE754 bit positions, and done.

Now, when the exponent is zero, you have a denorm number, so there is a slight subtle modification to normalize it.

If everything worked perfectly, those loops wouldn't be in there. Most of the time they never execute, but in testing, I ran into some weird 'off by 1' errors in exponent calculation that required adjustment.

-Ray

[Lex Spoon]

On Thu, Mar 20, 2008 at 4:44 AM, Ray Cromwell <cromw...@gmail.com> wrote:

BTW, since the code doesn't document exactly why it is supposed to work, here's my 'theory' on how this algorithm works to help review. To tell the truth, I'm not sure how correct it is around the margins, since I coded it according to the theory, and then 'fixed' it when some tests didn't work as expected.

In general, Ray, such notes are very helpful for getting your contribution reviewed.  So by all means share... or even put some key parts into the comments.  :)

-Lex

[Ray Cromwell]

 
No problem, if this patch is approved, I will submit a patch to Math.java that implements IEEEremainder, nextAfter/nextUp, ulp, getExponent, a fixed hypot(), along with mathematical explanations of how each is supposed to work. (Theory and reality are sometimes different :) )

-Ray

[Freeland Abbott]

(Patched in at r2201, if it matters for line numberings...)

Style question, I at least prefer capital "L" to indicate long constants, as lowercase 'l' looks like 1... does GWT code insist one way or the other?

doubleToLongBits,

• l.56, the "special cases" tests should be together... in particular, the test for Infinity should be near the test for NaN, and before the (irrelevant to that case) call to Math.abs().
• l.82, what's the "&& guess > 1024" clause for?  I get that we're checking whether our putative mantissa is too small (fractional) and we should ramp down the exponent to raise it to at least one, but where does the 1024 come from?  If the double is (fraction x 2^512), why wouldn't we want to adjust to (2*fraction x 2^256)?
• generally, we probably want symbolics for most of the constants.  I know why most are there, but balked at 1024 and at 10231... oh, duh, that last is 1023L, and that's also where 1024 comes from.  But I shouldn't have had to figure that out. ;-)
  

longBitsToDouble,

• name the input something other than 'l', for the same looks-like-a-1 reason I prefer constants ending in L,
• l.134 typo "out" -> "our" ?

Otherwise, I'm pretty convinced.

[Freeland Abbott]

Oh, and the test case should include more normal numbers, particularly something with negative exponents, and tests with min/max values for the separate fields but not min/max double.

In fact, we might fully permute {positive, negative} x { min/max/zero/arbitrary mantissa} x { min/max/arbitrary/zero exp}, plus NAN, +INF, and -INF.  (That's without thinking about which of the suggested permutations collapse into NaN or Infinity.)

[Ray Cromwell]

On Fri, Mar 21, 2008 at 1:58 PM, Freeland Abbott <gwt.team...@gmail.com> wrote:

(Patched in at r2201, if it matters for line numberings...)

Style question, I at least prefer capital "L" to indicate long constants, as lowercase 'l' looks like 1... does GWT code insist one way or the other?

100% agreed. I was too lazy to hit the shift key there. :)

doubleToLongBits,

• l.56, the "special cases" tests should be together... in particular, the test for Infinity should be near the test for NaN, and before the (irrelevant to that case) call to Math.abs().

Agreed, I'll fix.

• l.82, what's the "&& guess > 1024" clause for?  I get that we're checking whether our putative mantissa is too small (fractional) and we should ramp down the exponent to raise it to at least one, but where does the 1024 come from?  If the double is (fraction x 2^512), why wouldn't we want to adjust to (2*fraction x 2^256)?

 If you recall that the exponent is biased by +1023, this is effectively saying

while(absV < 1 && exponent > 0) { absV *= 2; exponent--; }

I agree it is unclear and in explaining this to you, I just noticed a bug. 1024 is not the exponent 0, 1023 is! I'll change this to ZERO_EXPONENT and set it to 1023.
The goal of that loop BTW is that if the leading decimal isn't 1, to keep shifting left by 2 (and decrementing the exponent) until it is so. However, if you shift left and exponent becomes zero,
then you have a denormalized number.

• generally, we probably want symbolics for most of the constants.  I know why most are there, but balked at 1024 and at 10231... oh, duh, that last is 1023L, and that's also where 1024 comes from.  But I shouldn't have had to figure that out. ;-)

Agreed.

• 
  

longBitsToDouble,

• name the input something other than 'l', for the same looks-like-a-1 reason I prefer constants ending in L,
• l.134 typo "out" -> "our" ?

Agreed (and typo). In general, I'm going to revise the comments and include my 'theory' post  on what each part does with respect to the theory.

Thanks for taking the time to review it.
-Ray

[Ray Cromwell]

Freeland,
  I agree. I also think we need to test all variations of denorms, plus negative-zero (bits for zero, but with negative sign bit), and mantissas with all bits set vs mantissas with 1 bit, but with non-max exponents. Most of my uncertainty over whether this is 100% correct revolves around whether log2(double) to get the exponent, and dividing off the exponent to get the mantissa doesn't truncate any bits.  For example, if you have a double value of 1.23456789 * 2^100 (let's say, exactly representable in the mantissa), and you do a division by 2^100, will the Javascript/Java IEEE division leave the mantissa intact and only affect the exponent? The answer seems to be yes in my tests, but I'm not 100% sure. Perhaps I could write a loop to run through a few hundred million samples comparing against the JDK version to see.

Even so, we don't neccessarily need our web-mode version to be bit-for-bit compatible with the JDK version, as long is the error is tolerable, say, converting to and from long values loses 1 bit of precision in the mantissa.

-Ray

[Ray Cromwell]

Freeland,
  Here is a revised version of the patch,  abstracted to work with arbitrary bit sites (you could cut/paste this code into Float to get a similar implementation). All former magic constants are now constant identifiers. I fixed some checkstyle issues, some bugs (mentioned in last post), and added extensive documentation so that what it's doing can be followed line by line for verification.

  One thing I was not able to do was to make the MAX_EXPONENT / MIN_EXPONENT constants dependent upon my code because gwt's checkstyle rules force private static fields to come after public, and that would create an illegal forward reference, otherwise, MAX_EXPONENT could be (1L << EXPONENT_BITSIZE) - 1 - EXPONENT_BIAS - 1 (e.g. 2048 - 1 - 1023 - 1 = 1023, a value of 1024 is reserved for NaN/Infinity). Note: If this patch is approved, then I will be submitting a patch to Math.java which includes a working getExponent() method which could also be used (e.g. getExponent(MAX_VALUE) )

  I did not update the tests yet, because I want to lock down the implementation first, then I'll go add extensive permutations.

-Ray

[Freeland Abbott]

Sorry to leave you hanging... this got pushed off from top-of-stack, and I apologize.

LGTM, although now I have to check where we are w.r.t. lock down for 1.5m2 before submit!

On Wed, Mar 26, 2008 at 3:12 PM, Ray Cromwell <cromw...@gmail.com> wrote:

damn it, I always forget to attach, here it is.

On Wed, Mar 26, 2008 at 12:12 PM, Ray Cromwell <cromw...@gmail.com> wrote:

Ok, I made it explicit so that the code is easier to understand. Here's an updated patch, major changes since last time you looked:

1) explicit assignment of absV = 0 in Infinite case
2) All magic numbers/constants changed to static final fields
3) static final fields are package protected so upcoming java.lang.Math patch can access them
4) added "Theory" comments everywhere
5) recognized special cases according to your comments
6) literal longs l -> L

This code works in the JDK, but I never tested it in GWT web-mode, because I haven't used a build with Long Emulation yet. However, if Scott's Long Emulation is broken, this
is the patch that will exercise it. :)

-Ray

On Wed, Mar 26, 2008 at 12:00 PM, Freeland Abbott <fab...@google.com> wrote:

This is a corner of a corner case; I'm not sure that I care.  That said, if I'm forced to pick one, I think the explicit assignment seems cleaner; I don't think the performance cost is significant, and it makes the code much more "as expected" rather than "read the fine print."

On Wed, Mar 26, 2008 at 2:49 PM, Ray Cromwell <cromw...@gmail.com> wrote:

I just looked at my working previous patch and it has the same "bug". On testing, it's not a bug because of implicit conversions of NaN/INF.

That is, if x = Infinity

mantissa = (x % 1) * Math.pow(2,52)

(x % 1) = NaN
NaN * Math.pow(2,52) = NaN

But NaN & 0xANYBITS = 0, that is, Infinity and NaN get coerced to zero, so although I could add an explicit "x = 0" in that case, it would be a no-op.

Should I force the variable to be 0, or just comment on the implicit nature of it working.

I did rename NUMBER_BITSIZE back to SIZE, I'll attach a new patch based on your response.

-Ray

On Wed, Mar 26, 2008 at 7:52 AM, Freeland Abbott <gwt.team...@gmail.com> wrote:

Glad I hadn't gotten to it, then!

On Wed, Mar 26, 2008 at 2:54 AM, Ray Cromwell <cromw...@gmail.com> wrote:

Freeland,
  I just reviewed my code again, I really fubared that patch. The case where d == Infinity fails, because of a missing line assigning absV = 0 (so that the mantissa becomes 0 , which is how Infinity is represented) Also, NUMBER_SIZE should be just SIZE which is the official JDK variable, I refactored this by accident. I'll send you a fixed version tommorow, hopefully the last one unless you find some issues. :)

-Ray

On Tue, Mar 25, 2008 at 8:33 AM, Ray Cromwell <cromw...@gmail.com> wrote:

Yeah, but I made the mistake of changing the subject line, here it is: http://groups.google.com/group/Google-Web-Toolkit-Contributors/browse_thread/thread/419e9675df541dfd

On Tue, Mar 25, 2008 at 8:12 AM, Freeland Abbott <gwt.team...@gmail.com> wrote:

Did you attach the patch?

[Freeland Abbott]

Submit at r2421.

[Ray Cromwell]

Coolio. Freeland, thanks for the review, I know you guys are probably
busy right now trying to push out 1.5 releases especially with Google
I/O coming up, and this patch isn't that critical, so much
appreciated.

-Ray

On Wed, Apr 9, 2008 at 11:37 AM, Freeland Abbott

[Freeland Abbott]

Turns out we leaked some Java 1.6 identifiers, though (Double.MIN_NORMAL isn't defined in hosted mode tests). :-/