[CLN-list] Floating point overflow discrepancy

[Michael Miller]

On 02/22/2011 02:52 AM, Richard B. Kreckel wrote:
> On 02/21/2011 07:05 AM, Michael Miller wrote:
>> The following code (to calculate x^100000):
>>
>>     cl_R x="1.0L100000", y=1;
>>     for (int i=1; i<=100000; i++) y=x*y;
>>     cout<<  y<<  "\n";
>>     cout<<  exp(100000*ln(x))<<  "\n";;
>>
>> results in the output
>>
>>     1.00000000000000302104L10000000000
>>     terminate called after throwing an instance of
>>       'cln::floating_point_overflow_exception'
>>     what():  floating point overflow.
>>     Aborted
>>
>> I know that these are very large numbers, but the two calculations are
>> computing the same value.  If the first doesn't overflow, then why does
>> the second?
> 
> It is not only the final result that matters. When an intermediate
> result overflows an exception is thrown, too. You may wish to have a
> look at the code in src/float/transcendental/.

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?

MM