I wonder, how did you came to this test? Its result looks paradoxical at first glance, but makes sense after thinking about it more.

That was just a coverage test, that triggers underflow inratio. The hard part, as you noted, was to figure out why it's true. I've just checked the exact answer in the Diofant (probably will work in the SymPy too), e.g. for the real part:

In [43]: x = Symbol('x', extended_real=True)In [44]: y = Symbol('y', positive=True, finite=True)In [45]: expr = re(expand((x + I*oo)*(y-I/y)))In [46]: expr.subs({x: Dummy(real=True)})Out[46]: ∞In [47]: expr.subs({x: Dummy(positive=True)})Out[47]: ∞In [48]: limit(expr, x, -oo)Out[48]: ∞

BTW, Mathematica does more brave things, but hardly it's correct:

In[14]:= Hold[(x + I*Infinity)*(2^10-I/2^10)]  (* no assumptions on x ! *)                                 10    IOut[14]= Hold[(x + I Infinity) (2   - ---)]                                       10                                      2In[15]:= ReleaseHold[%]         (1 + 1048576 I) InfinityOut[15]= ------------------------           Sqrt[1099511627777]In[16]:= ReleaseHold[%14/.{x->1-I*Infinity}]Infinity::indet: Indeterminate expression -I Infinity + I Infinity encountered.Out[16]= Indeterminate

However, it looks I should restore remark about corner cases. Not all are handled by this code (have no chance, because one component computed as non-nan). Here is an example:

>>>complex('(1+infj)')/complex(2**1000,2**-1000)(nan+infj)

This is correctly handled by C code both for gcc and clang, but it turns they use a different algorithm for tdiv (not Smith's method, rather something like Priest's algorithm; see full _Cdivd() code before handling nan+nanj case). Maybe it worth to adopt this, but probably it's out of the scope of this pr.

Edit:
An example, where current algorithm fails (bad imaginary part) is$(2^{1023}+i 2^{-1023})/(2^{677}+i 2^{-677})\approx 2^{346}-i 2^{-1008}$ (taken from M.Baudin: "A Robust Complex Division in Scilab"). gcc v12.2 (and, probably, clang), that follows to example in C standard, does produce the correct answer.

#include<complex.h>#include<stdio.h>#include<math.h>/* CPython's version */complextdiv_smith(complexa,complexb){constdoubleabs_breal=fabs(creal(b));constdoubleabs_bimag=fabs(cimag(b));if (abs_breal >=abs_bimag) {if (abs_breal!=0.0) {constdoubleratio=cimag(b) /creal(b);constdoubledenom=creal(b)+cimag(b)*ratio;returnCMPLX((creal(a)+cimag(a)*ratio) /denom,                         (cimag(a)-creal(a)*ratio) /denom);        }    }elseif (abs_bimag >=abs_breal) {constdoubleratio=creal(b) /cimag(b);constdoubledenom=creal(b)*ratio+cimag(b);returnCMPLX((creal(a)*ratio+cimag(a)) /denom,                     (cimag(a)*ratio-creal(a)) /denom);    }returnCMPLX(NAN,NAN);}complextdiv_smith2(complexa,complexb){constdoubleabs_breal=fabs(creal(b));constdoubleabs_bimag=fabs(cimag(b));if (abs_breal >=abs_bimag) {if (abs_breal!=0.0) {constlongdoubleratio= (longdouble)cimag(b) / (longdouble)creal(b);constdoubledenom=creal(b)+cimag(b)*ratio;returnCMPLX((creal(a)+cimag(a)*ratio) /denom,                         (cimag(a)-creal(a)*ratio) /denom);        }    }elseif (abs_bimag >=abs_breal) {constlongdoubleratio= (longdouble)creal(b) / (longdouble)cimag(b);constdoubledenom=creal(b)*ratio+cimag(b);returnCMPLX((creal(a)*ratio+cimag(a)) /denom,                     (cimag(a)*ratio-creal(a)) /denom);    }returnCMPLX(NAN,NAN);}intmain(){complexa=CMPLX(0x1p+1023,0x1p-1023);complexb=CMPLX(0x1p677,0x1p-677);complexz;z=a /b;printf("%e%+ei\n",creal(z),cimag(z));z=tdiv_smith(a,b);printf("%e%+ei\n",creal(z),cimag(z));z=tdiv_smith2(a,b);printf("%e%+ei\n",creal(z),cimag(z));return0;}

$ cc -std=c11 a.c -lm && ./a.out1.433437e+104-3.645561e-304i1.433437e+104+0.000000e+00i1.433437e+104-3.645561e-304i

It's quite odd to see NaNs disappearing in an arithmetic operation on operands that include NaNs. However, that seems unavoidable if we want to follow the C standard Annex G semantics, since its rules imply thatcomplex(inf, nan) should be regarded as a complex infinity, and then the division of a finite complex number by a complex infinity should be a complex zero (so can't have any NaN components).

It's quite odd to see NaNs disappearing in an arithmetic operation on operands that include NaNs.

I think that could be reasonably interpreted only as shown above, i.e.nan in components means "anything from the extended real line".

since its rules imply that complex(inf, nan) should be regarded as a complex infinity

I don't think there is such thing in the C standard. It has two types of infinities: when one component is finite, nonzero and other is infinite; other type (with few cases) is a directed infinity (phase is finite and absolute value is infinite). I.e.(1+infj) and(inf+infj). That brings an interesting clash with CAS like Mathematica, that has no notion like the first type if infinities. E.g.ArcSin[Infinity] is just-I Infinity c.f.cmath.asin(). Seempmath/mpmath#793

But CPython already adopts Annex G semantics in this respect, viacmath.isinf() andcmath.isfinite().

@mdickinson, what do you think about changing Smith's algorithm to something more closer to C standard?

@skirpichev

I don't think there is such thing in the C standard.

Take a look at Annex G. Quoting from §G.3 "Conventions" in C99, paragraph 1:

A complex or imaginary value with at least one infinite part is regarded as an infinity (even if its other part is a NaN).

The language is slightly changed in the most recent publicly available draft of C23 (n3096.pdf):

A complex or imaginary value with at least one infinite part is regarded as an infinity (even if its other part is a quiet NaN).