[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

["kargl at gcc dot gnu.org" <gcc-bugzilla@gcc.gnu.org>]

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

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:
> > 
> > Does anyone know what is meant by "Fortran rules"?  F66 does not
> > have any particular algorithm specified.  I'll look at F77 shortly.
> > 
> 
> Well, I hunted down the origins of -fcx-fortran-rules.
> 
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=29549
> 
> 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 = 1024*1024, dp = kind(1.d0)
   real(dp) re(n), im(n)
   complex(dp) x(n)
   integer i
   integer(8) t0, t

   t = 0
   do i = 1, 10
      call random_number(re)
      call random_number(im)
      x = cmplx(4 * re, 10 * im,8)
      t0 = rdtsc()
      x = x / x
      t = 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 
could be given a new -fcx-division=XXX option, where XXX is one of

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