Pulp 1.88.0
Cross-platform audio plugin framework
Loading...
Searching...
No Matches
pulp::signal::FastMath Struct Reference

Math helpers optimized for audio DSP. More...

#include <fast_math.hpp>

Static Public Member Functions

template<FastTrigProfile Profile>
static float sin_cycles (float phase_cycles) noexcept
 Computes sine for a normalized phase in [0, 1).
static float sin_cycles (float phase_cycles, FastTrigProfile profile) noexcept
 Runtime profile dispatch for setup, control, and non-inner-loop use.
static float tanh (float x)
 Fast tanh: an odd [5/4] rational, within 5.7e-5 of std::tanh for every float input, continuous, and bounded by ±1.
static float tanh_clamp (float x) noexcept
 tanh()'s argument clamp, propagating NaN.
static float tanh_accurate (float x)
 tanh within 3.0e-7 absolute of std::tanh when the compiler fuses multiply-adds, and within 3.6e-7 when it does not (-ffp-contract=off, MSVC, x86-64 without FMA), over the inputs the tests scan: every 64th non-negative finite float bit pattern, with odd symmetry covering the negatives.
static float tanh_accurate_rational (float x, float u)
 The unsaturated rational inside tanh_accurate(): x * P(u) / Q(u) with u = min(x^2, 81).
static float tanh_minimax (float x)
 Constrained-minimax [7/6] rational tanh: two more coefficients per polynomial than tanh()'s [5/4], and about 15x lower worst-case error.
static void tanh_minimax (const float *input, float *output, std::size_t count) noexcept
 Block form of tanh_minimax(); output may equal input.
template<typename T>
static T tanh_minimax_rational (T x) noexcept
 The rational inside tanh_minimax() for an argument already clamped to [-7, 7], generic over float and a GNU float vector so both forms round identically.
static float tanh_minimax_recursive (float x) noexcept
 tanh_minimax() for a serial recursion: the same fit, clamp and numerator scale, with each cubic evaluated in Estrin form (two independent halves joined by t^2), so the chain from x to the quotient is one multiply-add shorter than tanh_minimax_rational()'s Horner form.
static float sin (float x)
 Fast sin approximation (Bhaskara-style, max error ~0.001 for any x).
static float cos (float x)
 Fast cos approximation. Input in radians.
static float exp2 (float x)
 Computes 2^x with the platform's standard float exp2 semantics.
static float log2 (float x)
 Fast log2 approximation (max error ~0.007 for x > 0).
static float pow (float base, float exp)
 Fast pow(base, exp) via exp2(exp * log2(base)).
static float db_to_gain (float db)
 Fast dB to linear gain: 10^(dB/20) = 2^(dB * log2(10)/20).
static float gain_to_db (float gain)
 Fast linear gain to dB: 20 * log10(gain) = 20 * log2(gain) / log2(10).
static float log10_accurate (float x) noexcept
 log10(x) for finite, positive, normal x, within 4 ulps of the correctly rounded result (3.47 worst measured), with or without fused multiply-add.
static float exp10_accurate (float x) noexcept
 10^x within 1.2e-7 relative of the correctly rounded result, with or without fused multiply-add.
static float gain_to_db_accurate (float gain) noexcept
 20 log10(gain) for a finite, positive, normal gain (clamp to a floor first); within log10_accurate()'s bound.
static float db_to_gain_accurate (float db) noexcept
 10^(dB / 20), within exp10_accurate()'s bound.
static float rcp (float x)
 Fast reciprocal (1/x) via Newton-Raphson with SSE hint.
static float rsqrt (float x)
 Fast inverse square root (1/sqrt(x)).
static float clamp_unit (float x)
 Clamp to [-1, 1] without branching (for saturation).
static float soft_clip (float x)
 Soft clipping (polynomial saturation, no discontinuity).

Static Public Attributes

static constexpr float kTanhClamp = 5.24f
 tanh() clamps its argument to [-kTanhClamp, kTanhClamp], where its rational reaches 1.
static constexpr float kTanhMinimaxClamp = 7.0f

Detailed Description

Math helpers optimized for audio DSP.

RT contract: all functions are stateless scalar math helpers and allocate no memory.

Functions documented as approximations trade precision for speed. Functions without that label retain the standard-library numerical contract. Bounds quoted "with or without fused multiply-add" hold under either side of Pulp's contraction policy (tools/cmake/PulpFloatingPointContraction.cmake: Clang contracts within an expression, GCC does not).

Choosing a tanh — one family, four points on the cost/accuracy curve (worst absolute error over every float input; cost per element relative to tanh() in a 128-sample loop on an M5 Ultra, all measured in one run):

Function Worst error Cost Clamp
tanh() 5.7e-5 1x argument at 5.24, continuous
tanh_minimax() block 3.8e-6 1.17x at 7: slope-matched
tanh_minimax() scalar 3.8e-6 1.8x at 7: slope-matched
tanh_accurate() 3.0e-7 1.9x at 9, branch-free
std::tanh exact 10x none

The bare [7/6] Lambert Padé, which has no clamp, costs 1.01x but is not a tanh outside a small range: it is 1e-4 off by |x| = 5, exceeds 1 there and grows without bound (about x/28), so every bounded entry above pays for its clamp on purpose.

Prefer tanh_minimax() where tanh() would be used and the extra 15x of accuracy is worth its cost, tanh_minimax_recursive() for the same fit inside a per-sample recursion, and tanh_accurate() inside a resonant recursion whose loop gain amplifies saturator error.

float out = FastMath::tanh(input); // ~10x faster than std::tanh
float freq = FastMath::exp2(pitch);
float phase = FastMath::sin(angle); // ~5x faster than std::sin
float db = FastMath::log2(amplitude); // ~3x faster than std::log2
static float log2(float x)
Fast log2 approximation (max error ~0.007 for x > 0).
Definition fast_math.hpp:440
static float tanh(float x)
Fast tanh: an odd [5/4] rational, within 5.7e-5 of std::tanh for every float input,...
Definition fast_math.hpp:194
static float sin(float x)
Fast sin approximation (Bhaskara-style, max error ~0.001 for any x).
Definition fast_math.hpp:402
static float exp2(float x)
Computes 2^x with the platform's standard float exp2 semantics.
Definition fast_math.hpp:435

Member Function Documentation

◆ clamp_unit()

float pulp::signal::FastMath::clamp_unit ( float x)
inlinestatic

Clamp to [-1, 1] without branching (for saturation).

◆ cos()

float pulp::signal::FastMath::cos ( float x)
inlinestatic

Fast cos approximation. Input in radians.

◆ db_to_gain()

float pulp::signal::FastMath::db_to_gain ( float db)
inlinestatic

Fast dB to linear gain: 10^(dB/20) = 2^(dB * log2(10)/20).

◆ db_to_gain_accurate()

float pulp::signal::FastMath::db_to_gain_accurate ( float db)
inlinestaticnoexcept

10^(dB / 20), within exp10_accurate()'s bound.

◆ exp10_accurate()

float pulp::signal::FastMath::exp10_accurate ( float x)
inlinestaticnoexcept

10^x within 1.2e-7 relative of the correctly rounded result, with or without fused multiply-add.

Results beyond the float range saturate to +inf or 0 (through the subnormal range) as std::pow does, and NaN propagates.

10^x = 2^n * 10^d with n = round(x log2(10)) and d = x - n log10(2), |d| <= log10(2) / 2. The reduction is Cody-Waite in the decimal domain: log10(2) is split into a high part with few enough significant bits that n * hi is exact for every n used, so x - n * hi is exact (Sterbenz) and d carries no rounding from the product whether or not the compiler fuses multiply-adds. 10^d = e^(d ln 10) is the Taylor series through degree 7 (truncation below 6e-9 relative). 2^n is applied as two normal-range factors so n may run past the float exponent range and the final product over- or underflows naturally.

◆ exp2()

float pulp::signal::FastMath::exp2 ( float x)
inlinestatic

Computes 2^x with the platform's standard float exp2 semantics.

NaN propagates, positive infinity and finite overflow produce positive infinity, and negative infinity and finite underflow produce positive zero. With gradual underflow enabled, subnormal results and representable integer powers are preserved. Under Pulp's audio-callback ScopedFlushDenormals policy, the finite subnormal range may flush to positive zero exactly as the platform std::exp2 does in that ambient FP mode. This deliberately uses the standard implementation: a polynomial plus ldexp measured slower on supported Apple hardware while weakening these edge contracts.

◆ gain_to_db()

float pulp::signal::FastMath::gain_to_db ( float gain)
inlinestatic

Fast linear gain to dB: 20 * log10(gain) = 20 * log2(gain) / log2(10).

◆ gain_to_db_accurate()

float pulp::signal::FastMath::gain_to_db_accurate ( float gain)
inlinestaticnoexcept

20 log10(gain) for a finite, positive, normal gain (clamp to a floor first); within log10_accurate()'s bound.

◆ log10_accurate()

float pulp::signal::FastMath::log10_accurate ( float x)
inlinestaticnoexcept

log10(x) for finite, positive, normal x, within 4 ulps of the correctly rounded result (3.47 worst measured), with or without fused multiply-add.

Other inputs are outside the contract (callers clamp to a positive floor first).

x = m * 2^k with m in [sqrt(1/2), sqrt(2)); ln(m) = 2 atanh(s) with s = (m - 1) / (m + 1), |s| <= 0.1716, from the odd series 2s (1 + s^2/3 + s^4/5 + ...) truncated after s^8/9; the next term is below 2e-9 relative, under float resolution.

◆ log2()

float pulp::signal::FastMath::log2 ( float x)
inlinestatic

Fast log2 approximation (max error ~0.007 for x > 0).

◆ pow()

float pulp::signal::FastMath::pow ( float base,
float exp )
inlinestatic

Fast pow(base, exp) via exp2(exp * log2(base)).

◆ rcp()

float pulp::signal::FastMath::rcp ( float x)
inlinestatic

Fast reciprocal (1/x) via Newton-Raphson with SSE hint.

◆ rsqrt()

float pulp::signal::FastMath::rsqrt ( float x)
inlinestatic

Fast inverse square root (1/sqrt(x)).

◆ sin()

float pulp::signal::FastMath::sin ( float x)
inlinestatic

Fast sin approximation (Bhaskara-style, max error ~0.001 for any x).

Input in radians.

◆ sin_cycles() [1/2]

template<FastTrigProfile Profile>
float pulp::signal::FastMath::sin_cycles ( float phase_cycles)
inlinestaticnoexcept

Computes sine for a normalized phase in [0, 1).

The real-time profiles fold the already-bounded phase into a quarter cycle and avoid general range reduction. Callers with an arbitrary phase must wrap it before calling this function.

The two polynomial expressions are adapted from Lasse Schlör's "Fast MiniMax Polynomial Approximations of Sine and Cosine", pinned at https://github.com/publik-void/sin-cos-approximations/tree/d65178e684c7626b0fe7df6f261dbadc54403bce and used under the author's public permission: https://gist.github.com/publik-void/067f7f2fef32dbe5c27d6e215f824c91?permalink_comment_id=5556230#gistcomment-5556230

◆ sin_cycles() [2/2]

float pulp::signal::FastMath::sin_cycles ( float phase_cycles,
FastTrigProfile profile )
inlinestaticnoexcept

Runtime profile dispatch for setup, control, and non-inner-loop use.

◆ soft_clip()

float pulp::signal::FastMath::soft_clip ( float x)
inlinestatic

Soft clipping (polynomial saturation, no discontinuity).

◆ tanh()

float pulp::signal::FastMath::tanh ( float x)
inlinestatic

Fast tanh: an odd [5/4] rational, within 5.7e-5 of std::tanh for every float input, continuous, and bounded by ±1.

x·P(x²)/Q(x²), P quadratic and Q quadratic with Q(0) = 1, is the minimax fit of tanh in absolute error on [0, 5.24] under one equality constraint, f(5.24) = 1, so clamping the argument at ±5.24 (kTanhClamp) meets the saturation value with no step. Fit by direct minimisation of the maximum error over a dense grid, coefficients rounded to float. The clamp point balances the two error sources: the fit's equiripple error inside [0, 5.24] and the tail 1 - tanh(5.24) = 5.6e-5 beyond it are equal. The slope at the origin is 0.99976, not 1: the fit spends that freedom on the worst case. Over [-5, 5] this is within 5.7e-5 of tanh, where the [7/6] Lambert Padé of the same cost is 1.0e-4 off at 5 and unbounded beyond it.

The numerator is scaled by 1 - 2^-21, as in tanh_minimax_rational(), so every evaluation order stays strictly inside ±1; beyond the clamp the result is about ±0.9999996. Rounding can step the result down by up to 4.2e-7 (a few ulps) where the curve is nearly flat; otherwise it is monotone.

Measured in float over every input: worst absolute error 5.69e-5 with fused multiply-add and 5.69e-5 without, odd, |result| < 1. The clamp is one min and one max and propagates NaN (IEEE maximum/minimum on Clang, ordered selects elsewhere), so a NaN input returns NaN. No branches or calls, so loops over it vectorize: nine operations and a division, the same count as the unclamped [7/6] Padé.

◆ tanh_accurate()

float pulp::signal::FastMath::tanh_accurate ( float x)
inlinestatic

tanh within 3.0e-7 absolute of std::tanh when the compiler fuses multiply-adds, and within 3.6e-7 when it does not (-ffp-contract=off, MSVC, x86-64 without FMA), over the inputs the tests scan: every 64th non-negative finite float bit pattern, with odd symmetry covering the negatives.

Cheap enough for a saturator inside a per-sample feedback loop.

x * P(x^2) / Q(x^2) with both polynomials of degree 4 and P(0) = Q(0) = 1, fitted for minimax absolute error on [0, 9]; x^2 is clamped to 81 and the result to [-1, 1], so the result is exactly +/-1 for |x| >= 9.0001 (std::tanh itself rounds to 1 only above about 9.01). Properties the ladder relies on:

  • odd: tanh_accurate(-x) == -tanh_accurate(x) exactly, -0 stays -0;
  • unit slope at 0 and relative error below 3e-7 for small |x|, so the small-signal loop gain of a resonant filter is unchanged;
  • bounded for every non-NaN input (+/-inf give +/-1); NaN propagates, as with std::tanh;
  • monotone to within its error bound (float rounding can step down by a few ulps near saturation);
  • no branches or calls, so loops over it vectorize. Faster than std::tanh inside a ladder recursion.

◆ tanh_accurate_rational()

float pulp::signal::FastMath::tanh_accurate_rational ( float x,
float u )
inlinestatic

The unsaturated rational inside tanh_accurate(): x * P(u) / Q(u) with u = min(x^2, 81).

tanh_accurate() is this clamped to [-1, 1]; exposed so a vectorized caller can apply the same clamp with its own selects.

◆ tanh_clamp()

float pulp::signal::FastMath::tanh_clamp ( float x)
inlinestaticnoexcept

tanh()'s argument clamp, propagating NaN.

◆ tanh_minimax() [1/2]

void pulp::signal::FastMath::tanh_minimax ( const float * input,
float * output,
std::size_t count )
inlinestaticnoexcept

Block form of tanh_minimax(); output may equal input.

The clamp is a min/max pair rather than branches, so the whole evaluation is straight-line and runs four lanes at a time, four vectors per iteration (GNU vector extensions on Clang and GCC; a scalar loop the compiler may vectorise elsewhere). A clamped lane evaluates the rational at ±7, about 4e-7 short of ±1, so beyond the clamp this differs from the scalar form by that much and inside it not at all; both stay within the 3.8e-6 bound. The two clamp instructions cost a little over the unclamped [7/6] Padé, which has no clamp and grows without bound past |x| ~ 5. With auto-vectorisation disabled it holds its speed while the scalar forms slow down. The result for a NaN input is unspecified.

◆ tanh_minimax() [2/2]

float pulp::signal::FastMath::tanh_minimax ( float x)
inlinestatic

Constrained-minimax [7/6] rational tanh: two more coefficients per polynomial than tanh()'s [5/4], and about 15x lower worst-case error.

Truncating Lambert's continued fraction to the same [7/6] shape gives an error that is tiny near the origin and grows toward the clamp. This form spends the same coefficients on the worst case instead. They are the minimax fit of x·P(x²)/Q(x²) (P and Q cubic) to tanh on [0, 7] under three linear equality constraints: unit slope at the origin, f(7) = 1 and f'(7) = 0, so the clamp beyond |x| = 7 meets the curve with matching slope. Fit by bisection on the error bound over a linear-programming feasibility problem (the linearised rational minimax of Cheney & Loeb, Numer. Math. 3, 1961), coefficients rounded to float.

The numerator is then scaled by 1 − 2⁻²¹. Unscaled, rounding lets the rational overshoot 1 by up to two float ulps near the clamp, and how far depends on which multiply-adds a compiler contracts; the scale keeps every evaluation order strictly inside ±1 without an output clamp. Its cost is a 4e-7 step where the clamp takes over and a slope of 1 − 2⁻²¹ at the origin, both far below the fit error.

Measured in float over every representable input: absolute error below 3.8e-6 (the tail beyond the clamp contributes 1 − tanh(7) = 1.7e-6), odd, and bounded by ±1 exactly. The clamps are branches: inside a serial recursion such as a ladder they are predicted off the critical path, where min/max instructions would sit on it and cost 6-65% in a ladder (more under heavy drive, where a predicted branch skips the division). Over a buffer use the block overload below, which clamps branch-free. A NaN input returns NaN.

◆ tanh_minimax_rational()

template<typename T>
T pulp::signal::FastMath::tanh_minimax_rational ( T x)
inlinestaticnoexcept

The rational inside tanh_minimax() for an argument already clamped to [-7, 7], generic over float and a GNU float vector so both forms round identically.

◆ tanh_minimax_recursive()

float pulp::signal::FastMath::tanh_minimax_recursive ( float x)
inlinestaticnoexcept

tanh_minimax() for a serial recursion: the same fit, clamp and numerator scale, with each cubic evaluated in Estrin form (two independent halves joined by t^2), so the chain from x to the quotient is one multiply-add shorter than tanh_minimax_rational()'s Horner form.

It rounds differently from tanh_minimax(), by a few ulps, within the same 3.8e-6 bound (measured over every float, with and without fused multiply-add: 3.67e-6 and 3.70e-6).

Inside a polyphonic voice loop (oscillator call, envelope, ladder per voice), the Horner chain let the compiler's pre-register-allocation instruction scheduler emit an order in which consecutive voices stopped overlapping: the loop ran 2.6x slower than with the accurate saturator, with the core stalled on execution (processing-bottleneck share 82% against 49%; branch mispredicts and instruction delivery both under 1%). The Estrin form ran 0.71-1.0x the accurate saturator's time in every tested voice context, at -O2 and -O3. Over a buffer, where throughput rather than latency binds, Horner is cheaper by its one fewer multiply, which is why tanh_minimax() keeps it.

Member Data Documentation

◆ kTanhClamp

float pulp::signal::FastMath::kTanhClamp = 5.24f
staticconstexpr

tanh() clamps its argument to [-kTanhClamp, kTanhClamp], where its rational reaches 1.

◆ kTanhMinimaxClamp

float pulp::signal::FastMath::kTanhMinimaxClamp = 7.0f
staticconstexpr

The documentation for this struct was generated from the following file:
  • /home/runner/work/pulp/pulp/core/signal/include/pulp/signal/fast_math.hpp