#!/usr/bin/env python

# Jonathan R. Senning <jonathan.senning@gordon.edu>
# Gordon College
# January 27, 1999
# Converted to Python May 2008
#
# 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

from pylab import *

n = 5
x = 1e-6

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

print 'x      = %23.16e' % x
print 'sin(x) = %23.16e' % sinx

print ''
print '                   Nested Taylor Series    Regular Taylor Series'
print '                 -----------------------  -----------------------'

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

s = 1.0
for i in n + 2 - arange( 2 , n + 1 ):
    s = 1 - xx * s / ( ( 2 * i ) * ( 2 * i + 1 ) )

s = x * xx * s / 6.0

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

print 'Series solution: %23.16e  %23.16e' % ( s, t )
print 'Direct solution: %23.16e  %23.16e' % ( xsinx, xsinx )
print 'Difference     : %23.16e  %23.16e' % ( s - xsinx, t - xsinx )