Expected precision of Fortran FE constant evaluation

[Jakub Jelinek]

On Sat, Feb 09, 2019 at 01:00:20PM -0800, Steve Kargl wrote:
> > The Fortran standard doesn't specify in which order an expression is
> > evaluated, so intermediates in an expression might or might not overflow
> > depending on the order of evaluation (say, "a + b - c", if all numbers are
> > slightly above HUGE()/2). AFAIK the rules for compile-time evaluation are
> > no different than for runtime. So from that perspective it seems better to
> > me to range check only the results, and not every intermediate value,
> > although AFAICT both are standard conforming.
> > 
> 
> You're correct in that Fortran does not mandate anything about 
> order of evaluation.  gfortran does left to right, and range
> check the intermediate values in the evaluation of the RHS.
> 
> 
>   print *, (huge(1.) / 1.9) +  (huge(1.) / 1.9) - (huge(1.) / 1.9)
>   print *, (huge(1.) / 1.9) + ((huge(1.) / 1.9) - (huge(1.) / 1.9))
> end 
> 
> The first line causes an error.  The second line produces 1.79095E38.

For norm2, it can indeed do that in one pass (both for compile time and
runtime), but I think it can't check just the result but also needs to check
each operand before using **2 on it.  By using the power of two scales
(actually power of four) it can easily scale both the result which contains
sum of **2 and the operands before **2.  The result can be perhaps checked
only in certain iterations, e.g. by tracking the largest exponent seen so
far and assuming that every addition could have increased that exponent by
one.  So:
for (...) {
  temp = op[cur];
  temp *= scale1;
  if (temp >= constant (sqrt (huge()) / 64.0))
    {
      old_scale2 = scale2;
      frexp; adjust scale1 and scale2;
      result *= scale2 / old_scale2;
    }
  if (result is too large)
    {
      adjust scale1 and scale2;
      result *= scale2 / old_scale2;
    }
  result += temp * temp;
}
Though, perhaps the scale for the operands can be adjusted just once and by
fixed amount then, i.e. if the number is above say that 2**(__???_MAX_EXP__
/ 2 - 10) then make sure scale is 0.5**(__???_MAX_EXP__ / 2 + 6) or smaller,
then it will not overflow with any finite number.  The scale for the result
might need some tracking on where it would need adjustment.

	Jakub