Python Complex Number Arithmetic Implementation

This article explores how Python handles complex number arithmetic under the hood using its native complex data type. It covers the underlying C structures in CPython, the algorithms used for basic arithmetic operations like multiplication and division, and how Python preserves numerical precision while mitigating overflow and underflow risks.

The Underlying Data Structure

In standard Python (CPython), the built-in complex type is implemented at the C level as PyComplexObject. This object wraps a low-level C struct named Py_complex, defined as follows:

typedef struct {
    double real;
    double imag;
} Py_complex;

A Python complex instance stores two 64-bit IEEE 754 double-precision floating-point numbers: one for the real component and one for the imaginary component. Like Python floats, complex numbers are immutable; any arithmetic operation creates a new PyComplexObject rather than modifying existing operands.

Addition and Subtraction

Addition and subtraction are straightforward component-wise operations. When evaluating (a + bj) + (c + dj), Python directly computes:

• Real component: a + c
• Imaginary component: b + d

At the C level, CPython adds or subtracts the corresponding double fields directly:

Py_complex c_sum;
c_sum.real = a.real + b.real;
c_sum.imag = a.imag + b.imag;

These operations inherit the standard IEEE 754 floating-point rounding and edge-case behaviors (such as handling infinities and NaNs).

Complex Multiplication

Mathematically, complex multiplication expands as:

\[(a + bj) \times (c + dj) = (ac - bd) + (ad + bc)j\]

CPython implements this logic in its internal function _Py_c_prod. The computation is carried out using double-precision arithmetic:

• Real component: (a.real * b.real) - (a.imag * b.imag)
• Imaginary component: (a.real * b.imag) + (a.imag * b.real)

If any resulting component overflows the range of a standard double, the value saturates to infinity (inf), following standard floating-point rules.

Complex Division and Numerical Stability

Directly calculating complex division using the naive formula:

\[\frac{a + bj}{c + dj} = \frac{(ac + bd) + (bc - ad)j}{c^2 + d^2}\]

can easily result in premature overflow or underflow when squaring \(c\) and \(d\), even if the final result fits comfortably within floating-point limits.

To prevent this, CPython's internal _Py_c_diff function implements Smith's method (a variation of scaled complex division). It normalizes the operation depending on the relative magnitudes of the denominator's components:

1. If \(|c| \ge |d|\):

   • Compute the ratio \(r = d / c\).
   • Compute the denominator scale \(den = c + d \cdot r\).
   • Real part: \((a + b \cdot r) / den\)
   • Imaginary part: \((b - a \cdot r) / den\)

2. If \(|d| > |c|\):

   • Compute the ratio \(r = c / d\).
   • Compute the denominator scale \(den = c \cdot r + d\).
   • Real part: \((a \cdot r + b) / den\)
   • Imaginary part: \((b \cdot r - a) / den\)

By dividing by the larger component first, the ratio \(r\) remains within \([-1, 1]\), preventing intermediate overflow caused by squaring large values.

Absolute Value (Magnitude)

Calling abs(z) on a complex number computes its Euclidean norm: \(\sqrt{a^2 + b^2}\).

Instead of naive squaring, CPython delegates this operation to the standard C library function hypot(a, b). The hypot implementation internally rescales \(a\) and \(b\) to ensure intermediate squares do not overflow to inf or underflow to zero unnecessarily.

Advanced Mathematical Functions

While basic arithmetic (+, -, *, /, **, abs) is supported directly by the complex type via Python's operator slots, transcendental operations (such as square roots, logarithms, and trigonometric functions) are not methods on the object. Instead, they are provided by the cmath standard library module, which implements C99-compliant complex algorithms optimized for branch cuts and domain boundaries.