From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: by sourceware.org (Postfix, from userid 48) id 1B56C385223D; Fri, 18 Nov 2022 23:24:30 +0000 (GMT) DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org 1B56C385223D DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org; s=default; t=1668813870; bh=lplSvsZg8X6CSUD8M28d79lHt85RJo0XuY5wPLrM0r8=; h=From:To:Subject:Date:In-Reply-To:References:From; b=jJzNwfK4vz9+iwdqwcsndgPrY0aBXR+eoJJsixzXGvRWVr29waipWac+N18GzFpUX 5gHJxVXq5vx7kUV0Xeggms4upPk/LHiEMlactUTlZ27J0rJlUauoF8GQ9tg8FwS/Ax VIw/yrawTkCIqwDiyuuw5sRtE7YSxsB2lDTN3DSw= From: "sgk at troutmask dot apl.washington.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:24:29 +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: sgk at troutmask dot apl.washington.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 #5 from Steve Kargl -= -- On Fri, Nov 18, 2022 at 10:05:21PM +0000, kargl at gcc dot gnu.org wrote: > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=3D107753 >=20 > --- 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 Smit= h'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 perform= ance. > >=20 > > In any case, if the reporter wants to get robust results and in a porta= ble > > way, I would advise him to change/fix his algorithm accordingly. It ap= pears > > that a few other compilers behave here like gfortran. >=20 > It's likely coming from the middle-end where gcc.info has > the option >=20 > '-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. Does anyone know what is meant by "Fortran rules"? F66 does not have any particular algorithm specified. I'll look at F77 shortly. Tracking down what -fcx-fortran-rules does, one finds the eventually flag_complex_method is set to 1. The lower of complex division occurs in gcc/tree-complex.cc (expand_complex_division). If I use this patch % git diff gcc/tree-complex.cc | cat diff --git a/gcc/tree-complex.cc b/gcc/tree-complex.cc index ea9df6114a1..8051b7a3843 100644 --- a/gcc/tree-complex.cc +++ b/gcc/tree-complex.cc @@ -1501,6 +1501,7 @@ expand_complex_division (gimple_stmt_iterator *gsi, t= ree type, break; case 2: + case 1: if (SCALAR_FLOAT_TYPE_P (inner_type)) { expand_complex_libcall (gsi, type, ar, ai, br, bi, code, true= ); @@ -1508,7 +1509,6 @@ expand_complex_division (gimple_stmt_iterator *gsi, t= ree type, } /* FALLTHRU */ - case 1: /* wide ranges of inputs must work for complex divide. */ expand_complex_div_wide (gsi, inner_type, ar, ai, br, bi, code); break; to force gfortran through the C language code path, I get void doit (complex(kind=3D8) & restrict z) { 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; : _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>; _3 =3D __divdc3 (_7, _8, _9, _10); _11 =3D REALPART_EXPR <_3>; _12 =3D IMAGPART_EXPR <_3>; REALPART_EXPR <*z_5(D)> =3D _11; IMAGPART_EXPR <*z_5(D)> =3D _12; return; } with the result % gfcx -o z -fdump-tree-all a.f90 && ./z (1.79769313486231571E+308,1.79769313486231571E+308) (1.0000000000000000,0.0000000000000000) So, is -fcx-fortran-rules a relic of g77 past?=