This talk

Design custom math approximation functions
Focus
- Chebyshev polynomials
Also
- Taylor polynomials, Padé approximants
Tools
- Python & C++
- Claude assistance

About me

Ryan Robinson
- Independent developer
- Southern Oregon
Building apps & plug-ins since 2017
- AutoPad
  - iOS & Android
- Audio plug-ins at Sketch Audio, LLC
Spent five years teaching high school math

The plan

Motivation
Build a \sin x approximation
- Taylor, Padé
- Chebyshev polynomials
  - Refinement
Other functions
- \tan x
- \log_2 x, 2^x
- \tanh x

Along the way
- Error metrics
- Benchmarks
- Floating-point representation
- SIMD compatibility

Problem

In audio, we use a lot of special math functions
- Oscillators: sin, cos
- Filters: tan
- Nonlinearities: tanh
- Unit conversions: exp, log

Standard implementations
- Must handle any input
- Full precision expected
- Performance secondary to other requirements

Other considerations
- Inlinability
- Vectorizability

Solution

Find a polynomial that approximates f on [a, b]
May also use a rational function
Benefits
- Pure arithmetic
- Can control precision
Considerations
- Precondition: inputs must stay in [a, b]

\begin{aligned} p(x) &= p_0 + p_1 x + p_2 x^2 + \cdots + p_n x^n \\[14pt] R(x) &= \frac{p(x)}{q(x)} = \frac{p_0 + p_1 x + \cdots + p_n x^n}{q_0 + q_1 x + \cdots + q_m x^m} \end{aligned}

Building a \sin x approximation

Taylor polynomials

Match f and its first n derivatives at a

p_n(x) = \htmlClass{hl sum c-blue}{\sum_{k=0}^{n}} \htmlClass{hl coef c-blue}{\frac{f^{(k)}(a)}{k!}}\,\htmlClass{hl power c-blue}{(x - a)^k}

For \sin x about 0

\htmlClass{hl series c-blue}{\sin x \approx x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!}}

Differentiate c_3 x^3:

\begin{aligned} & c_3\, x^{\textcolor{#eb6834}{3}} \\ \tfrac{d}{dx} \;\to\;\; & \textcolor{#eb6834}{3}\, c_3\, x^{\textcolor{#eb6834}{2}} \\ \tfrac{d}{dx} \;\to\;\; & 3 \cdot \textcolor{#eb6834}{2}\, c_3\, x^{\textcolor{#eb6834}{1}} \\ \tfrac{d}{dx} \;\to\;\; & 3 \cdot 2 \cdot \textcolor{#eb6834}{1}\, c_3\, x^{\textcolor{#eb6834}{0}} = \htmlClass{hl fact c-blue}{3!}\, c_3 \end{aligned}

In general, at x = 0:

p^{(k)}(0) = \htmlClass{hl fact c-blue}{k!}\, c_k

Match f:

c_k = \frac{f^{(k)}(0)}{\htmlClass{hl fact c-blue}{k!}}

The idea with Taylor polynomials: match the value of the function and its first n derivatives at some central point a. Here's the formula. [click] Our polynomial is a sum [click] of powers of (x − a), [click] whose coefficients are the k-th derivative at a, divided by k!. [click] (Sidebar.) Where the factorials come from. Take one term, c3 x³, and differentiate it: each time the exponent comes down in front and drops by one, 3, then 2, then 1, so three derivatives leave 3·2·1·c3 = 3! c3. In general: write p(x) = c0 + c1 x + c2 x² + … and ask that its k-th derivative at 0 equal sin's k-th derivative at 0. Differentiating c_k x^k k times gives k·(k−1)·…·1·c_k = k!·c_k; every lower power has differentiated away, and every higher power still has an x in it, so it's zero at 0. So k!·c_k = f^(k)(0), and c_k = f^(k)(0)/k!. For sin, the derivatives at 0 cycle 0, 1, 0, −1: the even terms vanish and the odd terms alternate in sign. (Click) Back to the sin series: the first four terms.

Adding terms

Each term improves the fit

f(x) = \sin x

p_1(x) = x

p_3(x) = x - \frac{x^3}{3!}

p_5(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!}

p_7(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!}

Converges to \sin x
in the limit

From math to code

Taylor polynomial, degree 7

\sin x \approx x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!}

Factor out x

\sin x \approx \htmlClass{plain}{\htmlData{m=t1}{\htmlData{fly=f,drop=1}{x}} \htmlData{m=s3}{{}-{}} \htmlData{m=t3}{\frac{\htmlData{fly=f}{x}^{3}}{3!}} \htmlData{m=s5}{{}+{}} \htmlData{m=t5}{\frac{\htmlData{fly=f}{x}^{5}}{5!}} \htmlData{m=s7}{{}-{}} \htmlData{m=t7}{\frac{\htmlData{fly=f}{x}^{7}}{7!}}}

\sin x \approx \htmlClass{term}{\htmlData{m=t1}{\htmlData{fly=f,drop=1}{x}} \htmlData{m=s3}{{}-{}} \htmlData{m=t3}{\frac{\htmlData{fly=f}{x}^{3}}{3!}} \htmlData{m=s5}{{}+{}} \htmlData{m=t5}{\frac{\htmlData{fly=f}{x}^{5}}{5!}} \htmlData{m=s7}{{}-{}} \htmlData{m=t7}{\frac{\htmlData{fly=f}{x}^{7}}{7!}}}

\sin x \approx \htmlData{m=f}{x} \left(\htmlData{m=t1}{1} \htmlData{m=s3}{{}-{}} \htmlData{m=t3}{\frac{x^2}{3!}} \htmlData{m=s5}{{}+{}} \htmlData{m=t5}{\frac{x^4}{5!}} \htmlData{m=s7}{{}-{}} \htmlData{m=t7}{\frac{x^6}{7!}}\right)

Left with a polynomial in x^2

\sin x \approx x\,P(x^2)

P(u) = 1 - \frac{u}{3!} + \frac{u^2}{5!} - \frac{u^3}{7!}, \qquad u = x^2

From math to code, continued

Rewrite the polynomial for evaluation with Horner's method

P(u) = 1 - \frac{u}{3!} + \frac{u^2}{5!} - \frac{u^3}{7!}

Factor u out of the rest, again and again

P(u) = 1 \htmlClass{plain}{\htmlData{m=s3}{{}-{}} \htmlData{m=t3}{\frac{\htmlData{fly=f1,drop=1}{u}}{3!}} \htmlData{m=s5}{{}+{}} \htmlData{m=t5}{\frac{\htmlData{fly=f1}{u}^{2}}{5!}} \htmlData{m=s7}{{}-{}} \htmlData{m=t7}{\frac{\htmlData{fly=f1}{u}^{3}}{7!}}}

P(u) = 1 \htmlClass{term}{\htmlData{m=s3}{{}-{}} \htmlData{m=t3}{\frac{\htmlData{fly=f1,drop=1}{u}}{3!}} \htmlData{m=s5}{{}+{}} \htmlData{m=t5}{\frac{\htmlData{fly=f1}{u}^{2}}{5!}} \htmlData{m=s7}{{}-{}} \htmlData{m=t7}{\frac{\htmlData{fly=f1}{u}^{3}}{7!}}}

P(u) = 1 + \htmlData{m=f1}{u}\left(\htmlData{m=s3}{-}\htmlData{m=t3}{\frac{1}{3!}} \htmlClass{plain}{\htmlData{m=s5}{{}+{}} \htmlData{m=t5}{\frac{\htmlData{fly=f2,drop=1}{u}}{5!}} \htmlData{m=s7}{{}-{}} \htmlData{m=t7}{\frac{\htmlData{fly=f2}{u}^{2}}{7!}}}\right)

P(u) = 1 + \htmlData{m=f1}{u}\left(\htmlData{m=s3}{-}\htmlData{m=t3}{\frac{1}{3!}} \htmlClass{term}{\htmlData{m=s5}{{}+{}} \htmlData{m=t5}{\frac{\htmlData{fly=f2,drop=1}{u}}{5!}} \htmlData{m=s7}{{}-{}} \htmlData{m=t7}{\frac{\htmlData{fly=f2}{u}^{2}}{7!}}}\right)

P(u) = 1 + \htmlData{m=f1}{u}\left(\htmlData{m=s3}{-}\htmlData{m=t3}{\frac{1}{3!}} \htmlData{m=s5}{{}+{}} \htmlData{m=f2}{u}\left(\htmlData{m=t5}{\frac{1}{5!}} \htmlClass{plain}{\htmlData{m=s7}{{}-{}} \htmlData{m=t7}{\frac{\htmlData{fly=f3,drop=1}{u}}{7!}}}\right)\right)

P(u) = 1 + \htmlData{m=f1}{u}\left(\htmlData{m=s3}{-}\htmlData{m=t3}{\frac{1}{3!}} \htmlData{m=s5}{{}+{}} \htmlData{m=f2}{u}\left(\htmlData{m=t5}{\frac{1}{5!}} \htmlClass{term}{\htmlData{m=s7}{{}-{}} \htmlData{m=t7}{\frac{\htmlData{fly=f3,drop=1}{u}}{7!}}}\right)\right)

P(u) = \htmlClass{hl h4 c-orange}{1 + \htmlData{m=f1}{u}\left(\htmlClass{hl h3 c-orange}{\htmlData{m=s3}{-}\htmlData{m=t3}{\frac{1}{3!}} \htmlData{m=s5}{{}+{}} \htmlData{m=f2}{u}\left(\htmlClass{hl h2 c-orange}{\htmlData{m=t5}{\frac{1}{5!}} + \htmlData{m=f3}{u}\left(\htmlClass{hl h1 c-orange}{\htmlData{m=s7}{-}\htmlData{m=t7}{\frac{1}{7!}}}\right)}\right)}\right)}

Evaluate from the inside out

Written out, u³ is two multiplies, u² another, and so on. Instead: [click] look at P. [click] There's a u in these three terms, so factor it out: watch them slide to the front. [click] Now inside the bracket, there's a u in these two terms: [click] factor it out. [click] And a u in this one: [click] factor it out. That's Horner's method. [click] Now evaluate it, from the inside out: start with −1/7!. [click] Times u, plus 1/5!. [click] Times u, minus 1/3!. [click] Times u, plus 1. A number, times u, plus a number, over and over. [click] And that's the code: start with the top coefficient, and at each step multiply by u and add the next one. One variable, updated in place, and no powers at all. These are the actual double values, as the generated C++ will write them. What's nice: one multiply and one add per coefficient, and the compiler can fuse each pair into a single FMA (it does on arm64, and on x86 with -mfma). There are other schemes (Estrin's evaluates terms in parallel), but Horner is the one we'll use here.

Implementation

Putting it all together

Implementation, continued

Raising the degree

Measuring error

Absolute error
- |p(x) - f(x)|
- In dB: 20 \log_{10} |p(x) - f(x)|
Relative error
- |p(x) - f(x)| \,/\, |f(x)|
- Bits of precision: -\log_2(\text{relative error})
ULP distance
- |p(x) - f(x)| \,/\, \operatorname{ulp}(f(x))
- How many representable floats from the true value
Perceptual error

How did we do? Absolute error: how big the mistake is in the signal itself, in the output's own units. In dB it's the peak of the error you're adding to your waveform: −96 dB is about a 16-bit LSB. It's the worst case you'd see on a scope. (It doesn't care how big the answer is, which is also its blind spot near zero: keep that for the next slide.) Relative error: normalized, the error as a fraction of the answer, so it maps onto correct bits: one bit per halving, and a float holds 24, a relative error of about 6e−8. ULP distance: how far from the true answer, counted in representable floats. Within a binade it's about relative error × 2^23. Correctly rounded is at most half a ULP. Perceptual error: we may want an interpretation of the error in terms of what we hear. For sin, how loud are the harmonics? Their total energy (SINAD), or the amplitude of the largest one. The figure (top right): a toy float, 8 values per binade, so the gaps are visible. The true value f(x) = π/2 sits between floats; p(x) is a float, as every result is. Absolute: the gap, 0.30. Relative: zoom out to 0, the gap against |f(x)|, 0.19. ULP: zoom back in, the gap ruled off in the spacing at f(x), 2⁻³ here: about 2.4 floats away.

How did we do?

Taylor polynomial, degree 13, in float, on [-\pi, \pi]

Absolute error

Relative error

ULP distance

Harmonics

Binned max over samples

Log scale: the worst error at each x. [click] Absolute error: near 0 the polynomial is essentially exact, the floor is float rounding; the truncation error climbs as x^15/15! toward ±π, to about −93 dB. [click] Relative error: fine in the middle, but at f32(π), where sin is tiny, the polynomial is completely wrong: zero correct bits. [click] ULP distance: the same story counted in floats, 3 billion ULP at ±π (that one's over every float, from the C++ bench). [click] Harmonics: run it as an oscillator and the error repeats every cycle, so it shows up as harmonics, a broad spread around −120 dBc. The loudest single one (orange) is the 5th, at −117 dBc. SINAD adds them all up, every harmonic and any noise, against the fundamental: 105 dB, against 145 for float-rounded sin.

Scoreboard

Max over every normal float in [-\pi, \pi]

Apple M2 Pro, on power; clang 21, release, fast-math off. Time is ns per sample, one sample at a time over 65,536 sorted inputs, like an oscillator: the median of 100 passes, and of five launches (the next two slides show the program). 3× is the polynomial's own advantage. Let the compiler vectorize the loop and Taylor is 11×: std::sin is a library call and can't be. Sorted is std::sin's best case. Max ULP and abs error are exhaustive: every normal float in [−π, π], 2,140,282,808 inputs (about 2.1 billion, 1,070,141,404 each side; zero and subnormals left out). The check takes 11.4 s for all 9 sin methods on 12 threads (M2 Pro, on power): about 1.3 s per method. Loudest harmonic: one cycle in float, the single loudest harmonic in dBc. We'll add a row for every method.

Measuring speed

How fast is it? We can just call these functions, like std::sin, and time them. [click] Run a loop over 65,536 inputs, 100 times, [click] and take the median time per sample. [click] This is the loop we time. What's in it, and where we put escape, simulates different DSP scenarios. [click] escape: it forces the compiler to manifest everything we've written through this pointer, so it can't optimize our calculations away. If nothing reads the results, it deletes the loop (the first version of this measured 0.00 ns). escape is Chandler Carruth's, from his CppCon 2015 talk (about 45 minutes in): an empty asm statement that takes the pointer ("g": wherever it already is) and clobbers memory, so the compiler has to assume it reads and writes out. Every store has to happen, every pass. Google Benchmark's DoNotOptimize is the same idea. Standard C++ has no portable way to say this; on MSVC there's no inline asm, and it's a volatile write plus _ReadWriteBarrier instead. out has to come from outside: if the bench allocates it itself, the compiler can prove it doesn't overlap a lookup table, and in the vectorized loop (next slide) tables look 40% faster than they would in a plug-in. Why 65,536: std::sin branches on its input. Repeat 4096 random inputs and the branch predictor partly memorizes them; the time then varies from launch to launch (2.6 to 4.3 ns). The median, not the mean: it ignores the slow first pass and interrupts. 1000 runs gives the same numbers.

Three loops

Random inputs, sorted or unsorted

std::sin from Apple's libm

Estrin's scheme buys back in feedback

Three loops to put in the bench (bench/simple/loops.cpp). ns per sample, std::sin, Taylor 13, turns 13, over sorted inputs like an oscillator's sweep: std::sin's best case, so the speedups are conservative. Each loop: a click for the code, a click for its row, Taylor's speedup over std::sin, sorted and unsorted (median of five launches, slides/bench/simple-loops.txt). The table's empty to start; say it: Taylor's speedup over std::sin. Under it: MacBook Pro (M2 Pro) on power, Apple clang 21 -O3, no fast-math, std::sin from Apple's libm. [click] Buffer: the compiler may optimize within the buffer, and all we force is that the buffer is manifested. escape runs once, after the loop; inside it the compiler is free to vectorize, unroll and reorder across the whole buffer. clang inlines Taylor and swaps each scalar multiply and add for a NEON one, 4 lanes wide, and unrolls 4×: 16 samples per iteration, each one fmul for x², six fmla (Horner), one fmul by x. 1.8, 0.17, 0.23. [click] 10.3× sorted, 36.4× unsorted: std::sin's branches mispredict on random inputs, ours has none. About 11×, and 4× of that is the vector unit. std::sin is a scalar library call, so it doesn't vectorize. [click] Sample: the compiler may optimize within the sample, the way most per-sample DSP is written. escape after every store means each output has to be written, in order, before the next sample starts, so the compiler can only optimize the work inside one sample: no vectorizing. The CPU can still overlap neighboring samples, since each sample's work doesn't need the last one's: that's the out-of-order core, not the compiler. Taylor is the same 8 ops per sample, scalar; over 4 FP pipes that's about 2 cycles, 0.6 ns. std::sin gets the same overlap (1.85 ns here against 10.5 in feedback): it just has more instructions per sample, a call, branches, conversions to double. 1.85, 0.62, 0.81. [click] 3.0× sorted, 9.5× unsorted. This is the column every scoreboard quotes (sorted). About 3× and 2×: the polynomial's own advantage, and the number every scoreboard quotes. [click] Feedback: the previous output must be known; each input waits on it. 0 * y keeps the dependency, because we build without fast-math; fast-math would fold it away. Now it's latency, the longest dependent chain: about 10.5 ns for all three. Horner's 7 dependent multiply-adds cost what sinf's path does: no speedup, sadly. Feedback compressors, tanh inside a filter loop, FM feedback live here, and this is where a different evaluation scheme pays: Estrin's shortens the chain, 8 dependent ops to 5, 7.4 ns. Why is std::sin as fast as ours here? Apple's sinf (disassembled, notes/sinf-latency.md): convert to double, n = round(x · 2/π), r = x − n·π/2 in one fmsub, then a degree-7 sin or degree-8 cos polynomial on [−π/4, π/4], already Estrin-shaped, flip the sign by quadrant, convert back. About 10 dependent ops, against our 9: range reduction buys a low degree, Estrin a short chain. Its branches (size of x, sin or cos) are predicted, sorted inputs make that easy, so the core speculates past them and only the arithmetic chain counts. Random inputs (unsorted, same program): std::sin goes to 15.1 ns in feedback, Taylor and turns stay at 10.4 and 10.1: 1.45× for us, from the branches alone. (Sample loop: 5.8 against 0.61, 9.5×.) In short: buffer, approximations unlock vectorization; sample, they win on total work, the hardware hides the latency; feedback, only the dependent chain counts. [click] 1.0× sorted, 1.5× unsorted. The scoreboards quote the sample numbers. [click] And feedback is where we can buy some back: same coefficients in Estrin's scheme, a tree instead of a chain. 10.7 to 7.4 ns, 31% off. It costs about 20% in the other two loops, so Horner by default, Estrin on a feedback path (bench/simple/estrin.cpp, notes/latency-horner-estrin.md).

Padé approximants

Set the Taylor series equal to a rational function

{\textcolor{#2a78d6}{c_0 + c_1 x + c_2 x^2 + \cdots}} = \frac{\textcolor{#eb6834}{p_0 + p_1 x + \cdots + p_m x^m}}{\textcolor{#1baf7a}{1 + q_1 x + \cdots + q_n x^n}} + O\!\left(x^{m+n+1}\right)

Bring up the denominator

(\textcolor{#1baf7a}{1 + q_1 x + \cdots + q_n x^n})\,(\textcolor{#2a78d6}{c_0 + c_1 x + \cdots}) = \textcolor{#eb6834}{p_0 + p_1 x + \cdots + p_m x^m} + O\!\left(x^{m+n+1}\right)

Multiply out and match coefficients
- Linear system
- m + n + 1 equations and unknowns

\begin{alignedat}{2} \color{#8a8985}{x^0:}&\quad & p_0 &= c_0 \\ \color{#8a8985}{x^1:}&\quad & p_1 &= c_1 + q_1 c_0 \\ \color{#8a8985}{x^2:}&\quad & p_2 &= c_2 + q_1 c_1 + q_2 c_0 \\ & & &\;\;\vdots \\ \color{#8a8985}{x^{m+1}:}&\quad & 0 &= c_{m+1} + q_1 c_m + \cdots + q_n c_{m+1-n} \end{alignedat}

Start from the Taylor series and ask for a rational function, numerator degree m, denominator degree n, that agrees with it as far as possible. Fix q0 = 1 so the answer is unique. O(x^(m+n+1)) means "plus terms of degree m+n+1 and up": it's what turns "approximately equal" into "equal". The two sides agree through x^(m+n), and the O term is everything above that. Multiply both sides by the denominator: now both sides are polynomials, and the O term stays put: Q times O(x^(m+n+1)) is still O(x^(m+n+1)), since Q starts with 1 and can only push powers up. It's the mismatch between the two sides, Q·T − P, moved to the right: the left side is Q times the whole Taylor series, so it has every power of x, but P stops at x^m and we only have m+n+1 unknowns. So we ask for agreement through x^(m+n) only, and the O term is everything above that, where they're allowed to differ. It's not the Taylor series' own tail: that's already in the "⋯". Since Q(0) = 1, dividing through by Q keeps it: f − P/Q = O(x^(m+n+1)), so P/Q agrees with f through x^(m+n), like Taylor degree m+n. Each power of x gives one linear equation in the p's and q's: up to x^m each one gives a numerator coefficient; past x^m the left side is 0, and those n equations pin down the q's. m+n+1 equations, m+n+1 unknowns: a computer solves it, exactly, with fractions. We matched m+n+1 Taylor coefficients: the same information as a degree m+n Taylor polynomial, spent on a numerator and a denominator, and that generally behaves better toward the ends.

Padé vs Taylor

Padé in the wild

JUCE's FastMathApproximations header
- sin: Padé [7/6]
- tanh: Padé [7/6]

How did we do?

Padé approximant [7/6], in float, on [-\pi, \pi]

Absolute error

Relative error

ULP distance

Harmonics

Improving worst-case error

Problem: Taylor and Padé match f at one point
- The error gets large near the edges of the input interval
Goal: the error spread evenly across [a, b]
Idea: sample f at N points across the interval
- N points determine a polynomial of degree \le N - 1
- Fit a polynomial through them
  - Called an interpolant
  - Lagrange interpolation
- The error forms an oscillating curve through f

Equally spaced points

Equally spaced angles

Space the points evenly around a half circle, then project onto the x-axis

Move the points

Chebyshev nodes

Of the first kind

\begin{aligned} x_k &= \cos\theta_k \\[2pt] \theta_k &= \frac{(k + \frac12)\,\pi}{N} \\[7pt] k &= 0, \dots, N - 1 \end{aligned}

Chebyshev polynomials

T_n(\htmlClass{hl theta c-blue}{\cos\theta}) = \htmlClass{hl ntheta c-orange}{\cos(n\theta)}

x position on
unit circle
at \theta

x position on
unit circle
at n\theta

T_1(x) = x

T_2(x) = 2x^2 - 1

T_3(x) = 4x^3 - 3x

T_4(x) = 8x^4 - 8x^2 + 1

T_5(x) = 16x^5 - 20x^3 + 5x

One way to think about the Chebyshev polynomials: they're the function where, when you pass in an x position on the unit circle, [click] cos θ, the x position on the unit circle at θ, you get back [click] the x position on the unit circle at n times that angle, cos nθ. The top circle carries θ; its x drops straight down. The side circle carries nθ, turned a quarter so its x reads as height. Where they meet is (x, T_n(x)). (There's a T_0 too, the constant 1: the orange point never moves. I left it out because the animation doesn't look that good.) [click] n = 1: sweep θ from 0 to π, and the pen draws a line. [click] Double the angle: a parabola. [click] Triple it: a cubic, and the orange point laps its circle. [click] Four. [click] Five. The way this works out, cos nθ is always a polynomial in cos θ, is trig identities (cos 2θ = 2cos²θ − 1); no need to go further into the weeds. And every one wiggles between −1 and 1 with all its bumps the same height (equiripple: it comes back on the next slides). One of each degree, so they're a basis: any polynomial can be written in Chebyshev polynomials, and back again (NumPy's convert does it, in Implementation).

Adding Chebyshev terms

Each term improves the fit

f(x) = \sin x

Converges to \sin x
on [-\pi, \pi]

Like we did for sin with Taylor, let's build up an approximation, this time in Chebyshev polynomials. sin on [−π, π], so the T's are of x/π. [click] The line-ness: a_1 T_1. [click] Add the cubic-ness, a_3 T_3. [click] The quintic-ness. [click] And T_7: on top of sin across the whole interval, where Taylor was only good near 0. sin is odd, so the even terms are zero. (Secretly these are projections: the coefficients don't change as we add terms, so each partial sum is the best fit of its degree. Fitting each degree separately by interpolation would shift them a little.) Taylor at the same degree moves along in grey for comparison: great near 0, but peeling away toward ±π, where Chebyshev is still on sin. Keep adding terms and... [click] Pop-over: converges to sin x on [−π, π].

Chebyshev series

Any well-behaved f on [-1, 1]
can be written as
f(x) = \sum_{k=0}^{\infty} a_k\, T_k(x)
Find coefficients with the DCT
- Signal: the function to approximate
- Spectrum: the Chebyshev coefficients
Like audio, aliasing is an issue

\sin x on [-\pi, \pi]

Looks like a spectrum!

And as it turns out, in the same way the Taylor series converges to the true sin, the Chebyshev series converges as well. In fact any well-behaved function on the interval can be written as a sum of Chebyshev polynomials. (Precisely: Lipschitz continuous is enough for it to converge uniformly; Trefethen, ATAP, Theorem 3.1. The smoother it is, the faster the coefficients fall: like n^−k with k derivatives, geometrically for analytic functions like sin.) Below, the magnitudes of the coefficients, as a spectrum: the strength of each partial, each harmonic. sin is odd, so odd terms only, and they fall away fast. The shaded band is how many terms float needs, down to the dashed line at float resolution: about 13 for sin on [−π, π]. (log₂ on an octave needs about 8; that plot, cheb-coeffs-log2, is off the slide.) Why cut there: every T_k stays within ±1 with equal bumps, so cutting the series off leaves an error about the size of the first term dropped, spread evenly across the interval: the equiripple payoff. [click] Pop-over: looks like a spectrum! [click] And we find them the way we find a signal's harmonics. T_k of cos θ is cos kθ, so T_k really is the k-th harmonic, seen through x = cos θ. Sample and transform: the DFT wants samples equally spaced in time; the DCT wants them equally spaced in θ, which are exactly the Chebyshev nodes. (The DCT is the DFT of an even signal: only cosines.) NumPy's Chebyshev.interpolate does it in one call. [click] And like audio, aliasing comes along too. N samples can only see so many harmonics. [click] Say we interpolate at 10 points: that's our sample rate, and the dashed line is our Nyquist. Everything above it folds back, mirrored about N: a_11 lands on a_9, a_13 on a_7, and so on, and gets added into the coefficients we keep. [click] So do what audio does: oversample. Interpolate at a high degree, 128: Nyquist moves way out, where the coefficients are far below double precision, and the folded images go with it. [click] Then truncate at the same point, N = 10: keep a_1 … a_9, drop a_11 and up. Those are the clean coefficients of the degree-9 fit (for float we'd take N past 13, the band).

Implementation

Find the Chebyshev coefficients
Convert to powers of x

To finish the implementation, two steps. Find the Chebyshev coefficients: NumPy's Chebyshev.interpolate, with the domain doing the remap onto [−1, 1] for us, at a high degree, 128 (oversampled), and keep what we want: for degree 13, the first 14 coefficients. [click] Those are coefficients of Chebyshev polynomials; to evaluate with Horner we want powers of x. convert does both halves: the change of basis to powers, and the remap back from [−1, 1] onto [−π, π] (the series lives in t = x / π; t^k = x^k / π^k). Now it's the same shape as Taylor, ready for Horner, and these are the coefficients on the next slide. (Why not evaluate the Chebyshev form directly? You can, with Clenshaw's recurrence, but Horner is cheaper, and at degree 13 the conversion loses nothing.)

Chebyshev \sin x

Degree 13 on [-\pi, \pi], oversampled

How did we do?

Chebyshev, degree 13, in float, on [-\pi, \pi]

Absolute error

Relative error

ULP distance

Harmonics

Factor out the zeros

What the polynomial has to approximate

Zero at -\pi, 0 and \pi

\sin x \approx P(x)

Factor out the zero at 0

\sin x \approx x\,P(x^2)

Factor out the zeros at \pm\pi

\sin x \approx x\,\htmlClass{hl pifactors c-orange}{(x - \pi)(x + \pi)}\,P(x^2)

Example 04. For the next few slides we're going to invest some time, and a few operations, into this relative error problem, by looking at the structure of sin. The plot is what the polynomial has to approximate. [click] Straight at sin: it's zero at −π, 0 and π. So build those zeros into the form, and only approximate what's left. [click] Factor out the zero at 0: sin x over x, smooth and even (so a polynomial in x²), but still zero at ±π. [click] Factor those out too: what's left, sin x over x(x − π)(x + π), is a smooth, even bump between −1/π² and −1/2π² with no zeros (the axis zooms in ten times): easy to fit. Six coefficients in x². And that's the point for relative error: the polynomial never has to approximate a value near zero, so it never has to produce a tiny answer by cancelling big terms. The factors carry the zeros. One more multiply and one more add than before (the factors cost two of each; P is one term shorter). [click] But there's a problem, and it's these two factors. The π in them is the π we typed, and that's a float: f32(π) isn't π. So our zero isn't quite where sin's is. (Factoring alone: 1.8e7 → 1.2e7 ULP, barely moved.) Let's zoom in. Why factoring helps, the long version: notes/factoring-zeros.md.

Let's zoom in

The zero at \pi

\begin{aligned} \text{Math: } & (x - \pi) \\ \text{Code: } & (x - \mathrm{f32}(\pi)) \end{aligned}

Our curve crosses zero at the wrong place!

Splitting \pi

x - \pi \;\approx\; (x - \pi_{hi}) - \pi_{lo}

Now it crosses zero at \pi

Split \pi in C++

Factored zeros and split π: Colin, "Approximating sin(x) to 5 ULP with Chebyshev polynomials," mooooo.ooo (2017)

Example 05 as the bench builds it; factoring the zeros and splitting π is Colin's approach (mooooo.ooo, 2017: say it). The polynomial is the factored one, six coefficients in x². The new part: [click] π in two constants, π_hi and π_lo, [click] and the last three lines, the two factors built from them, (x − π_hi) − π_lo and (x + π_hi) + π_lo. Float vs double, since it's a template: - Every constant is written as a double literal and T(...) rounds it. For float, T(3.141592653589793) is f32(π) = 3.14159274, 8.7e−8 past π. For double it's f64(π), 1.2e−16 short of π. - π_lo is what's left over, so it depends on T: for float, π − f32(π) ≈ −8.74e−8 (computed in double, then rounded to float); for double, π − f64(π) ≈ 1.22e−16. Hence the is_same_v switch. - Double needs π_lo just as much, in relative terms: without it, a double version returns 0 at f64(π), 100% wrong, and is 22% off one float below. With it, relative error stays flat at about 6e−9 right up to π. In absolute terms the miss is only 1e−16, so if absolute error is all you care about, double π alone is fine. - The coefficients are fitted in double and rounded to float by T(...). That rounding costs a little; searching for the best float coefficients is example 08. - And degree 13 is sized for float. In double it's only good to about 1e−9: double would want a higher degree. It only works if the compiler does exactly these operations, in this order. -ffast-math lets the compiler reassociate: in real arithmetic (a − b) − c = a − (b + c), so it folds the two constants together, and π_hi + π_lo rounds straight back to π_hi. The split is gone. Measured: max ULP 4.44 → 1.2e7; the fast-math build is identical to the factored version on every float. A lot of audio code ships with fast-math on. [click] Hence the guard: three macros, wrapped around the function. The pragma gives the function precise semantics, so the two subtractions stay in order. On MSVC that's per function, and MSVC won't inline a precise function into a /fp:fast caller: it emits a call per sample, scalar, 4.8× slower. __forceinline makes it inline anyway, and the inlined body keeps its precise semantics: vectorized, bit-identical to /fp:precise. (Observed on MSVC 19.44, not documented by Microsoft.) Clang marks each instruction, so the pragma alone survives inlining; always_inline just matches. Guard where the function is defined: a pragma on the caller does nothing, on either compiler. On MSVC every inline wrapper between this and the hot loop needs __forceinline too. GCC is similar to MSVC: per-function options and inlining don't mix. Checked on every float: MSVC and clang-cl /fp:fast, 4.90 ULP (x64; my adc-msvc repo, github.com/rygrob/adc-msvc); Apple clang 21 -ffast-math on the M2, 4.44, identical to the precise build, at the precise build's speed. Unguarded: 1.2e7. (Off the slide: an #else branch leaves GCC unguarded; it has neither pragma.) The bench's plain cheb_split_sin_13 is this without the guard; the slide gives the guarded one its name.

How did we do?

Split \pi, degree 13, in float, on [-\pi, \pi]

Absolute error

Relative error

ULP distance

Harmonics

Scoreboard

Max over every normal float input

Effects of degree

Max absolute error

Max ULP

Speedup over std::sin

So what does a higher degree actually buy? Split π from degree 5 to 17; one term is two degrees. Every float; the dashed line is degree 13, where we are. [click] Absolute error: about 5 to 6 bits per term. The dashed line is the fit alone, in double, and it just keeps falling: sin has no singularities to slow it down. float stops at about 22 bits, at degree 11 to 13. Past that, you're paying for nothing. [click] ULP falls the same way, to 4.44 at degree 13, and more degree doesn't beat it: 4.79 at 15 and 17, just rounding. And the plain Chebyshev fit, zeros not built in, stops near 1e7: that's the form, and no degree fixes it. [click] And the price: the same for every term. Almost nothing in a buffer (0.02 ns), a little per sample (0.07), and a lot in feedback (1.1), where degree 13 is about where we stop beating std::sin (1.04× on sorted inputs). Sample: still 2×. So pick the form first, then raise the degree until float stops you. (notes/degree-effects.md: the same story for log₂ and tanh; a term costs the same, but log₂ gets 2.7 bits per term and tanh as a polynomial about 0.9, because of where their singularities are.)

Bonus: turns

\sin 2\pi x for x in [0, 1]

\sin 2\pi x \approx -u\,x\,(x - 1)\,P(u^2)

u = x - \tfrac{1}{2}

Zeros 0, ½, 1 are exact floats

Example 06. The problem was never the polynomial: it was π, which isn't a float. So change the input. Measure the angle in turns, so one period is [0, 1]. Now the zeros are at 0, ½ and 1, all exact floats, and the factors are built from x itself: x is exact, x − 1 is exact for x ≥ ½ by Sterbenz, u = x − ½ likewise. Nothing to split, nothing for fast-math to fold. And it's what an oscillator's phase accumulator holds anyway: phase in cycles, wrapping at 1. Underneath, its scoreboard row (max over every float in [0, 1]): the same accuracy as split π, quieter harmonics, slightly faster, and nothing the compiler can break. If your phase is in turns, this is the one to ship. But only if the phase stays in turns all the way in. Be aware of rounding wherever a transcendental constant like 2π sneaks back in. Two ways to throw it away (bench, max ULP over every input): - std::sin(2π·x): the multiply puts π back in, rounded, and std::sin is exact on the wrong input: 6.7e6 ULP next to x = ½. Our turns polynomial on the same x: 4.97. - A 32-bit integer phase accumulator, converted to float: near the top of the range the conversion rounds the phase up to a full turn. 1.6e9 for std::sin, 1.7e7 for our polynomial. Taking each factor straight from the integer (x, x − 1 and u, each exact) keeps 5.7.

Rounding out the toolkit

\tan x: poles
\log_2 x and 2^x: bit hacking
\tanh x: asymptotes

Building a \tan x approximation

\tan x

Zero at 0
Poles at \pm\pi/2

Dealing with the poles

\tan x \approx P(x)

Factor out the zero

\tan x \approx x\,P(x^2)

Divide out the poles

\tan x \approx \frac{x\,P(x^2)}{(\tfrac{\pi}{2} - x)(\tfrac{\pi}{2} + x)}

Put x in units of \pi

\tan \pi x \approx \frac{x\,Q(x^2)}{(\tfrac{1}{2} - x)(\tfrac{1}{2} + x)}

Implementation

\tan \pi x \approx \frac{x\,Q(x^2)}{(\tfrac{1}{2} - x)(\tfrac{1}{2} + x)}

Python
C++

The sin lesson again. Divide the poles out, fit the smooth quotient, and put the poles on exact floats. The Python: q is what we hand to Chebyshev.interpolate, the pole factor times tan(πx) over x. Plain tan(πx) / x is no good at 0: 0/0, and an interpolation node lands right there. So write tan as sine over cosine; numpy's sinc is sin(πx) / (πx), safe at 0, so tan(πx) / x = π sinc(x) / cos(πx). Then truncate, convert, and every other coefficient is Q's, in x². The radians version has to compute π·x first, and that rounding is as big as the whole distance to the pole for the last inputs. Splitting π/2 doesn't help: the error is already in t. The filter already has x = fc/fs, which is turns: Nyquist is ½, an exact float, and ½ − x is exact near it (Sterbenz). No π anywhere. Degree 9 is vsl::tan_pi, bit for bit. (Radians numbers sampled in Python; turns exhaustive.)

Scoreboard: \tan \pi x

Every normal float cutoff below Nyquist, at 48 kHz

Bonus: degree 5

Building \log_2 x and 2^x approximations

\log_2 x

Zero at 1
Challenges
- No symmetry, no periodicity
- Large domain, large range

Floating-point representation

Normal floats: scientific notation in base 2

x = \pm\,10^{e}\,m, \quad 1 \le m < 10 \qquad 440 = 10^{2} \times 4.4

x = \pm\,2^{e}\,m, \quad 1 \le m < 2 \qquad 440 = 2^{8} \times 1.71875

Binary representation

The log of a float

x = 2^{e}\,m, \quad 1 \le m < 2

\log_2 x = \log_2 \big(\htmlClass{term}{2^{e}\,m}\big)

\log_2 x = \htmlClass{term}{\log_2 2^{e}} + \log_2 m

\log_2 x = e + \log_2 m

e: an integer, from the exponent bits

\log_2 m: only need to approximate it on [1, 2)

Approximation strategy

Take e and m from the bits
\log_2 m: a polynomial on [1, 2)
Write the polynomial in m - 1
- \log_2 m = u\,P(u), u = m - 1
- P(u) = \dfrac{\log_2 (1 + u)}{u} = \dfrac{\ln (1 + u)}{u \ln 2}

log₂ goes on forever, but all we really need to approximate is this one region, 1 to 2. Every other octave is the same piece, stretched by 2^e and lifted by e. Pull the exponent field out as e, and reset the exponent to 0 to get m in [1, 2): bit manipulation, no arithmetic error at all (the Split the float slide, after the code). Then log₂ m on one octave. [click] And log₂ m is zero at m = 1, the same as sin's zeros: build it in. Write the polynomial in u = m − 1, which is exact (Sterbenz: m is within a factor 2 of 1), and factor out the u: log₂ m = u · P(u). What the polynomial actually fits is P, the orange curve: log₂ m / u, smooth, from 1/ln 2 ≈ 1.44 at m = 1 down to 1 at m = 2, nowhere near zero. A polynomial is good at that. The fit is the same call as sin's, with the quotient as a lambda: interpolate high, on [0, 1], truncate. (u = 0 is never evaluated: the Chebyshev nodes are inside the interval. log1p because 1 + u rounds near u = 0.) And near m = 1 the answer is u times something about 1.44, so it's as accurate relative to itself as P is: the zero is exact, and the values next to it keep their bits. Fit log₂ m directly instead and the polynomial has to land on zero by cancellation; its constant term's rounding alone is a huge relative error next to 1. (Maybe you've seen this: Mitchell, 1962, the original exploiter of this fact. Read all 32 bits as an integer, divide by 2^23, subtract 127: e + (m − 1), a straight line standing in for log₂ m, a chord across each octave. A rough log for one multiply and one add, about ±52 cents at worst with the offset tuned.)

\log_2 x in C++

Degree 7 on [1, 2)

Split the float

std::bit_cast and the IEEE layout

One function for float and double, with some type-trait shenanigans. std::bit_cast (C++20) reads the bits without undefined behavior, and compiles to nothing. The layout as traits, every member as wide as T: the unsigned and signed integers the bits fit in, how many significand bits are stored (one more is implicit, so float has 24 bits of precision and double 53), the exponent width, and the bias: 127 for float. split_exponent is then written once, and every type flows through at the float's own width. [click] bit_cast does the heavy lifting: the float's bits as an integer, free. Shift the exponent field down by the significand width (23), mask its 8 bits, subtract the bias: e. Keep the 23 fraction bits and set the exponent field to the bias, 2^0: m in [1, 2). The masks are built from the widths, (1 << 23) − 1: all ones over each field. Everything comes from three numbers: 23, 8, 127. pow2i, for 2^x, is the same traits run the other way. Positive normal x only: zero, subnormals, inf and NaN read the wrong e and m.

How did we do?

Degree 7 on [1, 2), in float; max over every positive float

Absolute error

Relative error

ULP distance

Move the octave

Any [c, 2c) will do

x = 2^e\,m, \quad m \in [1, 2)

m \ge 1.5:

\quad m \to m / 2

\quad e \to e + 1 \qquad 2^{e+1}\,\tfrac{m}{2} = 2^{e}\,m

m \in [0.75, 1.5)

Folding at 1.5: David Goldberg, "Fast Approximate Logarithms," eBay Tech Blog (2015)

We don't have to split at [1, 2). (Folding at 1.5 is straight from David Goldberg's eBay posts: say it.) The trouble is that log₂ m is zero at m = 1, right at the edge of the octave, so the floats just below 1 land at the far end, near 2, as −1 + 0.99999. Any octave from half a value to twice it works, so add a little logic to move it. [click] Take the top of the octave, m ≥ 1.5, [click] and halve it: it lands on [0.75, 1). On the curve it drops by exactly 1, since log₂(m/2) = log₂ m − 1. [click] So to keep x the same, bump the exponent: 2^(e+1) · m/2 = 2^e · m. One comparison and two selects. [click] Now the polynomial approximates log₂ m on [0.75, 1.5), through zero at 1 in the middle of it (u = m − 1 in [−¼, ½)): everything near 1 reads e = 0 and keeps its precision. (Folding at √2 centres it, slightly better: 20.7 ULP against 24.5. vsl::log2 folds at 1.5.)

\log_2 x in C++

Degree 7, folded at 1.5

How did we do?

Folded at 1.5, degree 7, in float; max over every positive float

Absolute error

Relative error

ULP distance

Scoreboard: \log_2 x

Every positive normal float

Bonus: 2^x in C++

Degree 5 on [-\tfrac{1}{2}, \tfrac{1}{2}]

As a bonus, the same technique works for 2^x, run backwards; we won't redo the algebra. Split x into an integer and a fraction: 2^x = 2^n · 2^f, n = round(x), f = x − n. The integer part, 2^n, is free (the exponent field, below). Only 2^f needs a polynomial, and it turns out we only have to approximate on [−½, ½]: every other unit interval is that piece times 2^n. f = x − round(x) is exact (Sterbenz: x and its rounding are within a factor 2), so f stays in [−½, ½] with no error. Round rather than floor: x − floor(x) isn't exact for negative x. Centring doesn't help the fit: 2^f on any unit interval is the same curve, scaled. (Schraudolph, 1999, is the 2^x version of Mitchell: write x · 2^23 into the bits directly, a chord per octave, ±51 cents.) The code: round, subtract, then 2^f = 1 + f·Q(f), with Q a Chebyshev fit again, oversampled, of (2^f − 1)/f: total degree 5. The 1 is built in, the way log₂'s zero is, so every integer x gives exactly 2^n. Multiply by 2^n. pow2i on the right, on the same ieee<T> traits as split_exponent: add the bias, shift it into the exponent field. nearbyint rounds in the current mode (round-to-nearest-even by default), without the halfway-away of std::round, and vectorizes: in a loop, clang emits frinti.4s, then fcvtzs.4s and shl.4s for pow2i, four lanes at once. Edges: below −126 or above 127 the exponent field wraps. Clamp if you need the full range; for mtof and dB it doesn't come up. Underneath, its scoreboard row: 3.2 ULP over every normal float in [−125, 126], 1.9× std::exp2 one sample at a time.

Building a \tanh x approximation

\tanh x

Zero at 0
Asymptotes at \pm 1
- Exactly 1 in float past

Polynomials and rationals run off

Polynomial, degree 9

Fits only on \pm

Padé [7/6]

Passes 1 at
Off by at 9

Further out

It keeps climbing

Top: the function. Bottom: the error, log scale, under the dashed line 1 − tanh x, what clamping at x would cost. [click] A polynomial can't level off. Fit degree 9 on ±2.80 and it's fine inside, then runs off the moment it leaves: to use it you clamp at 2.80 and pay 1 − tanh(2.80), −42 dB. [click] Padé [7/6], JUCE's. It looks right on top. Underneath, it's superb near 0 and climbs from about 3: it passes 1 at 4.97, and by 9 it's 5e−3 off. [click] And stretch the axis out to 30: it just keeps going, 1.37 at 30, growing like x. An odd over an even rational can't level off either; the degrees differ by one. The best you can do is clamp around 5, about −80 dB. Neither reaches 9: we have to clamp, and to clamp at 9 we need something that holds up all the way there.

Build in the asymptote

\htmlClass{hl t1 c-blue}{\tanh x = \frac{\sinh x}{\sqrt{1 + \sinh^2 x}}}

\htmlClass{hl t2 c-blue}{y = \frac{s}{\sqrt{1 + s^2}}, \quad s = x\,R(x^2)}

\htmlClass{hl t3 c-blue}{R(x^2) = 1 + x^2\,Q(x^2)}

An algebraic sigmoid, bounded by construction. The same lesson as sin and tan: put the structure you know into the form. [click] tanh is sinh over cosh, and cosh is √(1 + sinh²): exact. [click] So approximate sinh with an odd polynomial s, x times R of x², and the formula does the saturating: whatever s is, s/√(1 + s²) stays inside ±1, smooth and monotonic (the pop-over: s = x itself, the plain algebraic sigmoid, beside tanh: the right shape, too slow to saturate). One square root. This is vsl::tanh's form (and Raph Levien's, 2018). [click] s stands in for sinh x, so R(x²) stands in for sinh x / x, whose Taylor series is 1 + x²/3! + x⁴/5! + …: even, so a polynomial in x², and it starts with exactly 1. So keep the 1 in the form and fit a polynomial to the rest, Q = (sinh x / x − 1) / x² (the fit, not the Taylor terms: Q's coefficients come out a little off 1/3!, 1/5!, … to spread the error). [click] In the pop-over, the blue curve moves onto tanh: s is now x R(x²), our fit. The exact 1 means y = x for tiny x: 24 ULP, where a fitted constant gives 82. [click] The code, the Python first. The fit: Q is the same Chebyshev fit as everywhere else, with one twist: only on x in [0, 3.4] (v is x², so 3.4²). Past that tanh is within 0.2% of 1, and an error in s barely moves y (dy/ds falls like sech³), so there's no point spending accuracy there; fit sinh all the way to 9 and its growth swamps the fit (−15 dB). Past 3.4 the polynomial keeps growing like sinh, so y still saturates. (3.4 was scanned; the optimum is flat.) −117 dB over every float, matching vsl's hand-tuned coefficients. [click] The C++: the clamp at ±9, two instructions (fmin, fmax): without it, s² overflows once |x| passes 556 and y collapses to 0, then NaN. With it, every float is safe, and the clamp costs nothing since tanh is 1 there anyway. Then Horner in x², the 1 last, and one square root. (If asked: std::sqrt is one correctly rounded instruction, and clang vectorizes it, fsqrt.4s. The rsqrt estimate plus Newton steps needs two steps on NEON to match, and is slower than sqrt and divide on the M2.) (Further, if asked: Remez, minimax weighted by the error in y, is about 4 dB better at the same cost, −121 dB. NumPy and SciPy have no general Remez (scipy.signal.remez is filter design): hand roll it (ours is about 30 lines, examples/common/remez.py), or as a linear program, CVXPY on a dense grid matches our Remez to about six digits. Off the shelf elsewhere: Sollya, Chebfun, baryrat.)

Comparison

Polynomial, clamp 2.80

Padé [7/6], clamp

Sigmoid, clamp 9

The dashed line is 1 − tanh x: what clamping at x costs. A method can be clamped where its error meets that line, if it's still accurate there. [click] The polynomial: clamped at 2.80, −42 dB. [click] Padé: excellent near 0, but it climbs through 1 just before 5, so clamp its input where its error meets the line, 4.79: −83 dB over every float. (JUCE clamps the output instead; past about 1e7, x⁶ overflows in both P and Q, inf / inf is NaN, and an output clamp passes NaN through. Clamping the input avoids it, and fmin / fmax is faster than std::clamp: one instruction per bound, not compare and select.) [click] The sigmoid stays around 1e−6 all the way out, so it reaches 9, where the clamp costs below float resolution: −117 dB over every float, the clamp included. As the fit gets better, the clamp moves out.

Scoreboard: \tanh x

Every normal float

Every normal float, so a blown-up input counts: each method clamps its input (the Padé at 4.79, −83 dB; the poly + clamp row is off the slide, still in the data). Clamped at 9, the sigmoid is −117 dB everywhere, no blow-up at the ends (the asymptote is built in), at 1.8× std::tanh one sample at a time, 6× vectorized (std::tanh branches on its input: 2.0 ns sorted, 3.4 on random inputs; the polynomials don't care). If you need more, raise the degree: Chebyshev degree 11 (fit on [0, 3.9]) is about −133 dB and 4 ULP for a little more time, or Remez, −121 dB at degree 9. Max ULP: 24, from R(0) = 1; for a saturator −117 dB is the number. (Off the slide: the unclamped sigmoid returns 0 past 556, then NaN; an exp2-based expm1 form reaches −136 dB at 1.2× (5× vectorized).)

Recap

\sin x
- Taylor, Padé, Chebyshev
- Error metrics and benchmarks
- Floating-point considerations
\tan x
- Divide out the poles
\log_2 x and 2^x
- Bit hacking, in float and double
\tanh x
- Sigmoid structure

Lessons

The toolkit
- Padé approximants
- Chebyshev polynomials
- Remez algorithm
Know your function
- Zeros, poles, asymptotes
Know your requirements
- What is your domain?
- How much accuracy? Does ULP matter?
Know your floats

Measure
- Worst-case error
- Realistic DSP scenarios

What's the toolkit? Padé, Chebyshev polynomials. And to take it a little further, the Remez algorithm: the true minimax fit, where Chebyshev gets close (tanh: −121 dB against −117 at the same degree). Know your function: its zeros (sin), poles (tan), asymptotes (tanh). Build them in, and fit whatever's left. Think about your requirements: how much accuracy do you need? Do you care about ULP distance, or is −120 dB enough? Sometimes that extra ULP improvement costs operations, and for a saturator or a meter you don't need it. Think about floating point representation: which steps are exact (Sterbenz), which constants aren't floats (π), and what the bits give you for free (the exponent). Think about the domain: reduce it (turns, one octave, [−½, ½]) and clamp it where the function is done (tanh at 9). And measure: the worst-case error, and a float has only 4 billion values, so check every one; and the speed in realistic DSP scenarios (buffer, sample, feedback loops), not just one.

Acknowledgements

Lloyd N. Trefethen, Approximation Theory and Approximation Practice
Jason Sachs, "Chebyshev Approximation…"
- embeddedrelated.com/showarticle/152.php
Colin, "Approximating sin(x) to 5 ULP…"
- mooooo.ooo/chebyshev-sine-approximation
David Goldberg, "Fast Approximate Logarithms," parts I–III (eBay)
Raph Levien, "A Few of My Favorite Sigmoids"
- raphlinus.github.io/audio/2018/09/05/sigmoid.html
Jatin Chowdhury, math_approx
- github.com/Chowdhury-DSP/math_approx

Thank you

functions.rygrob.com
github.com/rygrob/fun-with-functions