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

Hi,
I'm integrating mpc_pow with the fortran frontend and I'm getting one
regression related to raising a value to a huge integer.
FAIL: gfortran.dg/integer_exponentiation_4.f90 -O (test for errors, line
41)
FAIL: gfortran.dg/integer_exponentiation_4.f90 -O (test for errors, line
42)
FAIL: gfortran.dg/integer_exponentiation_4.f90 -O (test for excess errors)
The testcase has on lines 41/42 these statements:
print *, (2.0,-4.3)**huge(0) ! { dg-error "Arithmetic NaN" }
print *, (2.0,-4.3)**(-huge(0)) ! { dg-error "Arithmetic NaN" }
I.e. it expects a NaN for both cases. What I'm getting with MPC is
"Arithmetic overflow" on the first line because it overflows to Inf, and no
error for the second presumably because it underflows to zero.
I think perhaps expecting a NaN here is bogus and is an artifact of the
hand-written algorithm chosen to evaluate complex**int inside the fortran
frontend. The code is in fortran/arith.c:complex_pow(). I suspect the
algorithm is getting confused somewhere along the way.
E.g. if I write in C:
cpow (2 - 4.3i, INT_MAX);
I get (Inf - Inf I).
I propose that we change the lines in the testcase to instead say:
print *, (2.0,-4.3)**huge(0) ! { dg-error "Arithmetic overflow" }
print *, (2.0,-4.3)**(-huge(0))
Note this matches the non-complex behavior of pow (2**huge) and
pow(2**-huge) in the testcase. Also if I write e.g. (2.0,-4.3) ** (8.e9/3)
I also get arithmetic overflow for the non-integer calculation of a huge
number.
I don't know if it's worth fixing arith.c:complex_pow() since it'll be
removed when we convert to MPC. However the flip side would be that fortran
"failures" would appear until the switchover if I correct the testcase.
Also older versions of gfortran would still have the bug. I can open a PR
for this (like the complex division by zero issue.)
Thoughts?
--Kaveh
--
Kaveh R. Ghazi

Hi Kaveh,
Kaveh R. Ghazi wrote:
> The testcase has on lines 41/42 these statements:
>
> print *, (2.0,-4.3)**huge(0) ! { dg-error "Arithmetic NaN" }
> print *, (2.0,-4.3)**(-huge(0)) ! { dg-error "Arithmetic NaN" }
>
> I.e. it expects a NaN for both cases. What I'm getting with MPC is
> "Arithmetic overflow" on the first line because it overflows to Inf, and no
> error for the second presumably because it underflows to zero.
>
> I think perhaps expecting a NaN here is bogus and is an artifact of the
> hand-written algorithm chosen to evaluate complex**int inside the fortran
> frontend. The code is in fortran/arith.c:complex_pow(). I suspect the
> algorithm is getting confused somewhere along the way.
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.
FUNCTION f(z, e)
complex :: f, z, e
f = z**e
END
PRINT *, f((2.0, -4.3), huge(0))
END
should behave the same. With gfortran 4.1.2 this prints NaN.
Cheers,
- Tobi

From: "Tobias Schlüter"
> 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).
However when I print the two arguments inside the Fotran function I get:
( 2.0000000 , -4.3000002 )
( NaN, 2.0000000 )
NaN
So for some reason the exponent is turning into (NaN, 2) prior to
calculating the power value!
My Fortran is weak, here's what I had. Can you explain it?
Also why does it print one NaN instead of a complex NaN NaN?
Finally, if I use huge(0) - 1000000000, I get:
( 2.0000000 , -4.3000002 )
( 916.84369 , 2.0000000 )
+Infinity
Again, I can't explain the exponent values nor why it prints a real response
instead of a complex. (This is all with mainline on x86_64-linux-gnu).
FUNCTION f(z, e)
complex :: f, z, e
print *, z;
print *, e;
f = z**e
END
PRINT *, f((2.0, -4.3), huge(0))
END

Kaveh R. Ghazi wrote:
> From: "Tobias Schl=FCter"
>=20
>> 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.
>=20
> 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=20
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 =3D 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$
Cheers,
- Tobi
--=20
Tobias Schl=FCter
Am Coulombwall 1, Zi. 326
85748 Garching b. M=FCnchen
Tel.: +49/89/289-14139

From: "Tobias Schlüter"
Kaveh R. Ghazi wrote:
> From: "Tobias Schlüter"
>
>> 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.
Hi Tobias,
Let me amend my comment to say, the frontend should agree with the
library and/or optimizers >> assuming the library and/or optimizers are
correct.

From: "Steve Kargl"
> 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).
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.
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.)
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. :-)
Other ideas?
--Kaveh

On Thu, Jun 25, 2009 at 11:28:57AM -0700, Kaveh R. Ghazi wrote:
> From: "Steve Kargl"
>
> >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