From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 7330 invoked by alias); 1 Jun 2009 10:19:38 -0000 Received: (qmail 7218 invoked by uid 48); 1 Jun 2009 10:19:17 -0000 Date: Mon, 01 Jun 2009 10:19:00 -0000 Message-ID: <20090601101917.7217.qmail@sourceware.org> X-Bugzilla-Reason: CC References: Subject: [Bug fortran/40318] Complex division by zero in gfortran returns wrong results In-Reply-To: Reply-To: gcc-bugzilla@gcc.gnu.org To: gcc-bugs@gcc.gnu.org From: "dominiq at lps dot ens dot fr" Mailing-List: contact gcc-bugs-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Archive: List-Post: List-Help: Sender: gcc-bugs-owner@gcc.gnu.org X-SW-Source: 2009-06/txt/msg00009.txt.bz2 ------- 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