The old round off problem?

[sam]

Hello all, I am taking a class in scientific programming at the local
college. My problem is that the following difference produces round off
errors as the value of x increases. For x >= 19 the diference goes to
zero.I understand the problem, but am curious as to whether their
exists a solution. I have tried various graphing programs,and they all
exihibit this problem.

thanks in advance Sam Schulenburg
f(x) = cosh^2(x) - sinh^2(x) = 1
[color=blue][color=green][color=darkred]
>>> from math import *
>>> for x in range(20):[/color][/color][/color]
print "x= %2d Sinh^2(x) = %20.3f f(x) = %2.10f
"%(x,pow(cosh(x ),2),pow(cosh(x ),2)- pow(sinh(x),2))


x= 0 Sinh^2(x) = 1.000 f(x) = 1.0000000000
x= 1 Sinh^2(x) = 2.381 f(x) = 1.0000000000
x= 2 Sinh^2(x) = 14.154 f(x) = 1.0000000000
x= 3 Sinh^2(x) = 101.358 f(x) = 1.0000000000
x= 4 Sinh^2(x) = 745.740 f(x) = 1.0000000000
x= 5 Sinh^2(x) = 5507.116 f(x) = 1.0000000000
x= 6 Sinh^2(x) = 40689.198 f(x) = 1.0000000000
x= 7 Sinh^2(x) = 300651.571 f(x) = 0.9999999999
x= 8 Sinh^2(x) = 2221528.130 f(x) = 1.0000000000
x= 9 Sinh^2(x) = 16414992.784 f(x) = 1.0000000037
x= 10 Sinh^2(x) = 121291299.352 f(x) = 1.0000000298
x= 11 Sinh^2(x) = 896228212.033 f(x) = 0.9999998808
x= 12 Sinh^2(x) = 6622280532.961 f(x) = 1.0000019073
x= 13 Sinh^2(x) = 48932402357.710 f(x) = 0.9999923706
x= 14 Sinh^2(x) = 361564266073.36 9 f(x) = 0.9999389648
x= 15 Sinh^2(x) = 2671618645381.6 16 f(x) = 1.0000000000
x= 16 Sinh^2(x) = 19740740045670. 668 f(x) = 0.9921875000
x= 17 Sinh^2(x) = 145865435631864 .219 f(x) = 1.0000000000
x= 18 Sinh^2(x) = 107780788677879 9.250 f(x) = 1.0000000000
x= 19 Sinh^2(x) = 796398293927843 8.000 f(x) = 0.0000000000[color=blue][color=green][color=darkred]
>>>[/color][/color][/color]

[Paul Rubin]

Re: The old round off problem?

"sam" <samschul@pacbe ll.net> writes:[color=blue]
> Hello all, I am taking a class in scientific programming at the local
> college. My problem is that the following difference produces round off
> errors as the value of x increases. For x >= 19 the diference goes to
> zero.I understand the problem, but am curious as to whether their
> exists a solution. [f(x) = cosh^2(x) - sinh^2(x) = 1][/color]

The usual way is to notice that the difference goes to zero because
the lowest order terms of the Taylor series for cosh^2 and sinh^2 are
equal and cancel each other out. The solution is to write down the
series for (cosh^2(x) - sinh^2(x)) and add up a few non-cancelled terms.
All these series should converge very fast.

[David Treadwell]

Re: The old round off problem?

On Mar 4, 2006, at 4:33 PM, Paul Rubin wrote:
[color=blue]
> "sam" <samschul@pacbe ll.net> writes:[color=green]
>> Hello all, I am taking a class in scientific programming at the local
>> college. My problem is that the following difference produces
>> round off
>> errors as the value of x increases. For x >= 19 the diference goes to
>> zero.I understand the problem, but am curious as to whether their
>> exists a solution. [f(x) = cosh^2(x) - sinh^2(x) = 1][/color]
>
> The usual way is to notice that the difference goes to zero because
> the lowest order terms of the Taylor series for cosh^2 and sinh^2 are
> equal and cancel each other out. The solution is to write down the
> series for (cosh^2(x) - sinh^2(x)) and add up a few non-cancelled
> terms.
> All these series should converge very fast.[/color]

Here's my analysis:

First, keep in mind that the range of values for a Taylor expansion
of sinh() and cosh() is limited. Then...

Write down the series for (cosh^2(x) - sinh^2(x)), as Paul suggested.

c2 = taylor(cosh(x)^ 2,x,8)

s2 = taylor(sinh(x)^ 2,x,8)

The first method gives the expected answer:

eval(c2-s2) ==> 1

Write down the series for cosh(x), square that; write the series for
sinh(x), square that and subtract.

Compare the two (I just did this in Eigenmath). Higher-order terms
remain for the second method.. I don't know how many terms my cosh()
and sinh() functions use, but if I evaluate (again, in Eigenmath)

The second method, which is more akin to what Python is doing (unless
there is an explicit 'sinh**2(x)' function) is different.

c2 = (taylor(cosh(x) ,x,8))^2, which gives terms up to and including
(x**8)**2 = x**16

s2 = (taylor(sinh(x) ,x,9))^2, which gives terms up to and including
x**18, then subtract the two, the result is

eval(c2 - s2) ==> 1 - 1/1814400 * x**10 - ... higher terms.

This begs the question, "What is the algorithm used by my Python?"

And, the ever-important question, "Does it matter in my final
result?" This is important, because I do not, in general, trust the
result when subtracting two large numbers.

:--David T.

[sam]

Re: The old round off problem?

David I beg I beg
Can you answer the question?

Also thanks for the information on using the Taylor series.

Sam Schulenburg

[David Treadwell]

Re: The old round off problem?

I wish I knew!

So I asked Google. Here's what I learned:

Most implementations are based on, or similar to the implementation
in the fdlibm package.

sinh(x) and cosh(x) are both based on exp(x). See http://

https://www.netlib.org/cgi-bin/netlibfiles.pl?filename=/fdlibm/e_sinh.c

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

On Mar 4, 2006, at 8:27 PM, sam wrote:
[color=blue]
> David I beg I beg
> Can you answer the question?
>
> Also thanks for the information on using the Taylor series.
>
> Sam Schulenburg
>
> --
> http://mail.python.org/mailman/listinfo/python-list[/color]

[sam]

Re: The old round off problem?


David Treadwell wrote:[color=blue]
> 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[/color]
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.