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

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

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

--- Comment #11 from anlauf at gcc dot gnu.org ---
(In reply to Weslley da Silva Pereira from comment #7)
> More data for the discussion:
> 1. In a Ubuntu 18.04.5 LTS, using GNU Fortran 7.5.0, I tested optimization
> flags `-O` but still reproduce the wrong result for complex divisions with
> huge numbers. See

It is possible that gfortran's dependence on optimization level depends
on the version.  If one wants to test run-time behavior and avoid
compile-time simplification, it may be helpful to add:

  volatile :: x, y, z

I then get consistent results for -O0 / -O1.

> 4. My Ubuntu 20.04.5 LTS with compiler ifort 2021.7.1 computes the complex
> division `x/x` accurately even for the case of huge numbers. Scenarios
> tested:
>    - I tested the program in
> https://github.com/Reference-LAPACK/lapack/blob/master/INSTALL/
> test_zcomplexdiv.f and the one in https://godbolt.org/z/b3WKWodvn.
>    - I tested ifort with flags -fp-model precise and -fp-model fast. The
> latter enables more aggressive optimizations on floating-point data.
>    - I tested compilation with optimization flags -O0, -O, -O1, -O2, -O3. 

Intel might be fine, but at least some current llvm-based compilers (Nvidia,
AMD flang) show more or less similar behavior to gfortran.

E.g. nvfortran 22.11:

 (1.7976931348623157E+308,1.7976931348623157E+308)
 (8.9884656743115795E+307,8.9884656743115795E+307)
 (4.4942328371557898E+307,4.4942328371557898E+307)
 (NaN,0.000000000000000)
 (NaN,0.000000000000000)
 (1.000000000000000,0.000000000000000)

As a sidenote: we are really discussing borderline cases here, valid
but only rarely occuring in normal code execution.  If I replace

  x = cmplx( huge(0.0d0), huge(0.0d0), dp )
  y = cmplx( b**(E-1), b**(E-1), dp )

by

  x = cmplx( nearest(huge(0.0d0),-1.d0), nearest(huge(0.0d0),-1.d0), dp )
  y = cmplx( nearest(b**(E-1),   -1.d0), nearest(b**(E-1),   -1.d0), dp )

then I get

               (1.0000000000000000,0.0000000000000000)
               (1.0000000000000000,0.0000000000000000)

instead of NaN.