fast_math.hpp

Go to the documentation of this file.

// sbl/dsp/fast_math.hpp — Fast analytical approximations (Audio Stack — Atoms)
//
// Lightweight math functions for hot audio paths where libm functions
// like sinf() are too expensive. All functions are pure arithmetic —
// no lookup tables, no flash cost.
//
// fast_sinf(x): 5th-order Taylor series for sin(x), valid for x in
//   [0, π/2]. Max error < 0.00016 at x = π/2. ~5 cycles on M7 FPU
//   vs ~500 cycles for newlib sinf().
//
// fast_exp2f(x): Fast 2^x for exponential modulation (1V/oct-style).
//   Integer part via IEEE 754 exponent shift, fractional part via
//   3rd-order minimax polynomial. ~8 cycles on M7 FPU. Accurate to
//   ~20 bits for |x| <= 4 (±4 octaves).
//
// Usage:
//   float w = pi * fc / (2.0f * fs);
//   float f = 2.0f * sbl::dsp::fast_sinf(w);  // SVF coefficient
//
//   // Exponential modulation (LFO → filter cutoff):
//   float mod_cutoff = cutoff * sbl::dsp::fast_exp2f(lfo * depth_oct);
 
#pragma once
 
namespace sbl::dsp {
 
/// Fast sine approximation for x in [0, π/2]
///
/// Uses 5th-order Taylor series: x - x³/6 + x⁵/120
///
/// Accuracy:
///   x = 0.0    → error = 0
///   x = 0.589  → error < 0.00016 (SVF max at 18 kHz / 48 kHz)
///   x = π/2    → error ≈ 0.00045 (theoretical max in domain)
///
/// @param x Input angle in radians, must be in [0, π/2]
/// @return Approximation of sin(x)
inline float fast_sinf(float x) {
    float x2 = x * x;
    return x * (1.0f - x2 * (1.0f / 6.0f - x2 * (1.0f / 120.0f)));
}
 
/// Fast tan(π·f) for normalized frequency f ∈ [0, 0.497]
///
/// [5,4] Padé approximant of tan(x) evaluated at x = π·f:
///
///   tan(x) ≈ x · (945 - 105x² + x⁴) / (945 - 420x² + 15x⁴)
///
/// Unlike a polynomial, the Padé rational function correctly models the
/// pole at f = 0.5 (Nyquist). The denominator goes to zero at x = π/2,
/// matching the true singularity of tan.
///
/// Accuracy vs true tan(π·f):
///   f = 0.10 (4.8 kHz)  → error < 0.001%
///   f = 0.35 (16.8 kHz) → error < 0.3%
///   f = 0.45 (21.6 kHz) → error < 2.4%
///   f = 0.497 (23.9 kHz) → error < 0.1%
///
/// The previous 5th-order polynomial (from MI stmlib) diverged badly
/// above f ≈ 0.35, giving 45% error at f = 0.45 — causing audible
/// artifacts (secondary resonant peaks) in the ZDF SVF.
///
/// Cost: ~8 FMA + 1 VDIV ≈ 20–25 cycles on M7 FPU (vs ~5 for the old
/// polynomial, vs ~50–100 for newlib tanf).
///
/// @param f Normalized frequency (freq_hz / sample_rate), must be < 0.497
/// @return Approximation of tan(π·f)
inline float fast_tan_pif(float f) {
    constexpr float pi  = 3.14159265f;
    constexpr float pi2 = pi * pi;
    constexpr float pi4 = pi2 * pi2;
    float f2 = f * f;
    float f4 = f2 * f2;
    float num = 945.0f - 105.0f * pi2 * f2 + pi4 * f4;
    float den = 945.0f - 420.0f * pi2 * f2 + 15.0f * pi4 * f4;
    return pi * f * num / den;
}
 
 
/// Fast 2^x approximation for exponential modulation
///
/// Decomposes x into integer and fractional parts. The integer part is
/// applied by shifting the IEEE 754 exponent field (exact). The fractional
/// part uses a 3rd-order minimax polynomial (accurate to ~20 bits).
///
/// This is the "exponential converter" primitive — the analog equivalent
/// of the circuit that makes 1V/oct work in a VCF or VCO.
///
/// Accuracy (vs std::exp2f):
///   |x| <= 1  → max relative error < 0.02%
///   |x| <= 4  → max relative error < 0.05%
///   |x| > 16  → clamped (returns 0 for x < -16)
///
/// @param x Exponent (e.g., ±2.0 for ±2 octave modulation)
/// @return Approximation of 2^x
inline float fast_exp2f(float x) {
    // Clamp to prevent overflow/underflow
    if (x < -16.0f) return 0.0f;
    if (x >  16.0f) x = 16.0f;
 
    // Decompose into integer and fractional parts
    // Use truncation toward negative infinity
    int i = static_cast<int>(x);
    float f = x - static_cast<float>(i);
    if (f < 0.0f) { f += 1.0f; --i; }
 
    // 4th-order polynomial for 2^f, f in [0, 1)
    // Remez-style minimax coefficients for improved accuracy over Taylor
    float p = 1.0f + f * (0.6931472f + f * (0.2402265f
            + f * (0.0558011f + f * 0.00898f)));
 
    // Apply integer exponent via IEEE 754 bit manipulation
    union { float fv; int32_t iv; } v;
    v.fv = p;
    v.iv += i << 23;
    return v.fv;
}
 
} // namespace sbl::dsp