The old round off problem?

[David Treadwell]

On Mar 5, 2006, at 1:01 AM, sam wrote:

>
> David Treadwell wrote:
>> exp(x) is implemented by:
>>
>> 1.  reducing x into the range |r| <=  0.5 * ln(2), such that x = k *
>> ln(2) + r
>> 2.  approximating exp(r) with a fifth-order polynomial,
>> 3.  re-scaling by multiplying by 2^k: exp(x) = 2^k * exp(r)
>>
>> sinh(x) is mathematically ( exp(x) - exp(-x) )/2 but it's not
>> calculated directly that way.
>>
>> My brain hurts at this point; it's late. I'll have another go at
>> hunting down the errors tomorrow. In the interim, does anybody out
>> there already know the answer?
>>
>> :--David
> I tried the exp(x) equivalent of sinh(x) and cosh(x) with the same
> results.
> I think I will write my own or steal the Taylor(f(x),x,n) function in
> both C++, and python.

After dreaming about this last night, there is a pretty  
straightforward answer.

A Python or C double precision real will deliver around 16 or 17  
meaningful decimal digits.

from math import *

for x_int in range(20):
         x_real = 1.0 + x_int
         sinh_x = sinh(x_real)
         cosh_x = cosh(x_real)
         c2s2_x = (cosh_x + sinh_x)*(cosh_x - sinh_x)
         print 'x:% *.2f,  cosh: % .4e,  sinh: % .4e,  c2-s2: % .4e'% 
(6, x_real, cosh_x, sinh_x,c2s2_x-1.0)


x:  1.00,  cosh:  1.5431e+00,  sinh:  1.1752e+00,  c2-s2:  0.0000e+00
x:  2.00,  cosh:  3.7622e+00,  sinh:  3.6269e+00,  c2-s2: -2.3315e-15
x:  3.00,  cosh:  1.0068e+01,  sinh:  1.0018e+01,  c2-s2: -2.4314e-14
x:  4.00,  cosh:  2.7308e+01,  sinh:  2.7290e+01,  c2-s2:  1.4766e-13
x:  5.00,  cosh:  7.4210e+01,  sinh:  7.4203e+01,  c2-s2:  1.4677e-12
x:  6.00,  cosh:  2.0172e+02,  sinh:  2.0171e+02,  c2-s2:  1.9333e-12
x:  7.00,  cosh:  5.4832e+02,  sinh:  5.4832e+02,  c2-s2: -3.6329e-11
x:  8.00,  cosh:  1.4905e+03,  sinh:  1.4905e+03,  c2-s2: -1.7109e-10
x:  9.00,  cosh:  4.0515e+03,  sinh:  4.0515e+03,  c2-s2:  3.1331e-09
x: 10.00,  cosh:  1.1013e+04,  sinh:  1.1013e+04,  c2-s2:  2.6562e-08
x: 11.00,  cosh:  2.9937e+04,  sinh:  2.9937e+04,  c2-s2: -1.2103e-07
x: 12.00,  cosh:  8.1377e+04,  sinh:  8.1377e+04,  c2-s2:  2.2313e-06
x: 13.00,  cosh:  2.2121e+05,  sinh:  2.2121e+05,  c2-s2: -4.2111e-06
x: 14.00,  cosh:  6.0130e+05,  sinh:  6.0130e+05,  c2-s2: -1.0882e-04
x: 15.00,  cosh:  1.6345e+06,  sinh:  1.6345e+06,  c2-s2:  1.2143e-04
x: 16.00,  cosh:  4.4431e+06,  sinh:  4.4431e+06,  c2-s2: -6.8998e-03
x: 17.00,  cosh:  1.2077e+07,  sinh:  1.2077e+07,  c2-s2: -1.0174e-02
x: 18.00,  cosh:  3.2830e+07,  sinh:  3.2830e+07,  c2-s2: -2.1590e-02

Consider just the characteristic, int(log10(), of cosh(x) and sinh 
(x). In the table, it varies from 0 to +7. Because it ends up getting  
squared, the characteristic of cosh(x)**2 varies from 0 to 15.

If you subtract that number from 16, you get roughly the precision of  
the result: The precision of the final answer for cosh(x)**2 - sinh(x) 
**2 - 1.0 is roughly 17 - 2 * int(log10(cosh(x)).

So, this isn't a rounding problem. It is a lack of precision problem.  
As I suspected last night, subtracting two very large and almost  
equal numbers from each other eliminates precision. It just doesn't  
look like it once the internal hex is converted to decimal.

Because the problem is mathematically correct, there must be a better  
way to state it so that you don't end up with this problem.

:--David
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/python-list/attachments/20060305/b58501b8/attachment.html>