Re: The old round off problem?

David Treadwell Sun, 05 Mar 2006 14:54:20 -0800

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

-- 
http://mail.python.org/mailman/listinfo/python-list

Re: The old round off problem?

Reply via email to