DspMath.h

// DSPark -- Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi -- MIT License
 
#pragma once
 
#include <algorithm>
#include <array>
#include <cmath>
#include <concepts>
#include <numbers>
 
namespace dspark {
 
// ============================================================================
// Concepts
// ============================================================================
 
template <typename T>
concept FloatType = std::floating_point<T>;
 
namespace detail {
// Shared coefficients of the existing reduced-range odd sine polynomial.
// Nonlinear offline kernels also use its exact slope and residual bounds.
template <FloatType T> inline constexpr std::array<T, 5> fastSinPolynomial{
    T(0.99999997408724855), T(-0.16666646026660671), T(0.0083328727116396326),
    T(-0.00019799239565814083), T(2.5871610835732768e-6)};
} // namespace detail
 
// ============================================================================
// Constants
// ============================================================================
 
template <FloatType T> inline constexpr T pi       = std::numbers::pi_v<T>;
 
template <FloatType T> inline constexpr T twoPi    = T(2) * std::numbers::pi_v<T>;
 
template <FloatType T> inline constexpr T invTwoPi = T(1) / twoPi<T>;
 
template <FloatType T> inline constexpr T halfPi   = std::numbers::pi_v<T> / T(2);
 
template <FloatType T> inline constexpr T sqrt2    = std::numbers::sqrt2_v<T>;
 
template <FloatType T> inline constexpr T invSqrt2 = T(1) / std::numbers::sqrt2_v<T>;
 
// ============================================================================
// Decibel / Gain Conversions
// ============================================================================
 
template <FloatType T>
[[nodiscard]] inline T decibelsToGain(T dB, T minusInfinityDb = T(-100)) noexcept
{
    // exp(dB * ln(10)/20): the same accuracy as pow(10, dB/20) (measured
    // 2.97e-7 vs 2.78e-7 relative in float over +-60 dB) at half the cost.
    return dB <= minusInfinityDb ? T(0) : std::exp(dB * T(0.11512925464970228420));
}
 
template <FloatType T>
[[nodiscard]] inline T gainToDecibels(T gain, T minusInfinityDb = T(-100)) noexcept
{
    // 20/ln(10) * ln(gain): log() is about twice as fast as log10() in float
    // at the same accuracy (glibc's log10f is not the optimised kernel).
    return gain > T(0) ? std::max(minusInfinityDb, T(8.6858896380650365530) * std::log(gain))
                       : minusInfinityDb;
}
 
// ============================================================================
// Interpolation and Mapping
// ============================================================================
 
template <FloatType T>
[[nodiscard]] inline T mapRange(T value, T inMin, T inMax, T outMin, T outMax) noexcept
{
    if (inMin == inMax) return outMin; // Safety: prevent NaN/Inf in audio pipeline
    return outMin + (outMax - outMin) * ((value - inMin) / (inMax - inMin));
}
 
template <FloatType T>
[[nodiscard]] inline T moveTowards(T from, T to, T maxDelta) noexcept
{
    // Return the target itself once within reach: from + (to - from) rounds
    // off it in about 10% of float cases, leaving a settled ramp one ulp away.
    const T delta = to - from;
    return std::abs(delta) <= maxDelta ? to : from + std::clamp(delta, -maxDelta, maxDelta);
}
 
// ============================================================================
// Fast Approximations
// ============================================================================
 
template <FloatType T>
[[nodiscard]] inline T fastTanh(T x) noexcept
{
    // std::clamp translates to SIMD min/max intrinsics (no branching).
    // Maintains continuity avoiding step-aliasing at the bounds.
    x = std::clamp(x, T(-3), T(3));
    
    const auto x2 = x * x;
    const auto x4 = x2 * x2;
    // Pade [5,4] approximant: max error < 0.05% for |x| <= 3
    return x * (T(945) + T(105) * x2 + x4) / (T(945) + T(420) * x2 + T(15) * x4);
}
 
template <FloatType T>
[[nodiscard]] inline T fastPow10(T x) noexcept
{
    // log2(10) folded into one constant: rounding it once is more accurate
    // than multiplying log2(e) * ln(10) at runtime, and it saves a multiply
    // (compilers may not reassociate x * a * b without fast-math).
    constexpr T kLog2Of10 =
        static_cast<T>(std::numbers::ln10_v<long double> / std::numbers::ln2_v<long double>);
    return std::exp2(x * kLog2Of10);
}
 
template <FloatType T>
[[nodiscard]] inline T fastExp(T x) noexcept
{
    return std::exp2(x * std::numbers::log2e_v<T>);
}
 
template <FloatType T>
[[nodiscard]] inline T fastTan(T x) noexcept
{
    constexpr T limit = T(1.520);  // ~pi/2 - 0.05
    if (std::abs(x) > limit) return std::tan(x);
    const T x2 = x * x;
    const T x4 = x2 * x2;
    // Pade [5,4]: tan(x) ~= x * (945 - 105*x^2 + x^4) / (945 - 420*x^2 + 15*x^4)
    return x * (T(945) - T(105) * x2 + x4) / (T(945) - T(420) * x2 + T(15) * x4);
}
 
template <FloatType T>
[[nodiscard]] inline T fastSin(T x) noexcept
{
    // Range-reduce to [-pi, pi] with a two-term Cody-Waite split of 2*pi so the
    // subtraction stays accurate in float even for arguments several periods out.
    const T k = std::floor(x * invTwoPi<T> + T(0.5));
    x -= k * T(6.28125);                  // high part (exactly representable)
    x -= k * T(0.0019353071795864769253); // low part of 2*pi
 
    // Fold into [-pi/2, pi/2] where the polynomial converges fast:
    // sin(x) = sin(pi - x) for x > pi/2 (and the odd mirror for x < -pi/2).
    if (x > halfPi<T>)       x = pi<T> - x;
    else if (x < -halfPi<T>) x = -pi<T> - x;
 
    const T x2 = x * x;
    // Endpoint-constrained Remez minimax coefficients for sin on [-pi/2, pi/2]
    // (equiripple error 3.73e-9, p(pi/2) == 1). Near-Taylor coefficients of
    // the same degree leave ~3e-6 of error: a -110 dB harmonic floor.
    constexpr auto c = detail::fastSinPolynomial<T>;
    return x * (c[0] + x2 * (c[1] + x2 * (c[2] + x2 * (c[3] + x2 * c[4]))));
}
 
template <FloatType T>
[[nodiscard]] inline T fastCos(T x) noexcept
{
    return fastSin(x + halfPi<T>);
}
 
template <FloatType T>
[[nodiscard]] inline T fastLog(T x) noexcept
{
    int e = 0;
    T m = std::frexp(x, &e);             // x = m * 2^e, m in [0.5, 1)
    // Normalise mantissa into [sqrt(0.5), sqrt(2)) so the series is centred.
    if (m < T(0.70710678118654752440)) { m *= T(2); --e; }
 
    // Classic atanh form (Cephes): ln(m) = 2*atanh(s), s = (m-1)/(m+1).
    // |s| <= 0.1716 so the odd series converges below 1e-7 with 4 terms.
    const T s  = (m - T(1)) / (m + T(1));
    const T s2 = s * s;
    const T lnm = T(2) * s * (T(1)
                + s2 * (T(1.0 / 3.0)
                + s2 * (T(1.0 / 5.0)
                + s2 *  T(1.0 / 7.0))));
 
    return lnm + static_cast<T>(e) * std::numbers::ln2_v<T>; // + e * ln2
}
 
// ============================================================================
// Utility
// ============================================================================
 
template <FloatType T>
[[nodiscard]] inline T wrapPhase(T phase) noexcept
{
    // floor() correctly handles negative values, preventing branching
    T wrapped = phase - twoPi<T> * std::floor(phase * invTwoPi<T>);
    // Rounding in the k * twoPi product can land wrapped just outside
    // [0, 2*pi) on either side: large phases where the product overshoots
    // the input leave a slightly negative result, and tiny negative phases
    // can round the sum up to exactly twoPi. Enforce the documented
    // half-open range. Order matters: the += branch can itself round up to
    // exactly twoPi, which the second branch then folds to 0.
    if (wrapped < T(0))      wrapped += twoPi<T>;
    if (wrapped >= twoPi<T>) wrapped -= twoPi<T>;
    return wrapped;
}
 
} // namespace dspark