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?
Re: Fortran pow (complex ** huge) returns NaN instead of Inf ? by Tobias Schlüter on
2009-06-25T08:01:00+00:00
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
Re: Fortran pow (complex ** huge) returns NaN instead of Inf ? by Kaveh R. Ghazi on
2009-06-25T16:22:39+00:00
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 ( 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
Re: Fortran pow (complex ** huge) returns NaN instead of Inf ? by Tobias Schlüter on
2009-06-25T16:35:01+00:00
Kaveh R. Ghazi wrote:
> From: "Tobias Schl=FCter" <tobias.schlueter@physik.uni-muenchen.de>
>=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 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
Re: Fortran pow (complex ** huge) returns NaN instead of Inf ? by Kaveh R. Ghazi on
2009-06-25T18:10:01+00:00
From: "Tobias Schlüter" <Tobias.Schlueter@physik.uni-muenchen.de>
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.
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. <<
I.e. we should not cater to bugs in the C library or the middle-end
optimizers if they exist. That your program prints (Inf, NaN) means to me
there's a bug in your libc cpowf. Note my glibc gives three different
values depending on whether one uses cpowf, cpow or cpowl. So certainly
there is a problem in the C library here.
cpowf (2 -4.3I, 3.40282347E+38) -> (inf nan)
cpow (2 -4.3I, 3.40282347E+38) -> (inf inf)
cpowl (2 -4.3I, 3.40282347E+38) -> (-inf -inf)
We have a long history when using e.g. MPFR to fold tan() that we go with
the correct answer rather than try to duplicate bugs in a particular C
library implementation. I believe the correct answer here does not involve
NaN in either the real or imaginary parts.
Now running your program side-by-side with the same values folded by the
fortran frontend gives an Arithmetic overflow at compile-time. So I'm not
even sure this example is meaningful to say the frontend and library should
be the same.
Do you agree?
print *, e
f = z**e
print *, f
end
print *, f((2.0, -4.3), huge(0.))
print *, (2.0, -4.3) ** huge(0.)
END
Re: Fortran pow (complex ** huge) returns NaN instead of Inf ? by Kaveh R. Ghazi on
2009-06-25T18:29:09+00:00
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).
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?
Re: Fortran pow (complex ** huge) returns NaN instead of Inf ? by Steve Kargl on
2009-06-25T19:15:20+00:00
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.