Module complex 

C’s complex arithmetic, as the generated code calls it.

What C’s operators actually do

Two rules decide everything in this module, and neither is what a naive reading of “the usual arithmetic conversions” suggests.

• A real operand stays real. C99 6.3.1.8 says the common type of a real and a complex operand is the complex one, but the result it prescribes is the same one the componentwise operation gives — and GCC and Clang both compute it componentwise, which is observable: 3.0 * z for z = (-0.0, -0.0) is (-0.0, -0.0) componentwise and (+0.0, -0.0) through the full product, because 3·(−0) − 0·(−0) is +0. The one exception is real / complex, which really is the full division with the imaginary part of the numerator taken as zero. Each mixed form therefore has a function of its own here: mul_real_f64, add_real_f64, sub_real_f64, real_sub_f64, div_real_f64 and real_div_f64, with the f32 twins beside them.

• Infinities are recovered. C99 Annex G.5.1 (which GCC implements by default — -fno-cx-limited-range) says that a product or a quotient involving an infinity must be infinite, even where the naive formula produces NaN + iNaN out of an ∞ − ∞ or an ∞ · 0. mul_f64 and div_f64 therefore compute the cheap answer first and only fix it up when both parts came out NaN, which is what makes the common case cost four multiplications and an addition. Only that answer and the NaN test are inlined into the caller, as GCC and Clang inline the product before their call of __muldc3: the fix-up is a separate #[cold] function, never inlined, in both operations. Inlined, it made LLVM pack a loop’s real and imaginary parts into one vector for the sake of the path never taken and put the shuffles on the one that is, costing a z = z * z + c loop a quarter of its time.

The algorithms

The product — in both widths — is the schoolbook one, (ac − bd) + i(ad + bc), followed by Annex G.5.1’s recovery: an infinite operand is reduced to a signed one or zero, a NaN in the other operand is reduced to a signed zero, and the product is recomputed scaled by infinity so that the sign of each part survives.

The quotient takes a different route in each width, because what is cheap and exact in one is not available in the other.

• div_f32 uses the closed form in the wider format: c² + d² cannot overflow or underflow in f64 for any f32 operands, so (ac + bd)/(c² + d²) and (bc − ad)/(c² + d²) are computed in f64 and rounded once.
• div_f64 has no wider format to hand, so it uses Smith’s algorithm — divide through by whichever of the divisor’s parts is larger, so that c² + d² cannot overflow when a single part can — with Baudin and Smith’s refinement for a subnormal ratio, where a · (d/c) has already lost bits that the equal (a/c) · d still has.

Both then apply the matching recovery: a zero divisor gives a signed infinity, an infinite dividend over a finite divisor gives an infinity, and a finite dividend over an infinite divisor gives a signed zero.

All of it is written from the standard’s description rather than taken from any implementation, and tests/complex.rs checks the result against the host’s own C compiler: every combination of 0, −0, ±1, ±2, ±10⁸, ±10⁻⁸, ±10¹⁵⁰, ±10⁻¹⁵⁰, ±∞ and NaN — a hundred and thirty thousand pairs — agrees exactly, and the product agrees for the 10³⁰⁰ and subnormal magnitudes too.

Where it is not the last word

libgcc goes further on the f64 quotient: it prescales both operands by powers of two before dividing, which buys accuracy in the last corner — operands whose exponents are hundreds apart, where an intermediate of Smith’s algorithm underflows even after the refinement above. GCC’s own gcc.c-torture/execute/ieee/cdivchkd.c is four such quotients, and two of its four are beyond what is here. C leaves the accuracy of complex arithmetic implementation-defined (Annex G.6), so this is a quality gap and not a conformance one; doc/gcc-torture.md records it.

Functions

add_real_f32

A float _Complex value plus a real one, which touches the real part only.

add_real_f64

A double _Complex value plus a real one, which touches the real part only.

conj_f32

The conjugate of a float _Complex value: conj, and GNU’s ~z.

conj_f64

The conjugate of a double _Complex value: conj, and GNU’s ~z.

div_f32

The quotient of two float _Complex values (C99 6.5.5, Annex G.5.1).

div_f64

The quotient of two double _Complex values (C99 6.5.5, Annex G.5.1).

div_real_f32

A float _Complex value divided by a real one, which C computes componentwise.

div_real_f64

A double _Complex value divided by a real one, which C computes componentwise.

eq_f32

Whether two float _Complex values are equal: C99 6.5.9p3, both parts equal.

eq_f64

Whether two double _Complex values are equal: C99 6.5.9p3, both parts equal.

mul_f32

The product of two float _Complex values (C99 6.5.5, Annex G.5.1).

mul_f64

The product of two double _Complex values (C99 6.5.5, Annex G.5.1).

mul_real_f32

A float _Complex value times a real one, which C computes componentwise.

mul_real_f64

A double _Complex value times a real one, which C computes componentwise.

narrow_f64

double _Complex narrowed to float _Complex.

ne_f32

Whether two float _Complex values differ — the negation of eq_f32 , and != in C.

ne_f64

Whether two double _Complex values differ — the negation of eq_f64 , and != in C.

nonzero_f32

Whether a float _Complex value is non-zero, which is what C’s conversion to _Bool and a controlling expression ask (C99 6.3.1.2).

nonzero_f64

Whether a double _Complex value is non-zero, which is what C’s conversion to _Bool and a controlling expression ask (C99 6.3.1.2).

proj_f32

The projection of a float _Complex value onto the Riemann sphere: C99 7.3.9.5’s cproj.

proj_f64

The projection of a double _Complex value onto the Riemann sphere: C99 7.3.9.5’s cproj.

real_add_f32

A real value plus a float _Complex one; see real_mul_f32 for why the order matters.

real_add_f64

A real value plus a double _Complex one; see real_mul_f64 for why the order matters.

real_div_f32

A real value divided by a float _Complex one.

real_div_f64

A real value divided by a double _Complex one.

real_mul_f32

A real value times a float _Complex one: the same product, with the operands in source order so that the generated code evaluates them where they were written.

real_mul_f64

A real value times a double _Complex one: the same product, with the operands in source order so that the generated code evaluates them where they were written.

real_sub_f32

A real value minus a float _Complex one, which negates the imaginary part.

real_sub_f64

A real value minus a double _Complex one, which negates the imaginary part.

sub_real_f32

A float _Complex value minus a real one, which touches the real part only.

sub_real_f64

A double _Complex value minus a real one, which touches the real part only.

widen_f32

float _Complex widened to double _Complex.