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;
}
