# A Glimpse into Floating-Point Accuracy 31

Posted by **Loren Shure**,

There are frequent posts on the MATLAB newsgroup as well as lots of questions posed to Technical Support about the floating point accuracy of MATLAB. Many people think they have found a bug when some seemingly simple arithmetic doesn't give the intuitive answer. The issue that arises has to do with the finite number of bits that a computational device uses to store and operate on values.

### Contents

### Calculator Demo

Have you ever impatiently waited for an elevator and continued pushing the `UP` button in the hopes that the elevator will arrive faster? If so, you may also have impatiently pushed buttons repetitively
on a calculator (assuming you know what a calculator is!). It's not just MATLAB that computes values in a way that seems to
defy our expectations sometimes.

```
clear
format long
```

Here's what I tried. I put the number 7 into my calculator. I think it calculates in single precision, but I'm not sure. It shows 7 digits after the decimal point when I take the square root, and doing it several times, I get

sqrts7inSingle = single([7 2.6457513 1.6265765 1.275373])

sqrts7inSingle = 7.0000000 2.6457512 1.6265765 1.2753730

Next, I take the final number and square it 3 times, yielding

squaringSqrt7inSingle = single([1.275373 1.6265762 2.6457501 6.99999935])

squaringSqrt7inSingle = 1.2753730 1.6265762 2.6457500 6.9999995

Compare the difference between the original and final values.

sqrts7inSingle(1)-squaringSqrt7inSingle(end)

ans = 4.7683716e-007

This number turns out to be very interesting in this context. Compare it to the distance between a single precision value of 7 and its next closest value, also know as the floating-point relative accuracy.

eps(single(7))

ans = 4.7683716e-007

We get the same value! So, we can demonstrate the effects of finite storage size for values on a handheld calculator. Let's try the equivalent in MATLAB now.

### Accuracy in MATLAB

Let's do a similar experiment in MATLAB using single precision so we have a similar situation to the one when I used the calculator. In this case, I am going to keep taking the square root until the value reaches 1.

num = single(7); count = 0; tol = eps('single'); dnum = sqrt(num); while abs(dnum-1)>=tol count = count+1; dnum = sqrt(dnum); end count

count = 23

So, it took 23 iterations before we reached the number 1. But now we're in a funny situation. There is no way to square exactly 1 and ever reach our original value of 7.

### Using Double Precision

Let's repeat the same experiment from the calculator in MATLAB using double precision.

num = 7; sqrtd = sqrt(num); sqrtd(2) = sqrt(sqrtd(1)); sqrtd(3) = sqrt(sqrtd(2)) prods(1) = sqrtd(end); prods(2) = prods(1)^2; prods(3) = prods(2)^2 finalNum = prods(3)^2

sqrtd = 2.64575131106459 1.62657656169779 1.27537310685845 prods = 1.27537310685845 1.62657656169779 2.64575131106459 finalNum = 7.00000000000001

Find the difference between `finalNum` and `num` and compare to the proper relative floating-point accuracy:

diffNum = finalNum-num acc = eps(num)*6 diffNum-acc

diffNum = 5.329070518200751e-015 acc = 5.329070518200751e-015 ans = 0

Why did I use 6 when I calculated `acc`? Because I performed 6 floating point operations to get the my final answer, 3 `sqrt` applications, followed by squaring the answers 3 times.

### Typical MATLAB Pitfall

The most common MATLAB pitfall I run across is when users check equality of values generated from the `:` operator. The code starts innocently enough. I've actually turned the logic around here a bit, so we will print out values in the
loop where the loop counter may not be exactly what the user expects.

format short for ind = 0:.1:1; if ind ~= fix(10*ind)/10 disp(ind - fix(10*ind)/10) end end

5.5511e-017 1.1102e-016 1.1102e-016

And then the question is, why do *any* of the values print out? It's because computers can't represent all numbers exactly given a fixed storage size such as double
precision.

### Conclusions

This was a very simplified explanation about floating point, but one, I hope, that is easily understood. Let me know.

### References

Here are some pointers to more resources.

Published with MATLAB® 7.2

**Category:**- Common Errors

### Note

Comments are closed.

## 31 CommentsOldest to Newest

**1**of 31

**2**of 31

**3**of 31

**4**of 31

**5**of 31

**6**of 31

**7**of 31

**8**of 31

#include <stdio.h> int main(int nargin, char *argv[]) { double a = 0.1; double b = 0.3; if ( (a+a+a) != b ) { fprintf(stdout, "Not equal\n"); } else { fprintf(stdout, "Equal\n"); } return(0); }

**9**of 31

#include <iostream> int main() { double a = 0.1; double b = 0.3; if ( (a+a+a) != b ) { std::cout << "Not equal\n"; } else { std::cout << "Equal\n"; } return(0); }

**10**of 31

**11**of 31

>> x1 = hex2dec('8D6F7897') x1 = 2.3729e+009 >> x2 = hex2dec('41C64E6D') x2 = 1.1035e+009 >> x3 = hex2dec('6073') x3 = 24691 >> x4 = x1.*x2+x3; >> x5 = uint64(x4) x5 = 2456dd29cafeba00--Loren

**12**of 31

**13**of 31

**14**of 31

**15**of 31

**16**of 31

**17**of 31

**18**of 31

**19**of 31

**20**of 31

```
% Tell MATLAB to display the numbers in hex format
format hex
% Look at the number 1
one = 1
% This is the next largest number
% The last bit of the mantissa has been incremented
nextLargest = one+eps(one)
% This is the next smallest number
% We decremented the last bit of the mantissa
% which required borrowing from the exponent
nextSmallest = one - (eps(one)/2)
% Is nextSmallest actually less than 1?
isSmaller = nextSmallest < one
% Reset the format back to the default
format
```

Cleve's newsletter article, linked at the end of Loren's blog posting, describes the hexadecimal format displayed by FORMAT HEX. **21**of 31

**22**of 31

**23**of 31

*The six coefficients are close to, but not exactly equal to, the power series coefficients 1/3!, 1/5!, …, 1/13! . They minimize the maximum relative error, | (sin() - p())/sin() |, over the interval. Six terms are enough to make this approximation error less than 2^(-52), which is the roundoff error involved when all the terms, and the sum, are less than one.*That article is titled "The Tetragamma Function and Numerical Craftsmanship" -- and there's definitely some craftmanship involved in the implementation of the trig functions.

**24**of 31

**25**of 31

**26**of 31

**27**of 31

**28**of 31

**29**of 31

**30**of 31

**31**of 31

## Recent Comments