[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

--- 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 Smith'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 performance.
> 
> In any case, if the reporter wants to get robust results and in a portable
> way, I would advise him to change/fix his algorithm accordingly.  It appears
> that a few other compilers behave here like gfortran.

It's likely coming from the middle-end where gcc.info has
the option

'-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.

Consider

program zdiv
   implicit none
   integer, parameter :: dp = kind(1.d0)
   real(dp) r
   complex(dp) :: y

!   r = huge(1.d0) / 2              ! yields (1,0)
!   r = nearest(huge(1.d0), -1.d0)  ! yields (1,0)
   r = huge(1.d0)                   ! yields (NaN,0)
   y = cmplx(r, r, dp)
   write(*,*) y
   call doit(y)
   write(*,*) y
   contains
      subroutine doit(z)
         complex(dp) z
         z = z / z
      end subroutine doit
end program zdiv

If you compile this with -fdump-tree-all, then one gets
(I've added annotation with <--- marker)


% more z-a.f90.241t.cplxlower0

__attribute__((fn spec (". w ")))
void doit (complex(kind=8) & restrict z)
{
  real(kind=8) D.4265;
  real(kind=8) D.4264;
  complex(kind=8) _1;
  complex(kind=8) _2;
  complex(kind=8) _3;
  real(kind=8) _7;
  real(kind=8) _8;
  real(kind=8) _9;
  real(kind=8) _10;
  real(kind=8) _11;
  real(kind=8) _12;
  logical(kind=1) _13;
  real(kind=8) _14;
  real(kind=8) _15;
  real(kind=8) _16;
  real(kind=8) _17;
  real(kind=8) _18;
  real(kind=8) _19;
  real(kind=8) _20;
  real(kind=8) _21;
  real(kind=8) _22;
  real(kind=8) _23;
  real(kind=8) _24;
  real(kind=8) _25;
  real(kind=8) _26;
  real(kind=8) _27;
  real(kind=8) _28;
  real(kind=8) _29;
  real(kind=8) _30;
  real(kind=8) _31;
  real(kind=8) _32;
  real(kind=8) _33;
  real(kind=8) _34;
  real(kind=8) _35;
  real(kind=8) _36;
  real(kind=8) _37;
  real(kind=8) _38;
  real(kind=8) _39;

  <bb 2> :
  _7 = REALPART_EXPR <*z_5(D)>;
  _8 = IMAGPART_EXPR <*z_5(D)>;
  _1 = COMPLEX_EXPR <_7, _8>;
  _9 = REALPART_EXPR <*z_5(D)>;
  _10 = IMAGPART_EXPR <*z_5(D)>;
  _2 = COMPLEX_EXPR <_9, _10>;
  _11 = ABS_EXPR <_9>;
  _12 = ABS_EXPR <_10>;
  _13 = _11 < _12;
  if (_13 != 0)
    goto <bb 4>; [50.00%]
  else
    goto <bb 5>; [50.00%]

  <bb 4> :
  _14 = _9 / _10;  <--- Should be one
  _15 = _9 * _14;  <--- huge(1.d0)
  _16 = _15 + _10; <--- huge(1.d0) + huge(1.d0) = inf
  _17 = _7 * _14;  <--- huge(1.d0)
  _18 = _17 + _8;  <--- huge(1.d0) + huge(1.d0) = inf
  _19 = _8 * _14;  <--- huge(1.d0)
  _20 = _19 - _7;  <--- huge(1.d0) - huge(1.d0) = 0
  _21 = _18 / _16; <--- inf / inf = NaN
  _22 = _20 / _16; <--- 0 / inf = 0
  _36 = _21;
  _37 = _22;
  goto <bb 3>; [100.00%]

  <bb 5> :
  _23 = _10 / _9;
  _24 = _10 * _23;
  _25 = _24 + _9;
  _26 = _8 * _23;
  _27 = _26 + _7;
  _28 = _7 * _23;
  _29 = _8 - _28;
  _30 = _27 / _25;
  _31 = _29 / _25;
  _38 = _30;
  _39 = _31;

  <bb 3> :
  # _34 = PHI <_36(4), _38(5)>
  # _35 = PHI <_37(4), _39(5)>
  _3 = COMPLEX_EXPR <_34, _35>;
  _32 = _34;
  _33 = _35;
  REALPART_EXPR <*z_5(D)> = _32;
  IMAGPART_EXPR <*z_5(D)> = _33;
  return;

}