[Bug fortran/40318] Complex division by zero in gfortran returns wrong results

[dominiq at lps dot ens dot fr]

------- Comment #3 from dominiq at lps dot ens dot fr  2009-06-01 10:19 -------
In http://gcc.gnu.org/ml/fortran/2009-06/msg00006.html, Dennis Wassel 
wrote:
> Complex numbers have a well-defined concept of infinity, which I like to
> visualise as the "infinite-diameter" ring around the finite complex
> numbers, so its argument is, of course, not well-defined.  This is
> conceptually somewhat different from the two points +/-Inf on the real
> line, and it is not signed.

My understanding of infinity in the complex plane is what is called (I
call?)  "directed inifinity": if abs((a,b)) goes to +Inf and atan2(a,b) has
a defined value in this limit, then (a,b) goes to infinity in the direction
given by atan2(a,b).  However (+/-Inf,+/-Inf) defines only four directions
and is unable to represent general "directed inifinity".  So I think that
from a "mathematical" point of view the problem is ill-posed,
(+/-Inf,+/-Inf) is in the same class of "undefined" values as 0/0 or
Inf/Inf, and should give NaNs.

Now atan2(a,b) is DEFINED such that atan2(+0,+/-0)=+/-0,
atan2(-0,+/-0)=+/-Pi, atan2(+Inf,+/-Inf)=+/-Pi/4,
atan2(-Inf,+/-Inf)=+/-3*Pi/4 (it seems that it is even built in the Intel
hardware).  With this definition of atan2, it is possible to give a
definition of (+/-Inf,+/-Inf) as the directions of the corners of the
"infinite square".  I have nothing against this defintion, except it should
be documented.

Part of the problem originates from the optimization of (a,b)/(c,0.0) as
(a/c,b/c), see pr30482, and it limit when c goes to 0 (the latter giving
(+/-Inf,+/-Inf) if b/=0.0 or has a unknown value at compile time,
(+/-Inf,+/-0) if b known to be zero at compile time, optimized as
(a/c,+/-0.0), or (+/-Inf,-/+Nan) if b==0 at run time).

Note that due to the definition of atan2(+/-0,+/-0),
(+/-0,+/-0)=0*exp((+/-k*Pi,+/-0)) (k being 0 or 1), then to be consistent
1.0/(+/-0.0,+/-0.0) should give 1/0.0*exp((-/+k*Pi,-/+0))=(+/-Inf,-/+Nan).

Last remark, if I remember correctly (I cannot find the right pointer) C99 
defines precisely how complex multiplications should behave in the 
exception limits, with the drawback of a large increase of the computation 
cost avoided with the -fcx-fortran-rules option.

As far as I can tell:

(1) Without the IEEE module, using or producing Inf or NaN makes the code 
non (fortran) standard comforming, so the "processor" can give any answer.

(2) With the IEEE module, Inf and NaN are part of the numerical model, 
however I did not find any definition of the values that intrinsincs 
should return in such cases.

To conclude, I'll repeat that in my opinion the only relevant answer for
overflow, divide by zero, or invalid argument exceptions is "call abort()".
Since I know that it is a lost battle, I think at least the behavior shall
not increase the computation time with complex numbers (or provide a way to
avoid the penalty).  Now if the fxxxx standard does not define the expected
result, but the C99 does, the best answer to this ill-posed problem is
probably to follow the C99 standard if possible (I have no idea of what
does mpc).  In this case the results at compile and run times should be the
same ("least surprising approach").


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=40318