0

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

Martin Fehrs
  • 792
  • 4
  • 13
  • Try `%.20f`. Double precision values are accurate to about 15 decimal places, even when not exact. – Mike Seymour Jan 08 '15 at 12:05
  • [Read this first please](http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html). – πάντα ῥεῖ Jan 08 '15 at 12:06
  • Could it be the compiler noticing the division and multiplication and thus making it void (or rearranging the order of the operation) for you? – David Jan 08 '15 at 12:11
  • @David: strange as it sounds, when I check the assembly output of MinGW g++ in Windows it really does the computations, yet produces exact result (checked by subtraction). thus, the end result is evidently rounded up instead of down. but then i don't know any reason why it should be rounded down. – Cheers and hth. - Alf Jan 08 '15 at 12:16
  • using %.20f doesn't change the results. It's really exact! Seems like the result is coincidentally rounded to the correct result. What I'm looking for is an example where doubles fail. – Martin Fehrs Jan 08 '15 at 12:19
  • Is the Lisp code and the associated results in your question from *your* interpreter, or from some 'standard' one? Even though "doubles cannot store an infinit number of digits", no one has said that Lisp implementations have to use doubles to represent decimal numbers. – Joshua Taylor Jan 08 '15 at 13:08

1 Answers1

2

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.

Patricia Shanahan
  • 25,849
  • 4
  • 38
  • 75
  • Thank you very much for this detailed explanation. I understand now what's happening. But I'm still searching for a negative example where reproducing the original number fails. – Martin Fehrs Jan 08 '15 at 12:50
  • @user3624086 What are your constraints for the counter-example? For example, can the divisor be smaller than 1? – Patricia Shanahan Jan 08 '15 at 13:28
  • It should be compareable with a ratio calculation to visualise the differences. So, no! The divisor should be a whole number. Something like (* 3 (/ 1 3)) vs (* 3.0 (/ 1.0 3.0)). (/ 1 3) produces as rational result of 1/3 in my interpreter. You can test it on remote-lisp.spdns.de. But it's not very comfortable at the moment. Well, of course it can be more complex as long as it's still understandable. I had no luck so far. Lets relax the contraints. The result should be any whole number. – Martin Fehrs Jan 08 '15 at 13:40
  • @user3624086 The easiest way to demonstrate rounding error is by repeated addition, rather than directly multiplying. Would you consider that? – Patricia Shanahan Jan 08 '15 at 13:45
  • Of course. The kind of operation doesn't matter. But how much nubmers I have to add before the error isn't rounded away? Sadly I don't have a loop construct implemented at the moment. Maybe I should focus on this. – Martin Fehrs Jan 08 '15 at 13:48
  • @user3624086 I've edited in an illustration of accumulated rounding error. Another way of getting rounding error quickly is to add two numbers of very different magnitude, and then subtract the larger or something similar to it. – Patricia Shanahan Jan 08 '15 at 13:53
  • Well, your example works if I set the visible precision very high. If i set the visible prevision to 15 which is the supported precission of a double on my system, the error is always rounded correctly away and I get the exact result. Seems like doubles are really good at handling periodic numbers. Far better than I thought. – Martin Fehrs Jan 08 '15 at 15:12
  • Finally I did it. When calulcating the sum of 21 instances of 0.1 I get a noticable error. Yay! – Martin Fehrs Jan 08 '15 at 15:31