Problems with numerical overflows/underflows in some functions

[Andreas113]

Hi Erik,
I have run numerous numerical tests over the past few days. I did not encounter any new parsing errors but I observed some incorrect results with overflow and underflow of input values in some mathematical functions like follows:

Function: hyperbolic cosine Cosh(..)
Example:

Cosh(Y*X):
----------
X            = -4987654.6211234567891011
Y            = -0.100000000000000009999995

MPA_Float    = 6.122507439387169799822007412143668689493139250589157593313679373744026E+216610 --> Exact!
QuadDouble   = NAN  --> WRONG! --> Correct for DoubleDouble range: +Inf
DoubleDouble = NAN  --> WRONG! --> Correct for QuadDouble   range: +Inf

Function: hyperbolic sine Sinh(..)
Example:

Sinh(Y*X):
----------
X            = -4987654.6211234567891011
Y            = -0.100000000000000009999995

MPA_Float    = 6.122507439387169799822007412143668689493139250589157593313679373744026E+216610 --> Exact!
QuadDouble   = NAN  --> WRONG! --> Correct for DoubleDouble range: +Inf
DoubleDouble = NAN  --> WRONG! --> Correct for QuadDouble   range: +Inf

Sinh(Y*X):
----------
X            = -4987654.6211234567891011
Y            = +0.100000000000000009999995

MPA_Float    = -6.122507439387169799822007412143668689493139250589157593313679373744026E+216610
QuadDouble   = NAN  --> WRONG! --> Correct for DoubleDouble range: -Inf
DoubleDouble = NAN  --> WRONG! --> Correct for QuadDouble   range: -Inf

I cannot judge which result the original C functions return: It would therefore be possible that there is an error there in the C-code, but also in the implementation of the +/- infinity values in Delphi is possible.

I will be testing further and reporting to you.

Kind regards,
Andreas

[neslib]

This happens because both DoubleDouble and QuadDouble have the same range as Double (10E+308) and the results of these calculations require a much larger range (10E+216610).

As you mentioned, the library should return +Inf or -Inf in this case (as would happen if you would use regular Double values).
However, the original C implementation returns NAN in this case, and there is not much I can do about that. But whether the result is NAN or INF, it means that the calculation cannot be performed.

So in short, there is not much I can do about this (other than modifying the C code, which I don't want to do since I am not that familiar with the underlying algorithms and math involved).

[Andreas113]

Hi Erik,
this is very disappointing and it is unfortunately not just a blemish, but a gross error in the C library. Let's assume we want to calculate the following:

Example:

1/Exp(Y*X):
-----------
X            = -800.0
Y            = -0.9
MPA_Float    = 2.032230802424293152866633766414812296724228166012459847311333592365251E-313
QuadDouble   = NAN  --> WRONG! --> Correct for QuadDouble range  : 0.0
DoubleDouble = NAN  --> WRONG! --> Correct for DoubleDouble range: 0.0

And this applies to all functions with asymptotic behaviour for +Inf or -Inf. A reliable mathematical library must be able to recognize such cases surely and to calculate them correctly!
Maybe I should write to professor David H. Bailey... Let's see.

Thank you & kind regards,
Andreas

[Andreas113]

Hi Erik,
Sorry, only today I had time to write to professor Bailey. He kindly answered me right away:

One difficulty here is that unfortunately, the person who wrote this software (approx 20 years ago) is no longer available to do any maintenance. Professor Bailey distributes the code, but is not sufficiently familiar with it to make any major changes on his own. But he will discuss this issue with one of his colleagues, check and respond.
So, we still have some hope for an update.

Regards, Andreas

[neslib]

Hi Andreas,

Thanks for going after this. If there is an update, I will be happy to apply the changes to this library as well!

Regards,
Erik

[neslib]

Hi Andreas,

I am closing this issue for now.
If there is a library update, I will be happy to implement it.