Complex arithmetic semantics

Three types, one plane¶

Where this sits: the C article covers the standard lineage that made these types optional in C11, and the reference page lists the headers involved. This page is about what Annex G actually requires of the arithmetic.

C99 gave the language two families of floating types where Fortran and C++ have one. Alongside the real types there are the complex types, and — only in an Annex G implementation — the imaginary types. An imaginary value is a single real number understood as a coefficient of i; it is not a complex number whose real part happens to be zero, and the difference is not cosmetic. 0.0 + iy and iy occupy the same point of the plane but do not behave the same way under multiplication, because one of them has a real part that participates in every partial product and the other has no real part at all.

The imaginary types are optional even where Annex G applies: an implementation signals them with __STDC_IEC_559_COMPLEX__, and portable code that wants the effects described here must be prepared for a translation unit where _Imaginary does not exist. The complex types themselves became optional in C11, guarded by __STDC_NO_COMPLEX__.

The three arithmetic type families of an Annex G implementation. Source: ISO/IEC 9899:2017, 6.2.5 and G.2.

Back to top

The type algebra of the operators¶

The multiplicative and additive operators disagree about which combination of operand types stays out of the complex types, and the disagreement is the whole reason imaginary types earn their keep. Under * and /, imaginary is closed with real: two imaginary operands give a real result and a mixed pair gives an imaginary one, so a chain of scalings never widens to two components. Under + and - it is the mixed pair that widens, because a real part and an imaginary part added together are exactly what a complex number is.

Read the two grids against each other rather than in turn. The four corners agree; the two off-diagonal cells swap; and the middle cell — imaginary against imaginary — is real under multiplication and imaginary under addition. Every practical consequence of having imaginary types in the language follows from those three differences.

Back to top

The operator formula tables¶

For every combination except complex against complex, G.5.1 and G.5.2 do not describe the result approximately: they give the exact expression the implementation must evaluate, and the floating-point exception behaviour is whatever that expression raises. Nothing is left to the implementation, and in particular nothing may be computed by promoting the operands to complex and running the general formula — that would introduce products with a zero part, and products with a zero part are exactly how an infinity turns into a NaN.

The plus-or-minus in Figure 5 is the operator itself, not an unspecified sign: the same table serves + and -, and a leading ±u means the real part is u for addition and -u for subtraction.

Back to top

What the tables cost in real operations¶

The formula tables are not only a semantic specification; they are a cost model. Every cell names the exact real operations the implementation performs, so the price of losing an imaginary type to an implicit conversion is countable. A rotation written as z * I with an imaginary I is two moves and a sign flip. The same rotation written against a complex constant is four multiplications, two additions, and a NaN-recovery branch that has to run whenever both parts of the result come out NaN.

Counts are of arithmetic operations named in the formula, not of instructions: a fused multiply-add may merge two of them, which is why both worked examples in the annex open with #pragma STDC FP_CONTRACT OFF.

Back to top

The infinity properties¶

Paragraph 4 of G.5.1 states four properties that hold for all real, imaginary and complex operands — including the complex-by-complex case the formula tables decline to cover. Footnote 383 is explicit that this is not a restatement of the tables: the tables already imply these results where they apply, and paragraph 4 extends them to every case the tables leave open. In practice this is the contract that forces the recovery code in the next section to exist.

Multiplication is symmetric here, so the first row covers both orders of operand. Division is not, and the three division rows are three genuinely different rules.

The four infinity properties. “An infinity” means a value with at least one infinite part, which by G.3 includes a value whose other part is a NaN. Source: ISO/IEC 9899:2017, G.5.1 paragraph 4.

Neither grid pins a single cell in a NaN row or column, and that is deliberate rather than an omission: paragraph 4 is about infinities, and a NaN operand is governed by the ordinary IEC 60559 rules for the real operations the formula tables name. The gap is where the trouble starts, because the ordinary rules turn ∞ × 0 into a NaN and there is no complex-level rule that puts the infinity back.

Back to top

Multiplication: recovering the infinities¶

Take z = ∞ + i∞ and w = 2 + i0. Both are perfectly good values, the product is unambiguously an infinity, and the textbook formula computes it as a NaN in both parts. The reason is that ∞ × 0 is a NaN under IEC 60559 and two of the four partial products contain one. Once a NaN enters a partial product it survives the subtraction and the addition, and the infinity that was there in the operand is gone.

The reference implementation¶

G.5.1 Example 1 gives an implementation of double _Complex multiplication that meets the paragraph 4 contract. Its shape is worth reading closely: the recovery path is guarded by a single isnan(x) && isnan(y) test. && short-circuits and an ordinary product has a non-NaN real part, so the finite path pays for one isnan and nothing else — which is exactly the cost the annex claims for it.

#include <math.h>
#include <complex.h>

/* Multiply z * w ... */
double complex _Cmultd(double complex z, double complex w)
{
    #pragma STDC FP_CONTRACT OFF
    double a, b, c, d, ac, bd, ad, bc, x, y;

    a = creal(z); b = cimag(z);
    c = creal(w); d = cimag(w);
    ac = a * c; bd = b * d;
    ad = a * d; bc = b * c;
    x = ac - bd; y = ad + bc;

    if (isnan(x) && isnan(y)) {
        /* Recover infinities that computed as NaN + iNaN ... */
        int recalc = 0;
        if (isinf(a) || isinf(b)) {              /* z is infinite */
            /* "Box" the infinity and change NaNs in the other factor to 0 */
            a = copysign(isinf(a) ? 1.0 : 0.0, a);
            b = copysign(isinf(b) ? 1.0 : 0.0, b);
            if (isnan(c)) c = copysign(0.0, c);
            if (isnan(d)) d = copysign(0.0, d);
            recalc = 1;
        }
        if (isinf(c) || isinf(d)) {              /* w is infinite */
            c = copysign(isinf(c) ? 1.0 : 0.0, c);
            d = copysign(isinf(d) ? 1.0 : 0.0, d);
            if (isnan(a)) a = copysign(0.0, a);
            if (isnan(b)) b = copysign(0.0, b);
            recalc = 1;
        }
        if (!recalc && (isinf(ac) || isinf(bd) ||
                        isinf(ad) || isinf(bc))) {
            /* Recover infinities from overflow by changing NaNs to 0 ... */
            if (isnan(a)) a = copysign(0.0, a);
            if (isnan(b)) b = copysign(0.0, b);
            if (isnan(c)) c = copysign(0.0, c);
            if (isnan(d)) d = copysign(0.0, d);
            recalc = 1;
        }
        if (recalc) {
            x = INFINITY * (a * c - b * d);
            y = INFINITY * (a * d + b * c);
        }
    }
    return x + I * y;
}

Source: ISO/IEC 9899:2017, G.5.1 Example 1, reformatted. The annex notes that this achieves the required treatment of infinities at the cost of one isnan test in ordinary finite cases, and that it is less than ideal in that undue overflow and underflow may occur.

The three recovery conditions of Example 1. copysign does the boxing, so the magnitude collapses to 1 or 0 while every sign survives — which is what lets the final INFINITY * reconstruct the correct quadrant. Source: ISO/IEC 9899:2017, G.5.1 Example 1.

Back to top

Division: scaling the denominator¶

Complex division has the same NaN problem as multiplication and one more of its own. The denominator c*c + d*d squares the operand, so it overflows whenever either part of the divisor exceeds the square root of the largest finite value, and flushes to zero whenever both parts fall below the square root of the smallest positive value. For double that means a divisor with a part above roughly 1.3 × 10¹⁵⁴ produces an infinite denominator and a zero quotient, and a divisor with both parts below roughly 2 × 10^-162 produces a zero denominator and an infinite quotient — in both cases from operands that were nowhere near the edge of the format.

Example 2 fixes this by scaling. It takes the binary exponent of the larger divisor part with logb, divides both parts by that power of two with scalbn, does the division on the scaled values, and scales the two result parts back. Because both scalings are by powers of two they are exact, which is why the annex notes that scaling with scalbn rather than with a division gives better roundoff characteristics. The numerator is left unscaled, and the annex says so plainly: this code does not defend against overflow and underflow there either.

#include <math.h>
#include <complex.h>

/* Divide z / w ... */
double complex _Cdivd(double complex z, double complex w)
{
    #pragma STDC FP_CONTRACT OFF
    double a, b, c, d, logbw, denom, x, y;
    int ilogbw = 0;

    a = creal(z); b = cimag(z);
    c = creal(w); d = cimag(w);

    logbw = logb(fmax(fabs(c), fabs(d)));
    if (isfinite(logbw)) {
        ilogbw = (int)logbw;
        c = scalbn(c, -ilogbw); d = scalbn(d, -ilogbw);
    }
    denom = c * c + d * d;
    x = scalbn((a * c + b * d) / denom, -ilogbw);
    y = scalbn((b * c - a * d) / denom, -ilogbw);

    /* Recover infinities and zeros that computed as NaN + iNaN;   */
    /* the only cases are nonzero/zero, infinite/finite, and        */
    /* finite/infinite ...                                          */
    if (isnan(x) && isnan(y)) {
        if ((denom == 0.0) && (!isnan(a) || !isnan(b))) {
            x = copysign(INFINITY, c) * a;
            y = copysign(INFINITY, c) * b;
        }
        else if ((isinf(a) || isinf(b)) && isfinite(c) && isfinite(d)) {
            a = copysign(isinf(a) ? 1.0 : 0.0, a);
            b = copysign(isinf(b) ? 1.0 : 0.0, b);
            x = INFINITY * (a * c + b * d);
            y = INFINITY * (b * c - a * d);
        }
        else if ((logbw == INFINITY) && isfinite(a) && isfinite(b)) {
            c = copysign(isinf(c) ? 1.0 : 0.0, c);
            d = copysign(isinf(d) ? 1.0 : 0.0, d);
            x = 0.0 * (a * c + b * d);
            y = 0.0 * (b * c - a * d);
        }
    }
    return x + I * y;
}

Source: ISO/IEC 9899:2017, G.5.1 Example 2, reformatted.

The three recovery branches, one per division property in G.5.1 paragraph 4. The comment in the annex is exact: these are the only cases that reach the branch, because every other route to NaN + iNaN is a genuine NaN. Source: ISO/IEC 9899:2017, G.5.1 Example 2.

Note that logbw, not c or d, is what the third branch tests. logb of an infinity is an infinity, so a single comparison identifies an infinite divisor without inspecting either part — and it also explains why the scaling is guarded by isfinite(logbw) in the first place.

Back to top

CX_LIMITED_RANGE and FP_CONTRACT¶

Everything above is what the implementation owes you by default. Two pragmas let you give some of it back, and both worked examples in the annex begin by switching one of them off.

Both are block-scoped and both revert with DEFAULT. Turning CX_LIMITED_RANGE on is a licence, not an instruction: the implementation may ignore it entirely.

The two pragmas that govern how literally the formula tables are followed. Footnote 383 attaches the infinity properties to the state where CX_LIMITED_RANGE is off. Sources: ISO/IEC 9899:2017, 6.5 and 7.3.4, and G.5.1 footnote 383.

Turning CX_LIMITED_RANGE on discards the infinity properties of G.5.1 paragraph 4 along with the cost of enforcing them. A product of two infinities may then be NaN + iNaN, and code that tests isinf(creal(z)) || isinf(cimag(z)) to detect overflow will stop detecting it.

Back to top

<complex.h> under Annex G¶

G.6 changes one macro definition, and the change is the practical payoff of the whole annex. Clause 7.3 defines I as _Complex_I, a complex constant whose real part is a zero. G.6 redefines it as _Imaginary_I, an imaginary constant with no real part at all. Notwithstanding 7.1.3, a program may undefine and then redefine imaginary; I itself is the macro that quietly decides whether x + I*y is correct.

This is why C11 added the CMPLX, CMPLXF and CMPLXL macros to <complex.h>. On an implementation without imaginary types there is no way to write a complex value with a given pair of parts using arithmetic alone, because x + I*y goes through the complex I and loses infinities and signed zeros. CMPLX(x, y) is defined to build the value directly, and it is the portable spelling.

/* Portable on any C11 implementation, whether or not it has imaginary types. */
double complex z = CMPLX(1.0, INFINITY);        /* 1.0 + i INF, exactly     */

/* Correct only where I is _Imaginary_I, that is, under Annex G. */
double complex w = 1.0 + I * INFINITY;          /* NaN + i INF otherwise    */

/* Detecting the two feature macros. */
#ifdef __STDC_NO_COMPLEX__
#  error "this translation unit needs the complex types"
#endif
#ifdef __STDC_IEC_559_COMPLEX__
   /* Annex G applies: imaginary types exist and the tables above are binding. */
#endif

Sources: ISO/IEC 9899:2017, 6.10.8.3 (conditional feature macros), 7.3.9.3 (the CMPLX macros), and G.6 paragraph 1.

Back to top

Branch cuts and the sign of zero¶

Every inverse of a periodic or multivalued function needs a cut in the plane where the principal value jumps. G.6 paragraph 3 states the rule that makes those jumps usable rather than merely survivable: the functions are continuous onto both sides of their branch cuts, taking into account the sign of zero. A signed zero in the imaginary part is not a rounding artefact to be normalised away — it is which side of the cut you are approaching from, and the annex gives csqrt(−2 ± i0) = ±i√2 as the worked case.

This is the practical reason a complex library must never normalise -0.0 to 0.0. Two values that compare equal under == sit on opposite sides of a cut and produce results that differ in sign. It is also why the additive table in Figure 5 matters: assembling x + i(-0.0) componentwise preserves the sign, and computing it as x + I*(-0.0) with a complex I does not.

Back to top

Identities and symmetry¶

Five of the <complex.h> functions have no special-value catalogue of their own. G.6 paragraph 7 specifies them implicitly, by an identity in terms of a function that does have one, and their special values follow from that identity rather than from a list. Two more, cabs and carg, are specified as real functions whose special cases live in Annex F.

The symmetry reduction¶

Paragraph 8 states the rule that lets each catalogue below be as short as it is. For a function f satisfying f(conj(z)) = conj(f(z)), the specifications for the upper half-plane imply those for the lower half-plane. If f is also even or odd, the first quadrant implies the other three. Every catalogue in G.6.1 to G.6.4 opens with exactly that line, so a list of a dozen cases is really a list of up to four dozen.

Which part of the plane each catalogue actually enumerates. The grids in the next section are filled out to the whole plane using these rules, so a cell for a negative real part is derived rather than quoted where the annex states only the positive case. Source: ISO/IEC 9899:2017, G.6 paragraph 8 and the opening line of each of G.6.1.1 to G.6.4.2.

Back to top

The special-value catalogue¶

Ten of the <complex.h> functions get an explicit list of special cases in G.6. The lists are exhaustive over the argument classes that matter, but they are lists: reading one tells you what happens at a point, not what the function does to the plane. Laid out as a grid over the real-part class and the imaginary-part class, the same information shows its shape, and the shape is consistent enough to be worth naming.

Every grid below shares the same axes. Columns are the class of the real part x, from −∞ on the left to NaN on the right; rows are the class of the imaginary part y, from +∞ at the top down to NaN at the bottom, so the top three rows read like the upper half-plane. Cells marked ordinary have no clause in the annex and take the mathematical value. A cell marked not specified is one the annex genuinely leaves open, not one omitted here. Cells derived from the symmetry rules of the previous section are marked in the notes table that follows the grids.

cacos¶

The inverse cosine. Its own entry in G.6.1; casin and catan are given by identity in G.6 paragraph 7 rather than enumerated. Specified in G.6.1.1.

cacosh¶

The inverse hyperbolic cosine, and the one function here whose cut is not symmetric about the origin. Specified in G.6.2.1.

casinh¶

The inverse hyperbolic sine. Odd, so the annex states only the first quadrant and the other three follow. Specified in G.6.2.2.

catanh¶

The inverse hyperbolic tangent, and one of the two functions in the annex whose catalogue can raise divide-by-zero. Specified in G.6.2.3.

ccosh¶

The hyperbolic cosine. Even, and one of the three functions whose infinite results are written with cis(y) = cos(y) + i sin(y). Specified in G.6.2.4.

csinh¶

The hyperbolic sine. Odd, and otherwise the mirror of ccosh. Specified in G.6.2.5.

ctanh¶

The hyperbolic tangent. Bounded, so its infinite arguments give finite results — including the memorable 1 + i0·sin(2y). Specified in G.6.2.6.

cexp¶

The exponential, with the largest catalogue in the annex: thirteen enumerated cases. Specified in G.6.3.1.

clog¶

The natural logarithm, and the other of the two: it raises divide-by-zero at the origin, from either side of the cut. Specified in G.6.3.2.

csqrt¶

The square root. Its top row is the only one in the annex that is constant across every real-part class, NaN included. Specified in G.6.4.2.

Exceptions and unspecified signs¶

Two kinds of caveat appear on the grids above and neither is visible in a return value: whether the call raises a floating-point exception, and whether a sign in the result is unspecified. Both matter to a caller — the first because a later fetestexcept will see it, the second because a test on the sign of a part is not portable. This table is generated from the same cells as the grids.

Every row is a case where the return value alone does not tell you what happened.

All 63 flagged cells across the ten catalogues. “Opt. invalid” is the annex’s optionally raises the invalid floating-point exception: an implementation may raise it and may not, so a caller must neither rely on it nor be surprised by it. Source: ISO/IEC 9899:2017, G.6.1.1 to G.6.4.2.

Back to top

What the catalogue adds up to¶

Summed over the ten grids, the annex is less a list of oddities than a statement about where the complex plane stops being ordinary. Of the 196 argument classes the ten grids cover, 30 take the ordinary mathematical value and every one of the remaining 166 is pinned to a specific result — and 4 combinations are left open by the annex altogether.

Enumerated cases are the bullets of each clause, excluding the symmetry identity that opens it; the two right-hand columns are counted from the grids above. The totals differ because one clause can cover several argument classes — csqrt(x + i∞) alone covers a whole row. Sources: ISO/IEC 9899:2017, G.6.1.1 to G.6.4.2, and the grids in the previous section.

Back to top

Type-generic macros and imaginary arguments¶

G.7 extends <tgmath.h> so that a type-generic macro accepts an imaginary argument, and it does not simply widen it to complex. The macro expands to an expression whose type is real, imaginary or complex depending on the function, which means a chain of type-generic calls on an imaginary value can stay one component wide from end to end.

The three outcomes. Six macros collapse an imaginary argument to a real result, eight keep it imaginary, and everything else widens. Source: ISO/IEC 9899:2017, G.7 paragraph 1.

Every right-hand side is a call to a real function, which is what makes the extension cheap: an implementation of tan for an imaginary argument is tanh on the coefficient, with the result reinterpreted as imaginary. No complex arithmetic is generated at all.

Back to top

How this reached the standard¶

The practical position today is that the complex types are widely available and Annex G is not. An implementation may define __STDC_IEC_559_COMPLEX__ only if it meets everything on this page, and most do not: imaginary types in particular are rarely provided, which is why CMPLX exists and why the third row of Figure 14 is a portability hazard rather than a curiosity. Code that needs the infinity behaviour without the annex has to reimplement the recovery paths of G.5.1 itself — which is exactly what the two worked examples are for.

Back to top

Sources¶