Search code examples
c++doublelisprational-numberieee

double precision and peridoc decimals


I'm working on a lisp interpreter and implemented rational numbers. I thought they have the advantage over doubles to be able to represent numbers like 1/3. I did some calculations to compare the results. I was surprised by the results

with doubles

(* 3.0 (/ 1.0 3.0)) -> 1
(* 3.0 (/ 4.0 3.0)) -> 4
(* 81.0 (/ 1.0 81.0)) -> 1

with ratios:

(* 3 (/ 1 3)) -> 1
(* 3 (/ 4 3)) -> 4
(* 81 (/ 1 81)) -> 1

Why are the results of the floating point operations exact? There must be a loss of precision. doubles cannot store an infinit number of digits. Or do I miss something?

I did a quick test with a small C-Application. Same result.

#include <stdio.h>

int main()   
{  
    double a1 = 1, b1 = 3;  
    double a2 = 1, b2 = 81;  

    printf("result : %f\n", a1 / b1 * b1); 
    printf("result : %f\n", a2 / b2 * b2);

    return 0;  
}  

Output is:

result : 1.000000
result : 1.000000

MFG

Martin


Solution

  • For the first case, the exact result of the multiply is half way between 1.0 and the largest double that is less than 1.0. Under IEEE 754 round-to-nearest rules, half way numbers are rounded to even, in this case to 1.0. In effect, the rounding of the result of the multiply undid the error introduced by rounding of the division result.

    This Java program illustrates what is happening. The conversions to BigDecimal and the BigDecimal arithmetic operations are all exact:

    import java.math.BigDecimal;
    
    public class Test {
      public static void main(String[] args) {
        double a1 = 1, b1 = 3;
    
        System.out.println("Final Result: " + ((a1 / b1) * b1));
        BigDecimal divResult = new BigDecimal(a1 / b1);
        System.out.println("Division Result: " + divResult);
        BigDecimal multiplyResult = divResult.multiply(BigDecimal.valueOf(3));
        System.out.println("Multiply Result: " + multiplyResult);
        System.out.println("Error rounding up to 1.0: "
            + BigDecimal.valueOf(1).subtract(multiplyResult));
        BigDecimal nextDown = new BigDecimal(Math.nextAfter(1.0, 0));
        System.out.println("Next double down from 1.0: " + nextDown);
        System.out.println("Error rounding down: "
            + multiplyResult.subtract(nextDown));
      }
    }
    

    The output is:

    Final Result: 1.0
    Division Result: 0.333333333333333314829616256247390992939472198486328125
    Multiply Result: 0.999999999999999944488848768742172978818416595458984375
    Error rounding up to 1.0: 5.5511151231257827021181583404541015625E-17
    Next double down from 1.0: 0.99999999999999988897769753748434595763683319091796875
    Error rounding down: 5.5511151231257827021181583404541015625E-17
    

    The output for the second, similar, case is:

    Final Result: 1.0
    Division Result: 0.012345679012345678327022824305458925664424896240234375
    Multiply Result: 0.9999999999999999444888487687421729788184165954589843750
    Error rounding up to 1.0: 5.55111512312578270211815834045410156250E-17
    Next double down from 1.0: 0.99999999999999988897769753748434595763683319091796875
    Error rounding down: 5.55111512312578270211815834045410156250E-17
    

    This program illustrates a situation in which rounding error can accumulate:

    import java.math.BigDecimal;
    
    public class Test {
      public static void main(String[] args) {
        double tenth = 0.1;
        double sum = 0;
        for (int i = 0; i < 10; i++) {
          sum += tenth;
        }
        System.out.println("Sum: " + new BigDecimal(sum));
        System.out.println("Product: " + new BigDecimal(10.0 * tenth));
      }
    }
    

    Output:

    Sum: 0.99999999999999988897769753748434595763683319091796875
    Product: 1
    

    Multiplying by 10 rounds to 1.0. Doing the same multiplication by repeated addition does not get the exact answer.