Re: Fortran pow (complex ** huge) returns NaN instead of Inf ?

[Steve Kargl <sgk at troutmask dot apl dot washington dot edu>]

On Thu, Jun 25, 2009 at 11:28:57AM -0700, Kaveh R. Ghazi wrote:
> From: "Steve Kargl" <sgk@troutmask.apl.washington.edu>
> 
> >troutmask:sgk[225] gfc4x -o z -fno-range-check c.f90
> >troutmask:sgk[226] ./z
> >(      -Infinity,      +Infinity)
> >(  2.0000000    , -4.3000002    )
> > 3.40282347E+38
> >(            NaN,            NaN)
> 
> Ah I didn't know about -fno-range-check.  FWIW, MPC passed these inputs 
> also yields (-Inf, Inf) at precision 24 for float.  So it won't be changing 
> gfortran's behavior for this case of pow (foo ** real4).

Without the ieee754 intrinsic module, the Fortran standard is
vague/silent on exceptional values.  gfortran's constant folding
takes the tack that it will issue errors on exceptional values,
and thus it does range checking.  -fno-range-check option was
added so people could generate Inf and NaN for testing, ie.,

real, paramete :: inf = 1. / 0., nan = 0. / 0.
...
 if (x == nan .or. x == inf) then
    do some damage control
 end if

> The difference that remains is my original example which is in the pow (foo 
> ** int) handler in the fortran.  I believe the frontend code to compute 
> this has a bug and yields NaN incorrectly.  Using MPC there fixes it, but 
> breaks the testcase because it incorporates the buggy value in its test.

I now understanding why gfortran gets NaN in this case.  It's easiest
to understand by looking at libgfortran/generated/pow_c4_i4.c.  In
particular, the comment is

/* Use Binary Method to calculate the powi. This is not an optimal but
   a simple and reasonable arithmetic. See section 4.6.3, "Evaluation of
   Powers" of Donald E. Knuth, "Seminumerical Algorithms", Vol. 2, "The Art
   of Computer Programming", 3rd Edition, 1998.  */

The main loops is

      for (;;)
        {
          if (u & 1)
            pow *= x;
          u >>= 1;
          if (u)
            x *= x;
          else
            break;
        }

For x = a + i b, x *= x becomes (a*a - b*b) + i 2 a b.  At some
point, the a*a and b*b are Inf, and the difference yields NaN.

The same algorithm is used in the constant foldings!

> So the options I see for using mpc_pow are:
> 
> 1.  Fix the testcase and live with the testcase failing for people who 
> don't use MPC > 0.6 (which has mpc_pow).
> 
> 2.  Leave the testcase as-is until the MPC hard-requirement switchover.  It 
> will fail for anyone who uses MPC > 0.6.
> 
> 3.  Comment out the offending lines in the testcase, put them back in after 
> the hard-requirement switch.  Nobody will see extra failures in this case. 
> (I'll add that to PR40302's todo list.)

I think this is the sanest way to move forward.  Options 1 and 2
will lead to numerous bug reports.

> 4.  Fix the testcase, and fix arith.c:complex_pow() so that it doesn't get 
> a NaN.  It's probably the "right" solution, but I don't volunteer to fix 
> the old code. :-)

I'll think about improvements to the algorithm.

-- 
Steve