Part D · Numerics › Nonlinearities without dividers

Module 9·Part D — Numerics·15 min

Nonlinearities without dividers

LUT, piecewise-linear, minimax polynomial and CORDIC — and why range reduction is the step that makes any of them work.

The core mental model

Division and transcendental functions are the expensive primitives in hardware. A divider is iterative, so its latency scales with the word width and it occupies real area; exp and tanh have no direct circuit at all. Everything in this module is a way to replace them with the operations that are cheap — add, shift, multiply, and table lookup — and the four standard techniques sit at different points on the same area-versus-latency curve. A lookup table is O(1)O(1) latency with area exponential in input width, viable to about ten or twelve bits. Piecewise-linear stores a slope and intercept per segment and costs one multiply-add, with area linear in segment count. A minimax polynomial costs dd multiply-adds for degree dd, with coefficients chosen to minimise worst-case rather than average error.

less area, more latencylookup tableone access — all area, no arithmeticpiecewise linearsmall table + one multiply-addpolynomial (minimax)a few multipliers, no tableCORDICshifts and adds only — one bit per iteration
Four ways to compute the same transcendental, ordered by what they spend. A lookup table is one memory access and all area; a polynomial trades area for multipliers; CORDIC needs neither but takes an iteration per bit. The right answer is set by which resource is scarce on your part — and on an FPGA with block RAM to spare, the unglamorous table usually wins.

The step that makes all four viable is range reduction, and it is the part most people skip. You never approximate a function over its whole domain. For exe^x, write x=kln2+rx = k\ln 2 + r with r12ln2\lvert r \rvert \le \tfrac{1}{2}\ln 2; then ex=2kere^x = 2^k \cdot e^r, where 2k2^k is an exponent adjustment — free in floating point, a shift in fixed point — and ere^r needs only a low-degree polynomial on an interval of width 0.69. A degree-3 minimax polynomial reaches about 10510^{-5} relative error there, where approximating exe^x directly over [20,20][-20, 20] would need an absurd degree and would suffer catastrophic cancellation at the ends. Range reduction converts an impossible approximation problem into an easy one, every time.

Softmax is the worked example that ties it together. Subtracting the maximum is presented as numerical stability, and in floating point that is what it is; in fixed point it is also a range reduction, because after subtraction every argument is 0\le 0, so ex(0,1]e^x \in (0, 1] and you can spend every bit of an unsigned fractional format on precision instead of reserving integer bits for a dynamic range you no longer have. The division by the sum is then the only division in the whole operation, it happens once per vector rather than per element, and you do not need a divider for it either: seed a reciprocal from a small table and run one or two Newton-Raphson steps, yn+1=yn(2xyn)y_{n+1} = y_n(2 - x y_n), which doubles the number of correct bits each iteration using two multiplies and a subtract.

Numbers worth memorizing

TechniqueLatencyAreaSweet spot
LUT1 access2n×m2^n \times m bits≤ 10–12 input bits
Piecewise linear1 MAC + lookupsegments × (slope, icpt)ML activations, 1e-3 to 1e-4
Minimax polynomialdd MACs (Horner)dd multipliersgeneral, after range reduction
CORDIC~WW iterationsshifters and adders onlyno multiplier available
QuantityValue
LUT size, 8-bit in / 16-bit out256 × 16 = 4 Kb — trivial
LUT size, 16-bit in / 16-bit out64 K × 16 = 1 Mb — a real budget
Minimax degree 3 for ere^r on [ln2/2,ln2/2][-\ln2/2, \ln2/2]~1e-5 relative
Minimax degree 5, same interval~1e-8, roughly fp32
CORDIC convergence1 bit per iteration
CORDIC gain KK≈ 1.64676 — pre-scale by 1/K1/K
Newton-Raphsondoubles correct bits per iteration
8-bit seed → after 1 / 2 iterations16 bits / 32 bits
tanh\tanh saturationx>4tanhx>0.9993\lvert x \rvert > 4 \Rightarrow \lvert \tanh x \rvert > 0.9993
Sigmoid saturationx>8σ>0.99966x > 8 \Rightarrow \sigma > 0.99966
Iterative divider latency~WW cycles
LUT + 1 Newton reciprocal~3–4 cycles

Identities that halve the work:

tanh(x)=2σ(2x)1,σ(x)=1σ(x),ex=2ker,  x=kln2+r\tanh(x) = 2\sigma(2x) - 1, \qquad \sigma(-x) = 1 - \sigma(x), \qquad e^x = 2^{k} e^{r},\; x = k\ln 2 + r

The first means you build one circuit, not two. The second means you tabulate half the domain. The third is range reduction itself.

Critical thinking

Implement softmax over 64 elements in fixed point with no divider. Give the latency in levels.

Six stages, and the interesting result is where the depth actually goes.

  1. Max — a comparator tree over 64 inputs: log264=\lceil \log_2 64 \rceil = 6 levels.
  2. Subtract — one level. Now every argument is 0\le 0, so the exponential lands in (0,1](0,1] and an unsigned fractional format spends all its bits on precision.
  3. Exponential — range-reduce x=kln2+rx = k\ln 2 + r (one multiply-add and a rounding to get kk), apply a degree-3 minimax polynomial in rr by Horner, then scale by 2k2^k with a shift: about 4–5 levels.
  4. Sum — an adder tree: 6 levels, and widen the accumulator by 6 bits per Module 8.
  5. Reciprocal — an 8-bit LUT seed plus one Newton step y(2xy)y(2 - xy): about 4 levels.
  6. Multiply — one level.

Total ≈ 22 levels.

The point worth extracting: the two tree reductions contribute 12 of those 22 levels, while the transcendental everyone worries about contributes 4–5. The reductions dominate, not the nonlinearity. If this were on a deadline, the productive move would be Module 3’s — attack the adder and comparator trees, not the exponential. It is a good reminder that intuitions about “the expensive operation” are usually about area or about software cost, and depth is a different accounting entirely.

When is a lookup table the wrong answer?

Four conditions, and the first is simply arithmetic.

  • Input width. Area is 2n×m2^n \times m bits. Eight bits in is 4 Kb and free; sixteen bits in is 1 Mb, which is a serious fraction of an FPGA’s block RAM for one function. Past about twelve bits, tables stop being the cheap option.
  • Wide dynamic range. A table spends entries uniformly, but functions like exe^x need resolution concentrated where they vary. Most of a naive table is wasted on flat regions — precisely the problem range reduction plus a small polynomial solves with a fraction of the area.
  • Many instances. If the design needs 64 parallel activation units, the table area multiplies by 64 while a piecewise-linear unit’s coefficient storage can often be shared.
  • Memory latency and porting. Block RAM is a 1–2 cycle read with a limited number of ports. Sixty-four parallel lookups means replicating the table or serialising the accesses, and serialising is exactly what you cannot afford.

The pattern to notice: LUTs win when the input is narrow and the instance count is low. As either grows, the answer shifts toward piecewise-linear, and then toward a polynomial with range reduction.

Is subtracting the max in softmax just numerical stability?

In floating point, yes — it prevents exe^x overflowing for large xx.

In fixed point it does something more valuable. Without it, the arguments span whatever range your logits happen to have, so exe^x spans an enormous dynamic range and your Q-format must reserve integer bits for the largest possible value. Those integer bits are taken from the fraction, so every value in the vector loses resolution because of one large element you were going to divide out anyway.

After subtracting the max, every argument is 0\le 0 and therefore ex(0,1]e^x \in (0, 1]. You can now choose an unsigned Q0.nn format — no integer bits at all — and spend the entire word on precision. On a 16-bit datapath that is the difference between roughly 10 fractional bits and 16, which is a factor of 64 in resolution.

So the operation converts a dynamic range problem into a static one, which is the recurring move of this whole module: range reduction lets you fix a format at design time instead of provisioning for a worst case that almost never occurs. Same instruction, different reason, considerably larger payoff.

CORDIC or polynomial, for a 20 ns budget at 250 MHz?

Twenty nanoseconds at 250 MHz is 5 cycles, and that settles it.

CORDIC produces one bit per iteration, so 16-bit output needs 16 iterations. Pipelined that is 16 cycles — over budget by 3×. Fully unrolled into combinational logic it is a chain of 16 dependent shift-add stages, which will not close timing anywhere near 250 MHz. You would also need to pre-scale by 1/K0.60731/K \approx 0.6073 to compensate the CORDIC gain.

Degree-3 minimax by Horner is three dependent multiply-adds. On DSP blocks that is roughly 3–4 cycles including the range reduction — inside budget, with room.

And if it were tight, Estrin’s scheme restructures the polynomial so the dependent depth is about log2d\log_2 d rather than dd, at the cost of more multipliers running in parallel. For degree 7 that is 3 levels instead of 7, bought with extra DSPs. Which is Module 4’s trade — area for depth — appearing inside a polynomial evaluation.

So: CORDIC is for area-constrained, latency-tolerant designs, or platforms with no multiplier at all. With DSP blocks available and a deadline, polynomials win, and Estrin is the lever when Horner is not fast enough.

How accurate does the nonlinearity actually need to be?

Ask what the consumer can resolve, not what you can achieve. This is where a lot of area gets wasted.

If the activation feeds an int8 quantised layer, the output is going to be rounded to 1/256 — about 3.9×1033.9\times10^{-3}. An approximation error below a quarter LSB, roughly 10310^{-3}, is invisible: it is smaller than the quantisation you are about to apply anyway. A piecewise-linear sigmoid with 32 segments comfortably achieves that, and a degree-5 polynomial delivering 10810^{-8} is spending multipliers to produce bits that get discarded one stage later.

The right procedure is to propagate the error budget backwards from the output tolerance, through each stage’s quantisation, and size every approximation to sit just under the noise floor its consumer imposes. Combined with Module 8’s accounting — where the accumulator width and the requantisation shift are already determined — this usually shows that the nonlinearity is nowhere near the accuracy bottleneck.

The failure mode is real and common: engineers optimise the approximation because it is the part with a clean mathematical objective, while the actual error is dominated by an 8-bit requantisation two stages downstream. Over-engineering the transcendental is one of the most frequent ways to spend area for nothing.

Self-check

Write the range reduction for exp and say why it matters.

Write x=kln2+rx = k\ln 2 + r with r12ln2\lvert r \rvert \le \tfrac{1}{2}\ln 2, so ex=2kere^x = 2^k e^r. The 2k2^k factor is an exponent adjustment — free in float, a shift in fixed point — and ere^r only needs approximating on an interval of width 0.69, where a degree-3 minimax polynomial reaches about 10510^{-5}. Without it the required degree over a realistic domain explodes and the ends suffer catastrophic cancellation.

Compute a reciprocal with no divider, and state the convergence rate.

Seed from a small lookup table, then iterate yn+1=yn(2xyn)y_{n+1} = y_n(2 - x y_n) — two multiplies and a subtract per step. Newton-Raphson doubles the number of correct bits each iteration, so an 8-bit seed gives 16 bits after one step and 32 after two. Total latency is roughly 3–4 cycles against ~WW cycles for an iterative divider.

In the 64-element fixed-point softmax, where does the latency actually go?

The two tree reductions dominate: a 6-level comparator tree for the max and a 6-level adder tree for the sum, 12 of about 22 total levels. The exponential — range reduction plus a degree-3 polynomial — is only 4–5, and the reciprocal about 4. If the design were latency-critical you would attack the trees (Module 3), not the transcendental, which is the opposite of most people’s first instinct.

Give the non-stability reason to subtract the max in fixed-point softmax.

It is a range reduction. Afterwards every argument is 0\le 0, so ex(0,1]e^x \in (0,1] and the format can be unsigned Q0.nn with no integer bits, putting the entire word into precision. On a 16-bit datapath that is roughly 10 fractional bits becoming 16 — a factor of 64 in resolution. It converts a dynamic-range problem into a static one fixed at design time.

Why does CORDIC lose to a polynomial under a tight deadline, and what is Estrin's scheme for?

CORDIC yields one bit per iteration, so 16 bits costs 16 dependent shift-add stages — too deep to pipeline within a few cycles and too long to close timing combinationally. A degree-3 Horner polynomial is three dependent multiply-adds, about 3–4 cycles with DSP blocks. Estrin’s scheme restructures a polynomial so dependent depth is roughly log2d\log_2 d instead of dd, buying depth with extra parallel multipliers — Module 4’s area-for-latency trade inside the arithmetic.

How do you decide the required accuracy of an activation approximation?

By the error budget of its consumer. Feeding an int8 layer means quantisation to about 3.9×1033.9\times10^{-3}, so approximation error below roughly 10310^{-3} is invisible and a 32-segment piecewise-linear fit suffices. Propagate the output tolerance backwards through each stage’s quantisation and size each approximation just under the noise floor it faces — over-engineering the transcendental while an 8-bit requantisation downstream dominates the error is a common and expensive mistake.