#!/usr/bin/octave -q

% Jonathan R. Senning <senning@gordon.edu>
% Gordon College
% January 27, 1999
%
% This program demonstrates the lost of precision that occurs when two
% numbers very close to each other are subtracted.  This example uses
% x and sin(x) for a small value of x.
%
% When the value of x - sin(x) is needed for small values of x, a good
% approach is to subtract the Taylor series for sin(x) about 0 from x.
% This program computes the value of the truncated Taylor series (n terms)
% in two different ways; as a nested polynomial and as a regular polynomial.
%
% Here is the output for the case where n = 5 and x = 1e-6:
%
%	x      =  9.9999999999999995e-07
%	sin(x) =  9.9999999999983330e-07
%
%	                   Nested Taylor Series    Regular Taylor Series 
%	                 -----------------------  -----------------------
%	Series solution:  1.6666666666665834e-19   1.6666666666665832e-19
%	Direct solution:  1.6665373237228359e-19   1.6665373237228359e-19
%	Difference     :  1.2934294374749207e-23   1.2934294374725133e-23

n = 5;
x = 1e-6;

xx = x * x;
sinx = sin( x );
xsinx = x - sinx;

fprintf( 'x      = %23.16e\n', x );
fprintf( 'sin(x) = %23.16e\n', sinx );

fprintf( '\n' );
fprintf( '                   Nested Taylor Series   Regular Taylor Series \n' );
fprintf( '                 ----------------------- -----------------------\n' );

% Compute x - sin(x) using nested multiplication of Taylor series

s = 1.0;
for i = n + 2 - [ 2 : n ]
    s = 1 - xx * s / ( (2*i) * (2*i+1) );
end
s = x * xx * s / 6.0;

% Compute x - sin(x) using regular evaluation of Taylor series
    
t = 0;
for i = 1 : n
    p = 1.0;
    for j = 1 : 2*i+1
	p = p * ( x / j );
    end
    if ( 2 * floor( i / 2 ) == i )
	t = t - p;
    else
	t = t + p;
    end
end

fprintf( 'Series solution: %23.16e %23.16e\n', s, t );
fprintf( 'Direct solution: %23.16e %23.16e\n', xsinx, xsinx );
fprintf( 'Difference     : %23.16e %23.16e\n', s - xsinx, t - xsinx );