From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: by sourceware.org (Postfix, from userid 48) id 11FAA3854563; Fri, 18 Nov 2022 22:05:22 +0000 (GMT) DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org 11FAA3854563 DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org; s=default; t=1668809122; bh=BQruRG57J9YrojCnjyTY69fmvrpygBVJzGyKI0h52eM=; h=From:To:Subject:Date:In-Reply-To:References:From; b=Tus9TZ6uTnAd9R/AqH3npf3oTGlou8M/1OZIEGsZJ3nQ6ez/4LlMjPWXeBGu5CXLY yaJaW0wmuQXF3zvs4qgOVWo6FHU4ubaFkj6fRsnmtzUYM4l5MPlqrjo8Nz0OHbvaZe jk2GDZKrUL7Hl3nyUd9Hrzkf4RXgeAMtgTQHz+kY= From: "kargl at gcc dot gnu.org" 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 22:05:21 +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: kargl at gcc dot gnu.org 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 #4 from kargl at gcc dot gnu.org --- (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. It's likely coming from the middle-end where gcc.info has the option '-fcx-fortran-rules' Complex multiplication and division follow Fortran rules. Range reduction is done as part of complex division, but there is no checking whether the result of a complex multiplication or division is 'NaN + I*NaN', with an attempt to rescue the situation in that case. Consider program zdiv implicit none integer, parameter :: dp =3D kind(1.d0) real(dp) r complex(dp) :: y ! r =3D huge(1.d0) / 2 ! yields (1,0) ! r =3D nearest(huge(1.d0), -1.d0) ! yields (1,0) r =3D huge(1.d0) ! yields (NaN,0) y =3D cmplx(r, r, dp) write(*,*) y call doit(y) write(*,*) y contains subroutine doit(z) complex(dp) z z =3D z / z end subroutine doit end program zdiv If you compile this with -fdump-tree-all, then one gets (I've added annotation with <--- marker) % more z-a.f90.241t.cplxlower0 __attribute__((fn spec (". w "))) void doit (complex(kind=3D8) & restrict z) { real(kind=3D8) D.4265; real(kind=3D8) D.4264; complex(kind=3D8) _1; complex(kind=3D8) _2; complex(kind=3D8) _3; real(kind=3D8) _7; real(kind=3D8) _8; real(kind=3D8) _9; real(kind=3D8) _10; real(kind=3D8) _11; real(kind=3D8) _12; logical(kind=3D1) _13; real(kind=3D8) _14; real(kind=3D8) _15; real(kind=3D8) _16; real(kind=3D8) _17; real(kind=3D8) _18; real(kind=3D8) _19; real(kind=3D8) _20; real(kind=3D8) _21; real(kind=3D8) _22; real(kind=3D8) _23; real(kind=3D8) _24; real(kind=3D8) _25; real(kind=3D8) _26; real(kind=3D8) _27; real(kind=3D8) _28; real(kind=3D8) _29; real(kind=3D8) _30; real(kind=3D8) _31; real(kind=3D8) _32; real(kind=3D8) _33; real(kind=3D8) _34; real(kind=3D8) _35; real(kind=3D8) _36; real(kind=3D8) _37; real(kind=3D8) _38; real(kind=3D8) _39; : _7 =3D REALPART_EXPR <*z_5(D)>; _8 =3D IMAGPART_EXPR <*z_5(D)>; _1 =3D COMPLEX_EXPR <_7, _8>; _9 =3D REALPART_EXPR <*z_5(D)>; _10 =3D IMAGPART_EXPR <*z_5(D)>; _2 =3D COMPLEX_EXPR <_9, _10>; _11 =3D ABS_EXPR <_9>; _12 =3D ABS_EXPR <_10>; _13 =3D _11 < _12; if (_13 !=3D 0) goto ; [50.00%] else goto ; [50.00%] : _14 =3D _9 / _10; <--- Should be one _15 =3D _9 * _14; <--- huge(1.d0) _16 =3D _15 + _10; <--- huge(1.d0) + huge(1.d0) =3D inf _17 =3D _7 * _14; <--- huge(1.d0) _18 =3D _17 + _8; <--- huge(1.d0) + huge(1.d0) =3D inf _19 =3D _8 * _14; <--- huge(1.d0) _20 =3D _19 - _7; <--- huge(1.d0) - huge(1.d0) =3D 0 _21 =3D _18 / _16; <--- inf / inf =3D NaN _22 =3D _20 / _16; <--- 0 / inf =3D 0 _36 =3D _21; _37 =3D _22; goto ; [100.00%] : _23 =3D _10 / _9; _24 =3D _10 * _23; _25 =3D _24 + _9; _26 =3D _8 * _23; _27 =3D _26 + _7; _28 =3D _7 * _23; _29 =3D _8 - _28; _30 =3D _27 / _25; _31 =3D _29 / _25; _38 =3D _30; _39 =3D _31; : # _34 =3D PHI <_36(4), _38(5)> # _35 =3D PHI <_37(4), _39(5)> _3 =3D COMPLEX_EXPR <_34, _35>; _32 =3D _34; _33 =3D _35; REALPART_EXPR <*z_5(D)> =3D _32; IMAGPART_EXPR <*z_5(D)> =3D _33; return; }=