From mboxrd@z Thu Jan  1 00:00:00 1970
Return-Path: <gcc-bugzilla@gcc.gnu.org>
Received: by sourceware.org (Postfix, from userid 48)
	id A54E53864A0F; Wed,  7 Dec 2022 21:50:14 +0000 (GMT)
DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org A54E53864A0F
DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org;
	s=default; t=1670449814;
	bh=8ntNTz6J20klDzZru69RMTnrpif3HJJ9BYFGswbJ2A8=;
	h=From:To:Subject:Date:In-Reply-To:References:From;
	b=B4JHBYuSqrG+93VH2+glKZaM/vFQ9aUr0LgWIovgz6UBVDuIaVcMOVz9R2VwRmpDh
	 kmCQfY28HTk0MJU1aA1VVDSw5VlSu1HJBo4882Nk29K1R4u0UskqqhToI7/51PN6wS
	 BU4yxZIOPrIF9uDcK6mjp0ozYLycT5vo9Pah3KyY=
From: "kargl at gcc dot gnu.org" <gcc-bugzilla@gcc.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: Wed, 07 Dec 2022 21:50:14 +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: P4
X-Bugzilla-Assigned-To: unassigned at gcc dot gnu.org
X-Bugzilla-Target-Milestone: ---
X-Bugzilla-Flags: 
X-Bugzilla-Changed-Fields: 
Message-ID: <bug-107753-4-4QBykbplCW@http.gcc.gnu.org/bugzilla/>
In-Reply-To: <bug-107753-4@http.gcc.gnu.org/bugzilla/>
References: <bug-107753-4@http.gcc.gnu.org/bugzilla/>
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: <gcc-bugs.sourceware.org>

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=3D107753

--- Comment #14 from kargl at gcc dot gnu.org ---
(In reply to anlauf from comment #13)
> (In reply to Steve Kargl from comment #12)
> > The optimization level is irrelevant.  gfortran unilaterally
> > uses -fcx-fortran-rules, and there is no way to disable this
> > option to user the slower, but stricter, evaluation.  One
> > will always get complex division computed by
> >=20
> > a+ib   a + b(d/c)     b - a(d/c)=20
> > ---- =3D ---------- + i ------------  |c| > |d|
> > c+id   c + d(d/c)     c + d(d/c)
> >=20
> > and similar for |d| > |c|.
> >=20
> > There are a few problems with this. d/c can trigger an invalid underflow
> > exception.  If d =3D=3D c, you then have numerators of a + b and b - a,=
 you
> > can get a invalid overflow for a =3D huge() and b > 1e291_8.
>=20
> I am wondering how slow an algorithm would be that scales numerator
> and denominator by respective factors that are powers of 2, e.g.
>=20
> e_num =3D 2. ** -max (exponent (a), exponent (b))
> e_den =3D 2. ** -max (exponent (c), exponent (d))
>=20
> The modulus of scaled values would be <=3D 1, even for any of a,... being
> huge().
> Of course this does not address underflows that could occur during scalin=
g,
> or denormalized numbers, which are numerically irrelevant for the result.
>=20
> Is there anything else wrong with this approach?

Comment #10 contains a simple timing measurement in from my Intel Core2 Duo
based system.  gfortran with its current method (ie., -fcx-fortran-rules) t=
akes
44.5 clock ticks for a complex division.  If I sidestep the option and forc=
e it
to use the C language method of evaluation, it takes 62 clock ticks.  I hav=
en't
looked at what algorithm C uses, but I suspect its along the lines you sugg=
est.
 The question is likely do we break backwards compatibility and remove
-fcx-fortran-rules or change when/how -fcx-fortran-rules applies (e.g., add=
 it
to -ffast-math?)=