Skip to content

Cfloat — Implementation

File Structure

cfloat/
├── tiny_cfloat.h        (24 lines — type definition + declarations)
└── tiny_cfloat.c        (70 lines — implementation)

Dependencies: tiny_math_config.h (via header), <math.h> (for sqrtf, cosf, sinf, atan2f, logf, fabsf).

Implementation Summary

All functions are pure C, pass-by-value, and return-by-value. There is no dynamic allocation or platform-specific dispatch — the module is fully portable.

Function Lines Key Implementation Detail
tiny_cf_add 3 Element-wise addition
tiny_cf_sub 3 Element-wise subtraction
tiny_cf_mul 3 Standard \((a+bi)(c+di) = (ac-bd) + i(ad+bc)\)
tiny_cf_div 20 Smith's robust algorithm (see below)
tiny_cf_conj 2 Negate imaginary part
tiny_cf_abs 2 sqrtf(re² + im²), guarded against zero
tiny_cf_arg 2 atan2f(im, re)
tiny_cf_from_polar 3 r·cosf(θ) + i·r·sinf(θ)
tiny_cf_sqrt 8 Principal branch via \(\sqrt{\frac{\|z\|+z_r}{2}} + i \cdot \operatorname{sgn}(z_i) \sqrt{\frac{\|z\|-z_r}{2}}\)
tiny_cf_log 7 \(\ln\|z\| + i \cdot \arg(z)\)

Key Algorithms

Division: Smith's Robust Method

Naive division computes \(|b|^2 = b_r^2 + b_i^2\) for the denominator. This can overflow for large coefficients. Smith's method scales the computation to avoid intermediate overflow:

tiny_cfloat_t tiny_cf_div(tiny_cfloat_t a, tiny_cfloat_t b)
{
    float abs_re = fabsf(b.re), abs_im = fabsf(b.im);
    float mag = (abs_re > abs_im) ? abs_re : abs_im;

    // Near-zero denominator guard
    if (mag < TINY_MATH_MIN_POSITIVE_INPUT_F32) {
        return sentinel_large_value(a);
    }

    if (abs_re >= abs_im) {
        float k = b.im / b.re;
        float den = b.re + k * b.im;
        r.re = (a.re + k * a.im) / den;
        r.im = (a.im - k * a.re) / den;
    } else {
        float k = b.re / b.im;
        float den = b.im + k * b.re;
        r.re = (a.re * k + a.im) / den;
        r.im = (a.im * k - a.re) / den;
    }
    return r;
}

Square Root: Principal Branch

tiny_cfloat_t tiny_cf_sqrt(tiny_cfloat_t a)
{
    if (a.re == 0.0f && a.im == 0.0f) return zero;
    float mag = tiny_cf_abs(a);
    r.re = sqrtf((mag + a.re) * 0.5f);
    r.im = (a.im >= 0.0f)
        ?  sqrtf((mag - a.re) * 0.5f)
        : -sqrtf((mag - a.re) * 0.5f);
    return r;
}

The sign of the imaginary part follows C99 convention: \(\operatorname{sgn}(\Im(\sqrt{z})) = \operatorname{sgn}(\Im(z))\), with the branch cut on the negative real axis.

Usage in Higher Modules

tiny_cfloat_t is the fundamental type used by:

  • decomp — Householder reflections require complex arithmetic for eigenvalue extraction
  • eigen — Eigenvalues are stored as tiny_cfloat_t[] (real eigenvalues have \(\Im(\lambda)=0\))
  • modal — Modal analysis operates on complex-valued eigenvalues and mode shapes
  • iterative — Arnoldi/Lanczos iterations produce complex Ritz values
// Example: eigenvalues from Francis QR are returned as complex numbers
tiny_cfloat_t eigs[10];
tiny_eig_general_f32(A, 10, eigs, workspace, ws_size);

for (int i = 0; i < 10; i++) {
    printf("λ%d = %f + i·%f\n", i, eigs[i].re, eigs[i].im);
}