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::math::fast_sinf(w);  // SVF coefficient
//
//   // Exponential modulation (LFO → filter cutoff):
//   float mod_cutoff = cutoff * sbl::dsp::math::fast_exp2f(lfo * depth_oct);
 
#ifndef SBL_DSP_MATH_FAST_MATH_HPP_
#define SBL_DSP_MATH_FAST_MATH_HPP_
 
#include <sbl/dsp/math/constants.hpp>
 
namespace sbl::dsp::math {
 
/// 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  = PI;
    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
    constexpr float C1 = 0.6931472f;  // ln(2)
    constexpr float C2 = 0.2402265f;  // ln(2)^2 / 2 (approx)
    constexpr float C3 = 0.0558011f;  // ln(2)^3 / 6 (approx)
    constexpr float C4 = 0.00898f;    // ln(2)^4 / 24 (approx)
    float p = 1.0f + f * (C1 + f * (C2 + f * (C3 + f * C4)));
 
    // 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;
}
 
/// Fast tanh approximation using [3,2] Padé approximant
///
/// tanh(x) ≈ x · (27 + x²) / (27 + 9·x²) for |x| ≤ 3
///
/// Same rational function as SoftLimiter::limit() in soft_limiter.hpp, but with
/// input clamping and NaN guard. SoftLimiter is deliberately unbounded
/// (approaches x/9 for large x); this function saturates to ±1.
///
/// Accurate to < 1% relative error for |x| < 3, which covers the
/// operating range of the Moog ladder filter (audio signals + moderate
/// drive + feedback). Inputs outside ±3 are clamped to ±1.
///
/// Cost: ~8 cycles on M7 FPU (2 MUL, 2 FMA, 1 VDIV).
///
/// @param x Input value (any range, clamped for |x| > 3)
/// @return Approximation of tanh(x)
inline float fast_tanhf(float x) {
    if (x < -3.0f) return -1.0f;
    if (x >  3.0f) return  1.0f;
    if (x != x) return 0.0f;  // NaN → 0 (defense in depth)
    // [3,2] Padé: coefficients 27 = 3^3 and 9 = 3^2
    float x2 = x * x;
    return x * (27.0f + x2) / (27.0f + 9.0f * x2);
}
 
/// A unipolar signal to a frequency, exponentially: min_hz at 0, min_hz · 2^octaves at 1.
/// The one su-to-Hz map for every cutoff and damping port.
inline float su_to_hz(float su, float min_hz, float octaves) {
    return min_hz * fast_exp2f(su * octaves);
}
 
/// 1 − e^(−x) by inverse Taylor, for a one-pole coefficient from a time
/// constant or a cutoff (x = 3/(τ·rate), 4.6/n, ω). Cold path; no expf.
inline constexpr float one_minus_exp_neg(float x) {
    const float denom = 1.0f + x * (1.0f + x * (0.5f + x * (1.0f / 6.0f)));
    return 1.0f - 1.0f / denom;
}
 
} // namespace sbl::dsp::math
 
#endif // SBL_DSP_MATH_FAST_MATH_HPP_