%% Demonstration of rounding error issues
%
% Example from the book N. Higham, Accuracy and stability of Numerical
% Algorithms, SIAM, 2002
%
% We are interested in evaluating the function
%
% $$ f(x)=\frac{1 - \cos x}{x^2}. $$
%
x = 1.2e-8;
c = cos(x);

fprintf('x = %20.16e \n',x);
fprintf('cos(x) = %20.16e (.999999999999999928)\n',c);
fprintf('1-cos(x) = %20.16e (.7200000000000000e-16)\n',1-c);
fprintf('(1-cos(x))/x^2 = %20.16e \n',(1-c)/x^2);

%%
% the result is calculated inaccurately and then amplified by
% $x^{-2}$.
%

%% Numerically stable evaluation
%
% Using
% 
% $$ \cos(x) = 1 - 2 \sin^2\frac{x}{2}, $$
%
% $f(x)$ can be evaluated in a numerically stable manner
%
y = x/2;
f = 0.5*(sin(y)/y)^2;


fprintf('0.5*(sin(x/2)/(x/2))^2 = %20.16e \n',f);
