Fast Math

Overview

Audio DSP calls transcendental functions constantly: an oscillator needs a sine every sample, a compressor converts to and from decibel (log10 / pow10) on every gain update, a soft clipper computes a tanh per sample. The standard-library versions are correctly rounded to the last bit, at a cost in cycles that a per-sample loop rarely needs to pay: the result is quantized, summed with noise and turned into sound. Q ships a family of fast approximations that trade a little accuracy for a large speed-up.

Every function comes in two tiers, distinguished by name:

  • fast_*: the accurate tier. Error is small enough to be inaudible in almost any audio use (a few parts in 10^-5 for the bounded functions). Use it by default.

  • faster_*: the cheap tier. It drops a polynomial term or two, so the error grows to the low percent range, in exchange for the fewest possible operations. Use it where the result feeds something forgiving, a modulation source or a rough envelope, and the cycles matter.

The approximations are IEEE-754 bit tricks and short polynomials, from Paul Mineiro’s fastapprox. They take and return float; there are no double overloads.

The error ripples as the input sweeps: each approximation is a smooth curve fitted to the true function, running slightly high then slightly low between the points where it crosses the exact value. fast_* fits closely, faster_* loosely, but neither introduces discontinuities, so both are safe to modulate. The worst-case sizes are tabulated in Accuracy and Speed.

Most of these functions are only valid over a restricted domain, listed per function below. The trigonometric approximations in particular assume the input is already range-reduced; feeding them an angle outside the stated range gives a wrong answer, not a clamped one.

fast_sqrt is an alias for std::sqrt, which lowers to a single hardware instruction on every target Q supports and is both exact and the fastest option. It is kept in the family so existing call sites need not change.
The fast_* trigonometric functions do no range reduction. fast_sin and fast_cos require -pi <= x <= pi; fast_tan requires -pi/2 <= x <= pi/2. Outside those ranges the polynomial diverges. Reduce the angle yourself (for example with a phase accumulator, which wraps by construction) before calling.

Include

#include <q/support/base.hpp>

Declaration

namespace cycfi::q
{
   // Exponential and logarithm
   inline float   fast_exp(float x);      inline float   faster_exp(float x);
   inline float   fast_log(float x);      inline float   faster_log(float x);
   inline float   fast_log2(float x);     inline float   faster_log2(float x);
   inline float   fast_log10(float x);    inline float   faster_log10(float x);
   inline float   fast_pow2(float x);     inline float   faster_pow2(float x);
   inline float   fast_pow10(float x);    inline float   faster_pow10(float x);

   // Trigonometric  (sin/cos: [-pi, pi];  tan: [-pi/2, pi/2])
   inline float   fast_sin(float x);      inline float   faster_sin(float x);
   inline float   fast_cos(float x);      inline float   faster_cos(float x);
   inline float   fast_tan(float x);      inline float   faster_tan(float x);

   // Hyperbolic tangent  (valid over all reals; saturates to +/-1)
   inline float   fast_tanh(float x);     inline float   faster_tanh(float x);
   constexpr float fast_rational_tanh(float x);   // Pade, -3 <= x <= 3

   // Reciprocal and division
   inline float   fast_inverse(float val);
   inline float   fast_div(float a, float b);

   // Square root (exact: alias for std::sqrt)
   inline float   fast_sqrt(float x);

   // Taylor-series exp, order 3..9  (small x)
   constexpr float fast_exp3(float x);    // ... through ...
   constexpr float fast_exp9(float x);

   // Fast integer RNG
   inline int     fast_rand();            // result in [0, 0x7FFF]
}

Expressions

Notation

x, a, b

float arguments.

Exponential and Logarithm

The logarithms require x > 0. The exponentials are valid over all reals; underflow is handled (a large negative input returns 0, not a denormal or NaN).

Expression Semantics Return Type

fast_exp(x)

e raised to x.

float

fast_log(x)

Natural logarithm of x.

float

fast_log2(x)

Base-2 logarithm of x.

float

fast_log10(x)

Base-10 logarithm of x (used by decibel).

float

fast_pow2(x)

2 raised to x.

float

fast_pow10(x)

10 raised to x (used by decibel).

float

Each row has a faster_ counterpart with the same signature and semantics but larger error (see Accuracy and Speed).

fast_pow10 computes pow2(x * log2(10)) using the exact log2(10), rather than routing through a generic pow that would scale by an approximate constant. Folding in the exact constant costs the same single multiply yet is markedly more accurate; against std::pow, fast_pow10 RMSE drops about 5x and faster_pow10 about 2x.

Trigonometric

Range-reduced only. fast_sin / fast_cos assume -pi <= x <= pi; fast_tan assumes -pi/2 <= x <= pi/2.

Expression Semantics Return Type

fast_sin(x)

Sine of x, x in [-pi, pi].

float

fast_cos(x)

Cosine of x, x in [-pi, pi].

float

fast_tan(x)

Tangent of x, x in [-pi/2, pi/2].

float

Each has a faster_ counterpart.

Hyperbolic Tangent

A staple nonlinearity for saturation and soft clipping.

Expression Semantics Return Type

fast_tanh(x)

Hyperbolic tangent, valid over all reals; saturates to +/-1 naturally, so no input clamp is needed. Exp-based; branchless and divide-free.

float

faster_tanh(x)

Cheaper, less accurate fast_tanh.

float

fast_rational_tanh(x)

Pade rational approximation, constexpr. Requires -3 <= x <= 3; outside that it does not saturate. Prefer fast_tanh unless you need a compile-time constant.

float

Reciprocal, Division, and Square Root

Expression Semantics Return Type

fast_inverse(v)

Approximate reciprocal 1/v by negating the IEEE-754 exponent. Coarse (a first-order estimate), but avoids a divide.

float

fast_div(a, b)

a * fast_inverse(b). Inherits `fast_inverse’s coarseness.

float

fast_sqrt(x)

Square root, x >= 0. Exact: an alias for std::sqrt.

float
fast_inverse (and therefore fast_div) is a low-accuracy estimate, good to only a few percent. Use it for control-rate math where the imprecision washes out, not where you need a faithful quotient.

Polynomial and Utility

Expression Semantics Return Type

fast_exp3(x) …​ fast_exp9(x)

constexpr Taylor-series exp, orders 3 through 9. Accurate only for small x; higher order widens the usable range. Use when a constexpr value is required.

float

fast_rand()

Fast integer pseudo-random number in [0, 0x7FFF], from a linear congruential recurrence with a fixed internal seed.

int
The bit-trick functions (fast_exp, fast_log*, fast_pow*, the trig and fast_tanh families, fast_inverse) are inline, not constexpr; they rely on runtime type punning. Only fast_rational_tanh and fast_exp3..fast_exp9 are constexpr.

Accuracy and Speed

Accuracy is the worst-case error over each function’s working domain: inaudible for the fast_* tier, coarse for faster_*. The figures below were measured with clang -O3 on an Apple-silicon Mac, scalar path; they shift with compiler, flags and target.

Maximum error over the working domain
Function fast_* faster_*

log2

~1.5e-4 (output bits)

~0.05 (output bits)

exp

~0.006% (relative)

~4% (relative)

sin

~4e-5 (absolute)

~9e-4 (absolute)

tanh

~3e-5 (absolute)

~2e-2 (absolute)

Speed depends on how the calls are arranged. Independent calls pipeline, and the CPU runs many at once (throughput); chained calls each wait for the last (latency). For these functions the two differ by up to an order of magnitude: the large savings are in throughput, and inside a dependency chain the pipeline latency dominates and the saving shrinks. Measure on your own target; the build-only micro-benchmarks under test/benchmark/ are for that.

decibel_bench.cpp reports both metrics for the log10 / pow10 pair used by decibel:

From decibel_bench (ns per call; clang -O3, Apple silicon)
Function Metric std fast_* faster_*

log10

throughput

5.8

0.31

0.16

latency

20

14

12

pow10

throughput

2.5

0.46

0.21

latency

23

19

11

In throughput fast_log10 is about twenty times cheaper than std::log10; in a dependency chain, about 1.4 times. log2_bench.cpp and sin_bench.cpp measure log2 and sin with a dependency-bound loop, nearer the latency column: std::sin runs about 12 ns against roughly 2 ns for both fast_sin and faster_sin, and std::log2 about 8 ns against 2.5 ns.

On a modern desktop FPU the faster_* tier is often no quicker than fast_*: both already collapse to a handful of instructions, so dropping a polynomial term is lost in pipeline latency. faster_* pulls ahead on cheaper hardware, where every multiply counts, and under vectorization. fast_* over the libc call is the saving that holds on every target. Rebuild any of these numbers on your own target by running the benchmarks.

Example

A cheap saturating soft clipper, one tanh per sample:

float drive(float x, float amount)
{
   return q::fast_tanh(x * amount);
}

Converting a linear amplitude to decibel and back in a gain stage, where log10 / pow10 would otherwise dominate the cost:

float g_db = 20.0f * q::fast_log10(amplitude);   // linear -> dB
// ... adjust g_db ...
float g = q::fast_pow10(g_db * 0.05f);           // dB -> linear

A range-reduced sine oscillator. fast_sin does no wrapping of its own, so the phase is accumulated in [-pi, pi) and folded back each sample to stay inside the valid domain:

float phase = 0.0f;                                  // in [-pi, pi)
float incr  = 2.0f * float(q::pi) * 440.0f / sps;    // radians per sample

// ... per sample:
float s = q::fast_sin(phase);
phase += incr;
if (phase >= float(q::pi))
   phase -= 2.0f * float(q::pi);                      // fold back into range
For production oscillators prefer the phase / phase_iterator accumulator, whose fixed-point counter wraps by construction, or the table-based Sine Wave Oscillator. The snippet above is here to make the domain restriction concrete.