[CLN-list] Floating point overflow discrepancy

[Richard B. Kreckel]

Hi!

On 02/24/2011 02:56 AM, Michael Miller wrote:
> I have looked at the code, and it seems that the problem is in
> scale_float, which throws a floating point exception if one tries to
> scale the exponent to more than 2^32-1.  This appears to be unnecessary.
>   For example, in the code below:
>
>     cl_I m="1073741824"; // 2^30
>     cl_LF x="1L0";
>     cout<<  3*m<<  " "<<  6*m<<  "\n";
>     cout<<  scale_float(scale_float(x, 3*m), 3*m)<<  "\n";
>     cout<<  scale_float(x, 6*m)<<  "\n";
>
> the results of the last two output lines should be identical, but I get
>
>     3221225472 6442450944
>     5.4667795091064729912L1939370979
>     terminate called after throwing an instance of
>     'cln::floating_point_overflow_exception'
>       what():  floating point overflow.
>     Aborted
>
> Now it's true that the exponent of the number being calculated (6*2^30)
> is larger than LF_exp_high (given as 2^32-1 in cl_LF.h), but that didn't
> interfere with the value computed by shifting twice by 3*2^30.  Is the
> comparison with LF_exp_high actually necessary?

Indeed, the routine scale_float in file cl_LF_scale_I.cc has not been 
adapted to 64-bit exponents on systems where uintD is 32 bit. This 
explains why it's broken on x86_32 only.

That comparison isn't the real problem. It's the assignments udelta = 
... The rhs is a 32 bit type, so truncation occurs! And of course the 
definition of IF_LENGTH is wrong, too. I suppose we should be able to 
fix all this by using the function cl_I_to_E for converting delta to an 
uintE, right?

Would you like to give it a try?

   -richy.
-- 
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>