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.