From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: by sourceware.org (Postfix, from userid 48) id 4E68738518B1; Fri, 18 Nov 2022 23:45:05 +0000 (GMT) DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org 4E68738518B1 DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org; s=default; t=1668815105; bh=tPtjzsGeRvv1apecDRqn2CbBp3a/V3ioSiUcyjuKo6I=; h=From:To:Subject:Date:In-Reply-To:References:From; b=bJIJEN8wZ6Zw7Snc96BZiCARmdlZ0DEUKmHJV+RYIkTzcyc8OQRvvTv4RIqA9nYWq wSfg4fNWgWe1KJ46YnM1lRi67Raiw1lJeGgRc1IE0P1YIseQPN1WKsC3tIuQwAb2x+ aQNb5OFTmdTgnrDtGKJeCLkMmu2KxszRd+KL8GJ4= From: "weslley.pereira at ucdenver dot edu" To: gcc-bugs@gcc.gnu.org Subject: [Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I) Date: Fri, 18 Nov 2022 23:45:04 +0000 X-Bugzilla-Reason: CC X-Bugzilla-Type: changed X-Bugzilla-Watch-Reason: None X-Bugzilla-Product: gcc X-Bugzilla-Component: fortran X-Bugzilla-Version: 12.2.0 X-Bugzilla-Keywords: X-Bugzilla-Severity: normal X-Bugzilla-Who: weslley.pereira at ucdenver dot edu X-Bugzilla-Status: WAITING X-Bugzilla-Resolution: X-Bugzilla-Priority: P3 X-Bugzilla-Assigned-To: unassigned at gcc dot gnu.org X-Bugzilla-Target-Milestone: --- X-Bugzilla-Flags: X-Bugzilla-Changed-Fields: Message-ID: In-Reply-To: References: Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable X-Bugzilla-URL: http://gcc.gnu.org/bugzilla/ Auto-Submitted: auto-generated MIME-Version: 1.0 List-Id: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=3D107753 --- Comment #7 from Weslley da Silva Pereira --- (In reply to anlauf from comment #3) > I guess the reporter assumes that gcc uses a clever algorithm like Smith's > to handle such extreme cases of complex division. Not sure if that one is > available by some compilation flag, and I think it would impact performan= ce. >=20 > In any case, if the reporter wants to get robust results and in a portable > way, I would advise him to change/fix his algorithm accordingly. It appe= ars > that a few other compilers behave here like gfortran. Thanks for the suggestion of changing the algorithm that needs such a divis= ion. What are the ranges for nominator and denominator where one should rely on = the intrinsic complex division? Maybe this is the good question to ask. Then we= can build algorithms that attend such requisites. LAPACK has cladiv and zladiv, which are routines for complex division that avoids unnecessary under and overflow. They are used in many parts of the code. This implementation is f= or sure less efficient than the intrinsic complex division, but we rely on it because of robustness. More data for the discussion: 1. In a Ubuntu 18.04.5 LTS, using GNU Fortran 7.5.0, I tested optimization flags `-O` but still reproduce the wrong result for complex divisions with = huge numbers. See https://github.com/Reference-LAPACK/lapack/issues/575#issuecomment-910616816 that used the code from https://github.com/Reference-LAPACK/lapack/blob/master/INSTALL/test_zcomple= xdiv.f. This is the test currently in LAPACK 3.11.0. 2. I have just reproduced what was reported in https://github.com/Reference-LAPACK/lapack/issues/575#issuecomment-91061681= 6 in my Ubuntu 20.04.5 LTS, using GNU Fortran 9.4.0. 3. I noticed that the optimization flag is unable to target divisions like `x/x` depending on where they are inside a program. 4. My Ubuntu 20.04.5 LTS with compiler ifort 2021.7.1 computes the complex division `x/x` accurately even for the case of huge numbers. Scenarios test= ed: - I tested the program in https://github.com/Reference-LAPACK/lapack/blob/master/INSTALL/test_zcomple= xdiv.f and the one in https://godbolt.org/z/b3WKWodvn. - I tested ifort with flags -fp-model precise and -fp-model fast. The la= tter enables more aggressive optimizations on floating-point data. - I tested compilation with optimization flags -O0, -O, -O1, -O2, -O3.=20 Here is the implementation of the complex division in LAPACK if it somehow helps the discussion: https://netlib.org/lapack/explore-html/d8/d9b/group__double_o_t_h_e_rauxili= ary_gad1c0279ec29e8ac222f1e319f4144fcb.html=