Here is a program demonstrating the round-off error associated with the
64-bit binary representation of the decimal number 0.1:
-------------------------------------------------------------
x=0.1
for n=1:7
s=0;
for i=1:10^n
s=s+x;
end;
expected = x * 10^n;
error = s - expected;
fprintf('Round-off error from summing %g 10^%d times is %.2e\n',x,n,error);
end;
---------------------------------------------------------------
Note that summing 0.1 only 10 times already leads to a sizable
round-off error! As an exercise, change x=0.1 to x=0.5, and you will
see that there is no error whatsoever! This is because 0.5 is a
negative power of 2.
As a little algebra exercise, I want to prove the binary representation
of the decimal number 0.1 that I have written in class:
0.1 (decimal) = 0.00011001100110011001100110011..... (binary)
This means that 0.1 = (1/2^4 + 1/2^5) + (1/2^8 + 1/2^9) + (1/2^12 +
1/2^13) +...
Note that this can be written as
0.1 = (1/2^4 + 1/2^5) + (1/2^4 + 1/2^5) / 24 + (1/2^4 + 1/2^5) /
28 +...
factoring out the first bracketed sum yields
0.1 = (1/2^4 + 1/2^5) (1 + 2^(-4) + 2^(-8) + 2^(-12) +... )
Now, the infinite series in the second set of brackets is a geometric
series with a ratio of 2^(-4), which can be easily summed, so:
0.1 = (1/2^4 + 1/2^5) / (1 - 1/2^4)
Multiplying the numerator and the denominator by 24, we get
0.1 = (1 + 1/2) / (2^4 - 1) = 1.5 / 15 = 0.1 QED!
------------------------------------------------------------------------
Here is also the second program - the function computing the cosine
using the Taylor series, and which checks for possible NaN and Inf
occurring during execution. Note the "elseif" construct: lookup "help
if" to see the full syntax of the "if" command. Now, try calling this
function with big (x > 70) and small values of the argument
function s = myCos(x)
s= 0;
tolerance = 1e-14;
for n = 0:1000
an = (-1)^n * x^(2*n) / factorial(2*n);
s = s + an;
if ( abs(an) < tolerance )
fprintf('myCos(%g)=%g, error=%.2e, %d terms used\n', x, s, s-cos(x), n);
break;
elseif ( abs(an) == Inf )
fprintf('Term number %d is too large: terminating\n', n);
break;
elseif ( isnan(an) )
fprintf('Term number %d is NaN: terminating\n', n);
break;
end;
end;
-------------------------------------------------
There is really no good reason we cannot use the condition (an == NaN)
instead of "isnan(an)" in above function - it's just a quirk of Matlab.