[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]

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

--- Comment #13 from anlauf at gcc dot gnu.org ---
(In reply to Steve Kargl from comment #12)
> The optimization level is irrelevant.  gfortran unilaterally
> uses -fcx-fortran-rules, and there is no way to disable this
> option to user the slower, but stricter, evaluation.  One
> will always get complex division computed by
> 
> a+ib   a + b(d/c)     b - a(d/c) 
> ---- = ---------- + i ------------  |c| > |d|
> c+id   c + d(d/c)     c + d(d/c)
> 
> and similar for |d| > |c|.
> 
> There are a few problems with this. d/c can trigger an invalid underflow
> exception.  If d == c, you then have numerators of a + b and b - a, you
> can get a invalid overflow for a = huge() and b > 1e291_8.

I am wondering how slow an algorithm would be that scales numerator
and denominator by respective factors that are powers of 2, e.g.

e_num = 2. ** -max (exponent (a), exponent (b))
e_den = 2. ** -max (exponent (c), exponent (d))

The modulus of scaled values would be <= 1, even for any of a,... being huge().
Of course this does not address underflows that could occur during scaling,
or denormalized numbers, which are numerically irrelevant for the result.

Is there anything else wrong with this approach?