% This script can be used to explain how a divergent sum can converge in
% finite precision arithmetic--in this case in 2-digit decimal
% floating-point arithmetic.  In the first run, the sum becomes
% stationary when it reaches a value of 3.9.  The second run prints out
% enough additional information to show what is going on.  The teacher
% can supply running comments.

format short e
format compact
SetFlapRlevel(2);
disp(' ')
disp('Summing 1 + 1/2 + 1/3 + 1/4 + ...  2-digits FP arithmetic.')
disp(' ')
sum = flap(0);
for n=1:25
   sum = sum + 1/n;
   disp(sum.d+10*eps)
end
pause
disp(' ')
disp('Same sum with more information')
disp(' ')
disp('   sum          fl(1/n)      sum+fl(1/n)  fl[sum+fl(1/n)]')
disp(' ')
pause
sum = flap(0);
for n=1:25
   term = flap(1/n);
   nextsum = sum + term;
   disp([sum.d, term.d, sum.d+term.d, nextsum.d+10*eps])
   sum = nextsum;
   pause
end