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 06:34:48PM +0200, Tobias Schlüter wrote:
> Kaveh R. Ghazi wrote:
> >From: "Tobias Schlüter" <tobias.schlueter@physik.uni-muenchen.de>
> >
> >>I agree with everything you say, but I'd like to point out that one thing
> >>you should make sure is, that compiled code (or code generated by the
> >>optimizers) and code evaluated by the frontend yield the same results,
> >>i.e.
> >>[...]
> >>should behave the same.  With gfortran 4.1.2 this prints NaN.
> >
> >I agree with that.  However when investigating *why* it returns NaN I get
> >confused.  I did -fdump-tree-original on your program and it appears in
> >function "F" it just calls __builtin_cpowf on the two arguments.  The
> >equivalent C program that calls f(2-4.3I, INT_MAX) where f() just calls
> >cpowf returns (inf inf).
> 
> Sorry, Kaveh, my Fortran was wrong, the exponent was not correctly 
> typed.  Here's a corrected program:
> Tobias.Schlueter@zuppc13 src$ cat t.f90
> FUNCTION f(z, e)
>  complex :: f, z
>  real e
>  print *, z
>  print *, e
>  f = z**e
>  print *, f
> end
> 
> print *, f((2.0, -4.3), huge(0.))
> END
> Tobias.Schlueter@zuppc13 src$ gfortran t.f90
> Tobias.Schlueter@zuppc13 src$ ./a.out
>  (  2.000000    , -4.300000    )
>   3.4028235E+38
>  (     +Infinity,           NaN)
>             NaN
> Tobias.Schlueter@zuppc13 src$
> 

Tobi,

Here's a slightly modified version where the interface matches.

troutmask:sgk[223] cat c.f90
program test
  print *, (2.0, -4.3)**huge(0.)
  print *, f((2.0, -4.3), huge(0.))
  contains
  function f(z, e)
    complex :: f, z
    real e
    print *, z
    print *, e
    f = z**e
  end function f
end program test
troutmask:sgk[224] gfc4x -o z c.f90
c.f90:2.22:

  print *, (2.0, -4.3)**huge(0.)
                      1
Error: Arithmetic overflow at (1)
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)

This is without MPC, so gfortran already gets inconsistent results.
gfortran constant folding generates 

    {
      static complex(kind=4) C.1545 = __complex__ ( -Inf,  Inf);

      _gfortran_transfer_complex (&dt_parm.2, &C.1545, 4);
    }

which is probably acceptable when one thinks of the polar representation
of (2,-4.3) = 4.74 exp(- i phi) where  -pi < phi <= pi (too lazy to compute
phi).

Furthermore, for the function, the dump shows

  {
    complex(kind=4) D.1542;

    D.1542 = *z;
    __result_f = __builtin_cpowf (D.1542, COMPLEX_EXPR <*e, 0.0>);
  }

I suspect __builtin_cpowf() is doing something wrong due to the
promotion of real *e to COMPLEX_EXP <*e, 0.>, but I haven't
tried to investigate any further.


-- 
Steve