(math) floating point rounding error on win32

Name: laC46010 Date: 10/30/97



The following example floating-point multiplication
produces different results on solaris and win32:

---------------------------------------------------------
class t {
  public static void main (String[] args) {
    double x = Double.longBitsToDouble(0x1e7ee00000000000L);
    double y = Double.longBitsToDouble(0x2180101010101010L);
    double z = x*y;
    System.out.println(Long.toBinaryString(Double.doubleToLongBits(z)));
  }
}
---------------------------------------------------------
1111011111110111111101111111011111110111111101111111 - solaris (sparc, x86)
1111011111110111111101111111011111110111111110000000 - win32

This is reproduced with different version of JDK (1.0.2 - 1.2)
It's interesting to note that resulting value is denormalized
floating-point number i.e. its exponent is equal to 0 (= -1023).

The result can be reproduced with other numbers:

class t2 {
  static String dup(String s, int n) {
   return n == 0 ? "" : s + dup(s, n-1);
  }
  static String bin(double d) {
   String s = Long.toBinaryString(Double.doubleToLongBits(d));
        return (dup("0", 64 - s.length()) + s);
  }

  static void mul(String sx, String sy, String sr) {
    double x = Double.longBitsToDouble(Long.parseLong(sx, 2));
    double y = Double.longBitsToDouble(Long.parseLong(sy, 2));
    double w = Double.longBitsToDouble(Long.parseLong(sr, 2));
    double z = x*y;
    if (z != w) {
     throw new RuntimeException("Rounding error in mantissa"
     + "\n Op1: " + bin(x)
     + "\n Op2: " + bin(y)
     + "\n Expected: " + sr
     + "\n Actual: " + bin(z) );
    }
  }
  public static void main (String[] args) {
    mul("0001111001111001000000000000000000000000000000000000000000000000"
       ,"0010000110000001000100010001000100010001000100010001000100010001"
       ,"0000000000001101010101010101010101010101010101010101010101010101");
  }
}

Output:

java.lang.RuntimeException: Rounding error in mantissa
      Op1: 0001111001111001000000000000000000000000000000000000000000000000
      Op2: 0010000110000001000100010001000100010001000100010001000100010001
 Expected: 0000000000001101010101010101010101010101010101010101010101010101
   Actual: 0000000000001101010101010101010101010101010101010101010101010110
        at t2.mul(t2.java:16)
        at t2.main(t2.java:24)



See also original examples below.


/*
 * @test
 * @bug <putNumberHere>
 * @summary Floating-point rounding error on win32
 *
 * Because test implemented with compile-time constants, bug also
 * "appears" on other platforms if running win32-compiled class.
 * (Same result with runtime variables - see original submission below.)
 *
*/
class RoundingFPBug {

  public static void main (String[] args) {
    
    final double roundError =
    ((247.0/256.0)
    / ( ((double) (1L << 50) * (double) (1L << 57))
     * ((double) (1L << 50) * (double) (1L << 57))
     * ((double) (1L << 50) * (double) (1L << 57))
     * ((double) (1L << 50) * (double) (1L << 57))
     * ((double) (1L << 50) * (double) (1L << 57))
     ))
   *
    ((double) 0x1010101010101L
    / ( ((double) (1L << 50) * (double) (1L << 57))
     * ((double) (1L << 50) * (double) (1L << 57))
     * ((double) (1L << 50) * (double) (1L << 57))
     * ((double) (1L << 50) * (double) (1L << 57))
     * ((double) (1L << 50) * (double) (1L << 57))
     )) ;
    final long mantissa = Double.doubleToLongBits(roundError) & 0xfffffffffffffL;
    final long expected = 0xF7F7F7F7F7F7FL;
    // 0b1111011111110111111101111111011111110111111101111111L;
  
    if (mantissa != expected ) {
     throw new RuntimeException("Rounding error in mantissa"
     + "\n Expected: " + Long.toBinaryString(expected)
     + "\n Actual: " + Long.toBinaryString(mantissa) );
    }
  }
}

/************************************ start original submission ***********
class RoundingTest {

  /* // restore comment here
    This is supposed to illustrate the effects of double rounding.
    Because Intel's "double precision" emulation does not reduce the
    size of the "double" mantissa when the exponent strays into the
    denormalized range of IEEE floating point numbers, two roundings
    occur. In this example, the first rounding occurs with a 53-bit
    mantissa, the second with a 52-bit denormalized mantissa. The
    slightly doctored output of this program on an intel machine is
    shown below.

    Note that the normalized-representation numbers (exponents > -1023)
    have an implicit leading one bit. Also note that IEEE rounding
    chooses the even (trailing zero mantissa) outcome in the event of
    a tie.

    Sgn=0, exp= -536, man=0b1110111000000000000000000000000000000000000000000000
    Sgn=0, exp= -487, man=0b0000000100000001000000010000000100000001000000010000
    Sgn=0, exp=-1023, man=0b1111011111110111111101111111011111110111111110000000
    EXACT MANTISSA 11110111111101111111011111110111111101111111011111110111
                            12345678901234567890123456789012345678901234567890123456
    53-bit rounding 11110111111101111111011111110111111101111111011111111
    2nd (52-bit) rounding 1111011111110111111101111111011111110111111110000000
    Single 52-bit rounding 1111011111110111111101111111011111110111111101111111

    Notice that the two 52-bit roundings differ. That is the problem.
    * / // <--------------- [restore comment here if using original submission]

  static void printDoubleRep(double d) {
    long l = Double.doubleToLongBits(d);
    int sign = (int) (l >> 63) & 1;
    int exponent = (int) (l >> 52) & 0x7ff;
    long mantissa = l & 0xfffffffffffffL;
    System.out.println("Sgn=" + sign +
", exp=" + Integer.toString(exponent-1023) +
", man=0b" + Long.toBinaryString(mantissa));
  }

  public static void main (String[] args) {
    
    double two_to_the_107 = (double) (1L << 50) * (double) (1L << 57);
    double two_to_the_535 =
      two_to_the_107 * two_to_the_107 *
      two_to_the_107 * two_to_the_107 * two_to_the_107;
    
    double factor1 = (247.0/256.0) / two_to_the_535;

    long mantissa = 0x1010101010101L; // 123456789012 = 49 bits total

    double factor2 = (double) mantissa / two_to_the_535;

    double product = factor1 * factor2;
    printDoubleRep(factor1);
    printDoubleRep(factor2);
    printDoubleRep(product);
  }
}


Hook 5(hook5): test

======================================================================