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 55F7E384F6F0; Sat, 19 Nov 2022 19:11:42 +0000 (GMT)
DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org 55F7E384F6F0
DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org;
	s=default; t=1668885102;
	bh=ruPjts3Amy5LAfXro83E8Ii5ecZCX8g+htkKM6lir98=;
	h=From:To:Subject:Date:In-Reply-To:References:From;
	b=wMSQ6cPUlisJeYX09G/1JeVIclLkmRExwsLTGY47jfXzvj++k0iC/CS6oqZ8lUfgE
	 CQCQ3H0ff+iBGl2kjgkJRdtOQyTWxY+xYjxCRi2ipbsrsSDG17IjaJYCs0/0JnOoiv
	 EHqWmCcpnkAKABhJiISk1T25N3t9Gi1qkzqUNsuU=
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: Sat, 19 Nov 2022 19:11:41 +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: priority
Message-ID: <bug-107753-4-aiTb8H4N0e@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

kargl at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           Priority|P3                          |P4

--- Comment #10 from kargl at gcc dot gnu.org ---
(In reply to Steve Kargl from comment #9)
> On Fri, Nov 18, 2022 at 11:24:29PM +0000, sgk at troutmask dot
> apl.washington.edu wrote:
> >=20
> > Does anyone know what is meant by "Fortran rules"?  F66 does not
> > have any particular algorithm specified.  I'll look at F77 shortly.
> >=20
>=20
> Well, I hunted down the origins of -fcx-fortran-rules.
>=20
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=3D29549
>=20
> So, it appears to be an optimization, where Smith's algorithm
> will fail for extreme values of the real and imaginary parts
> of the complex number.

So, I wrong a dirty little program to time complex division.


program foo

   use timerm, only : rdtsc

   implicit none

   integer, parameter :: n =3D 1024*1024, dp =3D kind(1.d0)
   real(dp) re(n), im(n)
   complex(dp) x(n)
   integer i
   integer(8) t0, t

   t =3D 0
   do i =3D 1, 10
      call random_number(re)
      call random_number(im)
      x =3D cmplx(4 * re, 10 * im,8)
      t0 =3D rdtsc()
      x =3D x / x
      t =3D t + (rdtsc() - t0)
   end do
   print '(G0,1X,G0)', x(1)
   print *, real(t,10) / 10 / n
end program foo


Compiled with gfortran with its current method of doing division
(i.e., -fcx-fortran-rules), I see roughly 44.5 clock ticks per
division.  If run with a patched gfortran that uses the method that
the C compiler uses, I get about 62 ticks per division.  So, using
the stricter method impacts performance.

I'll note that gfortran unilaterally enforces -fcx-fortran-rules,
i.e., -fno-cx-fortran-rules has no effect.  Perhaps, gfortran=20
could be given a new -fcx-division=3DXXX option, where XXX is one of

   naive --> what -ffast-math does   (flags_complex_method =3D 0)
   smith --> what -fcx-fortran-rules (flags_complex_method =3D 1)
   strict -> default C behavior      (flags_complex_method =3D 2)=