Solidity exp(): A Fixed-Point Exponential in 289 Gas
The exponential function is everywhere in DeFi: compounding, discounting, interest rates, and the pricing of options, binaries and futures.
In my DeFiMath library, exp() sits underneath most of that — it is a building block for European and binary options, for futures, and for other math functions like the standard normal cumulative distribution. Anything embedded that deeply has to be both cheap and accurate, because every function above it inherits whatever it costs and whatever it gets wrong.
Here's how I wrote it. It computes e^x in 289 gas, with a max relative error of 2.2e-14 — and 3.0e-16 absolute, the metric that matters when exp(x) < 1. We'll look at the code first, then work through how it gets there.
Exp Function Code
/// @notice Computes the natural exponential of x in 18-decimal fixed-point.
/// @dev Reverts with ExpUpperBoundError when x >= EXP_UPPER_BOUND (135e18);
/// returns 0 when x <= EXP_LOWER_BOUND (~-41.446e18).
/// Max relative error: < 2.2e-14 for any y >= 1e18.
/// Max absolute error: < 3.0e-16 for any y < 1e18.
/// @param x Signed input in 18-decimal fixed-point format.
/// @return y Result e^x in 18-decimal fixed-point format.
function exp(int256 x) internal pure returns (uint256 y) {
unchecked {
if (x >= 0) {
// check input
if (x >= EXP_UPPER_BOUND) revert ExpUpperBoundError();
uint256 x_ = uint256(x);
// The algorithm works in 3 steps:
// 1) reduce the range of X to [0, ln(2)/64]
// 2) approximate result in a narrow range
// 3) recover reduction
//
// Two-stage range reduction:
// stage 1 — split x = k·ln(2) + r with integer k, r ∈ [0, ln(2)).
// Then exp(x) = exp(k·ln(2)) · exp(r) = 2^k · exp(r).
// stage 2 — divide r by 64 (right-shift by 6): r' = r/64 ∈ [0, ~0.0108].
uint256 k = x_ / LN_2;
x_ -= k * LN_2;
x_ >>= 6;
// Next, we use Padé[3/3] approximant using the following formula:
// exp(x) ≈ (120 + 60x + 12x² + x³) / (120 - 60x + 12x² - x³)
uint256 x2 = x_ * x_;
uint256 even = 120e54 + x2 * 12e18; // 120 + 12x² (even powers of x)
uint256 odd = x_ * (60e36 + x2); // 60x + x³ (odd powers of x)
uint256 p = (even + odd) * 1e18; // numerator: 120 + 60x + 12x² + x³
uint256 q = even - odd; // denominator: 120 − 60x + 12x² − x³
assembly ("memory-safe") {
y := div(p, q)
}
// Finally, undo both reductions in reverse order:
// stage 2 — raise y to the 64th power to recover exp(r) from exp(r/64).
// stage 1 — left-shift by k to multiply by 2^k and undo the ln(2) factoring.
y = y * y;
y = y * y / 1e54;
y = y * y;
y = y * y / 1e54;
y = y * y;
y = y * y / 1e54;
y <<= k;
} else {
// check input
if (x <= EXP_LOWER_BOUND) return 0;
// Negative branch: exp(x) = 1 / exp(|x|). Runs the same algorithm as the
// positive branch on x_ = |x|
uint256 x_ = uint256(-x);
// range reduction (see positive branch)
uint256 k = x_ / LN_2;
x_ -= k * LN_2;
x_ >>= 6;
// Padé[3/3] approximant (see positive branch)
uint256 x2 = x_ * x_;
uint256 even = 120e54 + x2 * 12e18; // 120 + 12x² (even powers of x)
uint256 odd = x_ * (60e36 + x2); // 60x + x³ (odd powers of x)
uint256 p = (even + odd) * 1e18; // numerator: 120 + 60x + 12x² + x³
uint256 q = even - odd; // denominator: 120 − 60x + 12x² − x³
assembly ("memory-safe") {
y := div(p, q)
}
// undo reductions (see positive branch)
y = y * y;
y = y * y / 1e54;
y = y * y;
y = y * y / 1e54;
y = y * y;
y = y * y / 1e54;
y <<= k;
// reciprocate: exp(-|x|) = 1 / exp(|x|)
assembly ("memory-safe") {
y := div(1000000000000000000000000000000000000, y)
}
}
}
}Like the sqrt function, exp is split into two branches, one for positive and one for negative input x. The code is almost identical in both, and the identity exp(-x) = 1 / exp(x) folds the negative case back into the positive one — the only extra work is a single division at the very end.
Let's take a closer look at the positive branch.
3-step Algorithm
Computing an exponential is straightforward enough on paper — the Wikipedia article on the exponential function lists several ways to do it. But we are doing it in Solidity, where every iteration costs gas. A Taylor series gives us all the precision we could want and charges far too much for it. What we need is an approximation that stays simple and still lands inside our error budget.
Range reduction
Using the following identity, we can reduce the range in which we approximate, from the full input domain down to a single narrow interval:
exp(x) = exp(k·ln(2) + r) = 2^k · exp(r)
Every x splits into an integer multiple of ln(2) plus a remainder: k is the integer division x / ln(2), and r is what's left over, so r ∈ [0, ln(2)). That means we only ever have to approximate exp() on an interval of width 0.693, and the 2^k factor is recovered for free at the very end with a single left-shift.
That stage alone isn't enough though — 0.693 is still too wide for a cheap approximant to hit our precision target. So we reduce once more, this time using the fact that dividing the argument turns into taking a root of the result:
exp(r) = exp(r / 64)^64
Dividing by 64 is a right-shift by 6, and raising back to the 64th power is just 6 successive squarings — no pow, no loop, no lookup table. In exchange, the interval shrinks by a factor of 64, from 0.693 down to ln(2) / 64 ≈ 0.0108. That narrow interval is what lets a very cheap approximation stay accurate enough — and that approximation is what we'll look at next.
Put together, the whole two-stage reduction is three lines of code:
uint256 k = x_ / LN_2;
x_ -= k * LN_2;
x_ >>= 6;Both stages get undone at the very end, and that recovery turns out to be just as cheap — we'll come back to it.
Approximation
I have tested multiple approximants in the narrow range of [0, 0.0108] to find the most efficient one. I wanted to minimize gas cost while maintaining high precision. In the end, I settled for the Padé[3/3] approximant:
exp(x) ≈ (120 + 60x + 12x² + x³) / (120 - 60x + 12x² - x³)
On [0, 0.0108] this approximant is far more accurate than it needs to be: its own error peaks at 1.7e-19 relative, right at the top of the interval. Recovery then raises the result to the 64th power, and raising to a power multiplies relative error by that power — so about 1.1e-17 of it survives into the output.
That is three orders of magnitude below the 2.2e-14 the function guarantees end to end, which is the whole point of reducing the range first. The approximation is not what limits precision here. Fixed-point truncation is — chiefly the >>= 6, which discards up to 63 wei of the argument before the approximant ever sees it, and those 63 wei get multiplied by 64 on the way back out.
Recovery
Recovery runs in the reverse order of the reductions: first we raise the result to the 64th power, then we multiply it by 2^k. The 64th power costs six squarings, since 2⁶ = 64, and they come in pairs — a plain y * y, then a second one that also divides by 1e54 to get back to 1e18 base. Because y stays very close to 1e18 throughout, the intermediates never overflow. The final 2^k is a left-shift:
y = y * y;
y = y * y / 1e54;
y = y * y;
y = y * y / 1e54;
y = y * y;
y = y * y / 1e54;
y <<= k;Monotonicity
Does exp() ever step down? For a primitive sitting under option premiums and discount factors, a dip would be bad news — a premium falling while its input rises, a binary search that stops converging. So it is worth being precise about how much we can actually say.
The answer comes in two parts, and they carry different weight.
The approximation is provably increasing. Write it as N(x)/N(-x), with N(x) = x³ + 12x² + 60x + 120. Differentiate, and the mess cancels:
R'(x) = 24·(x⁴ - 30x² + 600) / N(-x)²
The bottom is a square. The top has no real roots — substitute u = x², and the discriminant is 900 - 2400 = -1500 — so with a positive leading coefficient it is positive everywhere. R'(x) > 0, always. And nothing downstream can undo that: truncating division, the squarings and the final shift all preserve order.
The reduction seams are checked, not proved. Each time k ticks up, r drops from just under ln(2) back to 0 while the result doubles. On paper the two sides match exactly; in integers they differ by a wei or two, and the algebra says nothing about which way. There are 255 such points, so test_MONO_expAtReductionSeams sweeps every one of them, wei by wei. All pass.
Together that gives non-decreasing, not strictly increasing — truncation makes neighbouring inputs share a wei, and exp(0) and exp(-1 wei) both return exactly 1e18.
Wrapping up
The interesting part of this function turned out not to be the approximation. The Padé[3/3] is accurate to 1.7e-19 on the reduced interval, and even after recovery multiplies that by 64 it still lands three orders of magnitude inside the error budget. What the gas actually buys is the range reduction: three lines going in, seven coming out. Those ten lines are what let a small rational function stand in for e^x across the entire domain — shrink the interval far enough, and the approximation problem mostly disappears.
Try it
npm install defimath-libThis exp is one of dozens of primitives in DeFiMath, a gas-optimized fixed-point math library for Solidity — expm1, ln, pow, sqrt and a full Black-Scholes options suite are all built the same way. The full reference for this function lives on the exp docs page. It's MIT-licensed, pure Solidity, with zero runtime dependencies. Source is on GitHub, with reference benchmarks against PRBMath, ABDK and Solady in defimath-compare.
If you want more of the same, I wrote up sqrt in 197 gas and the CLZ opcode that several of these primitives lean on. The approximation used here is a Padé approximant of order [3/3].