% This concerns cancellation errors.  And is related to the HW problem
% Chapter 1, problem 27.  
%
% Natasha asked, "How can you tell if you might have cancellation errors
% when you don't know exactly what the correct answer is?  I.e., if you
% were trying to calculate cos(2pi) and you didn't already know that
% cos(2pi) = 1..."  
%
% The answer is, "Well, you have a program that takes x and gives out y.
% You want to know what happens if you perturb x slightly and then give
% it to the program.  If small perturbations are magnified hugely then
% you may have cancellation errors.  In any case, if small perturbations
% are magnified hugely then you definitely have an ill-conditioned
% problem (this is how you define "ill-conditioned") and you need to
% find some other way of computing y."
%
% In the following, we take x to be 2pi, 4pi, 6pi, 8pi.  We first
% estimate how far in the power series we should go before the tail
% becomes of the size of round-off.  (eps is round-off).  We then
% compute the partial sum, adding up in reverse order to minimize
% round-off error.

x = 2*pi;
%x = 4*pi;
%x = 6*pi;
%x = 8*pi;
tol = 10*eps;

test = 10*tol;
n = 1;
while test > tol,
% want to see how far out we have to go until the tail of the
% series becomes small:
  test = abs(x^(2*n)/gamma(2*n+1));
  n = n+1;
end

my_sum = 0;
for i=0:n,
% do the sum from smallest to largest:
  ii = n-i;
  my_sum = my_sum + (-1)^ii*x^(2*ii)/gamma(2*ii+1);
end

% okay, so that gives us an approximation of cos(x).  Now we want to see
% what happens with small changes in x.

for j=1:99
  xx(j) = x + (j-50)*100*eps;
  my_sum(j) = 0;
  for i=0:n,
% do the sum from smallest to largest:
    ii = n-i;
    my_sum(j) = my_sum(j) + (-1)^ii*xx(j)^(2*ii)/gamma(2*ii+1);
  end
end

n
max(my_sum)-min(my_sum)
(xx(99)-xx(1))
(max(my_sum)-min(my_sum))/(xx(99)-xx(1))


% for x=2pi, n = 21, (max(my_sum)-min(my_sum))/(xx(99)-xx(1)) = 2.0e-2
% for x=4pi, n = 31, (max(my_sum)-min(my_sum))/(xx(99)-xx(1)) = 1.0e+1
% for x=6pi, n = 40, (max(my_sum)-min(my_sum))/(xx(99)-xx(1)) = 4.3e+3
% for x=8pi, n = 49, (max(my_sum)-min(my_sum))/(xx(99)-xx(1)) = 2.1e+6
%
% so for x=8pi, we have xx ranging over [8pi-1e-12,8pi+1e-12]
% and my_sum ranging over [1-2.3e-6,1+2.3e-6].  So a small perturbation
% in the input has a huge effect on the output.  This is a signal that
% there might be cancellation errors.  Which there are.