# WAD math for the solved-AMM strategy

## Format and scope

Use the five `WadMath` entrypoints below: internal pure Solidity 0.8.24 functions, signed `int256` input and output, scale **10^18**. For example, `ln(2e18)` returns approximately `693147180559945309`; `cdf(0)` returns `500000000000000000`. Inputs are real numbers multiplied by WAD, not unscaled integers. Solidity signed division truncates toward zero. Exp's binary downscaling and positive probability divisions truncate downward.

This is the bounded **math-wad** support item. It supplies math, measurements and integration guidance, not a strategy or compressed policy. Taylor (gtmcknight) is credited in the library header. None of the goal's replica rewards, real-engine scores or P/Q/F policy gates is claimed here.

| Routine | Domain safe from overflow/revert, in real units | Deliberate endpoint behavior |
|---|---|---|
| ln | Every positive representable input: raw integer 1 through `type(int256).max`, i.e. 10^-18 through `(2^255-1)/10^18` | Nonpositive input reverts `LnDomain()` |
| exp | Every representable input <= 130 | Input <= -42 returns zero; input > 130 reverts `ExpDomain()` |
| erf | Every `int256` | Saturates to ±1 for |x| >= 6; exact zero at zero |
| pdf | Every `int256` | Returns zero for |x| >= 10; exp/WAD quantization also produces zeros earlier |
| cdf | Every `int256` | Saturates to 0/1 for x <= -9 / x >= 9; erf saturation can also saturate it earlier |

Overflow safety is broader than the measured accuracy domains. `norm`, `reduce` and `scale` are implementation helpers; use them only as the supplied entrypoints do. In particular, `scale` is not a general overflow-checked multiplication API.

## Code: complete selected library

All five routines are in this one standalone library; no import, constructor, linked deployment, storage or dependency is required. The remaining Solidity block contains only comparison/measurement contracts and is not needed by a strategy.

```solidity
// SPDX-License-Identifier: MIT
pragma solidity 0.8.24;
// Math support for Taylor (gtmcknight)'s solved AMM policy; no policy code is ported here.
// Locally derived ln/exp and piecewise erfc interpolation. See result.md for provenance.
library WadMath {
    int256 internal constant W = 1e18;
    int256 internal constant L = 693147180559945309;
    int256 internal constant IS2 = 707106781186547524;
    int256 internal constant IS2PI = 398942280401432678;
    error LnDomain();
    error ExpDomain();
    function msb(uint256 u) private pure returns (uint256 r) { unchecked {
        if (u >= 1 << 128) {u >>= 128; r += 128;}
        if (u >= 1 << 64) {u >>= 64; r += 64;}
        if (u >= 1 << 32) {u >>= 32; r += 32;}
        if (u >= 1 << 16) {u >>= 16; r += 16;}
        if (u >= 1 << 8) {u >>= 8; r += 8;}
        if (u >= 1 << 4) {u >>= 4; r += 4;}
        if (u >= 1 << 2) {u >>= 2; r += 2;}
        if (u >= 1 << 1) r += 1;
    } }
    function norm(int256 x) internal pure returns (int256 m, int256 k) { unchecked {
        if (x <= 0) revert LnDomain();
        k = int256(msb(uint256(x))) - 59;
        m = k >= 0 ? int256(uint256(x) >> uint256(k)) :
                     int256(uint256(x) << uint256(-k));
        if (m < W) {m *= 2; --k;}
    } }
    function ln(int256 x) internal pure returns (int256) { unchecked {
        (int256 m, int256 k) = norm(x);
        if (m > 1414213562373095048) {m /= 2; ++k;}
        int256 z = (m-W)*W/(m+W);
        int256 q = z*z/W;
        int256 p;
        p = 86956521739130434;
        p = 95238095238095238 + p*q/W;
        p = 105263157894736842 + p*q/W;
        p = 117647058823529411 + p*q/W;
        p = 133333333333333333 + p*q/W;
        p = 153846153846153846 + p*q/W;
        p = 181818181818181818 + p*q/W;
        p = 222222222222222222 + p*q/W;
        p = 285714285714285714 + p*q/W;
        p = 400000000000000000 + p*q/W;
        p = 666666666666666666 + p*q/W;
        p = 2000000000000000000 + p*q/W;
        return k*L + z*p/W;
    } }
    function reduce(int256 x) internal pure returns (int256 r, int256 k) { unchecked {
        k = x / L;
        r = x - k*L;
        if (r > L/2) {r -= L; ++k;}
        if (r < -L/2) {r += L; --k;}
    } }
    function scale(int256 y, int256 k) internal pure returns (int256) { unchecked {
        return k >= 0 ? int256(uint256(y) << uint256(k)) :
                        int256(uint256(y) >> uint256(-k));
    } }
    function exp(int256 x) internal pure returns (int256) { unchecked {
        if (x <= -42*W) return 0;
        if (x > 130*W) revert ExpDomain();
        (int256 r, int256 k) = reduce(x);
        int256 p; int256 q;
        p = 1503126503126;
        q = 1503126503126;
        p = 63131313131313 + p*r/W;
        q = 63131313131313 - q*r/W;
        p = 1262626262626262 + p*r/W;
        q = 1262626262626262 - q*r/W;
        p = 15151515151515151 + p*r/W;
        q = 15151515151515151 - q*r/W;
        p = 113636363636363636 + p*r/W;
        q = 113636363636363636 - q*r/W;
        p = 500000000000000000 + p*r/W;
        q = 500000000000000000 - q*r/W;
        p = 1000000000000000000 + p*r/W;
        q = 1000000000000000000 - q*r/W;
        return scale(p*W/q, k);
    } }
    function erf(int256 x) internal pure returns (int256) { unchecked {
        if (x >= 6*W) return W;
        if (x <= -6*W) return -W;
        if (x == 0) return 0;
        bool neg = x < 0;
        int256 a = neg ? -x : x;
        uint256 j = uint256(a/W);
        int256 t = 2*(a-int256(j)*W)-W;
        int256 p;
        if (j == 0) {
            p = 322729;
            p = 9470031 + p*t/W;
            p = -23918910 + p*t/W;
            p = -359971579 + p*t/W;
            p = 1208750157 + p*t/W;
            p = 10403079024 + p*t/W;
            p = -47938662076 + p*t/W;
            p = -242844471010 + p*t/W;
            p = 1532711689771 + p*t/W;
            p = 4233225590655 + p*t/W;
            p = -39248392632403 + p*t/W;
            p = -42228205371650 + p*t/W;
            p = 781902902264501 + p*t/W;
            p = -228849629949995 + p*t/W;
            p = -11442481496555265 + p*t/W;
            p = 18307970394489505 + p*t/W;
            p = 109847822366930599 + p*t/W;
            p = -439391289467722415 + p*t/W;
            p = 479500122186953462 + p*t/W;
        }
        else if (j == 1) {
            p = 385504;
            p = 2141606 + p*t/W;
            p = -23358504 + p*t/W;
            p = -17344978 + p*t/W;
            p = 860103862 + p*t/W;
            p = -2192700871 + p*t/W;
            p = -18967271976 + p*t/W;
            p = 124242956993 + p*t/W;
            p = 90736602605 + p*t/W;
            p = -3339491537607 + p*t/W;
            p = 9229209457676 + p*t/W;
            p = 37055103186218 + p*t/W;
            p = -301971437482081 + p*t/W;
            p = 464571442273314 + p*t/W;
            p = 2787428653678825 + p*t/W;
            p = -17344000511778915 + p*t/W;
            p = 44598858458861014 + p*t/W;
            p = -59465144611814692 + p*t/W;
            p = 33894853524689273 + p*t/W;
        }
        else if (j == 2) {
            p = 85292;
            p = -173751 + p*t/W;
            p = -2724122 + p*t/W;
            p = 19639129 + p*t/W;
            p = -11081619 + p*t/W;
            p = -572699810 + p*t/W;
            p = 3437511025 + p*t/W;
            p = -2508254661 + p*t/W;
            p = -76954425765 + p*t/W;
            p = 488838626262 + p*t/W;
            p = -1018242866084 + p*t/W;
            p = -4237578599787 + p*t/W;
            p = 43726408703825 + p*t/W;
            p = -191167131674724 + p*t/W;
            p = 538898442404972 + p*t/W;
            p = -1043761193710643 + p*t/W;
            p = 1361427643970443 + p*t/W;
            p = -1089142115176355 + p*t/W;
            p = 406952017444959 + p*t/W;
        }
        else if (j == 3) {
            p = 2169;
            p = -29346 + p*t/W;
            p = 127117 + p*t/W;
            p = 171926 + p*t/W;
            p = -5635607 + p*t/W;
            p = 36828229 + p*t/W;
            p = -108342526 + p*t/W;
            p = -217142657 + p*t/W;
            p = 4532243569 + p*t/W;
            p = -29942905911 + p*t/W;
            p = 133824904060 + p*t/W;
            p = -454630864410 + p*t/W;
            p = 1214753849843 + p*t/W;
            p = -2566133819956 + p*t/W;
            p = 4232363177065 + p*t/W;
            p = -5286938719524 + p*t/W;
            p = 4724498430212 + p*t/W;
            p = -2699713388692 + p*t/W;
            p = 743098372341 + p*t/W;
        }
        else if (j == 4) {
            p = -100;
            p = 449 + p*t/W;
            p = 5 + p*t/W;
            p = -14196 + p*t/W;
            p = 135717 + p*t/W;
            p = -860590 + p*t/W;
            p = 4275123 + p*t/W;
            p = -17560797 + p*t/W;
            p = 60987893 + p*t/W;
            p = -180617068 + p*t/W;
            p = 456520535 + p*t/W;
            p = -980087992 + p*t/W;
            p = 1769207185 + p*t/W;
            p = -2641959458 + p*t/W;
            p = 3183936145 + p*t/W;
            p = -2981107620 + p*t/W;
            p = 2037719133 + p*t/W;
            p = -905652948 + p*t/W;
            p = 196616044 + p*t/W;
        }
        else if (j == 5) {
            p = 1;
            p = -5 + p*t/W;
            p = 23 + p*t/W;
            p = -97 + p*t/W;
            p = 378 + p*t/W;
            p = -1297 + p*t/W;
            p = 3963 + p*t/W;
            p = -10812 + p*t/W;
            p = 26227 + p*t/W;
            p = -56231 + p*t/W;
            p = 105742 + p*t/W;
            p = -172593 + p*t/W;
            p = 241080 + p*t/W;
            p = -282698 + p*t/W;
            p = 270898 + p*t/W;
            p = -203870 + p*t/W;
            p = 113071 + p*t/W;
            p = -41117 + p*t/W;
            p = 7358 + p*t/W;
        }
        // p approximates erfc(a). Clamp polynomial rounding at the boundaries.
        if (p < 0) p = 0;
        if (p > W) p = W;
        int256 y = W-p;
        return neg ? -y : y;
    } }
    function pdf(int256 x) internal pure returns (int256) { unchecked {
        if (x >= 10*W || x <= -10*W) return 0;
        return exp(-x*x/(2*W))*IS2PI/W;
    } }
    function cdf(int256 x) internal pure returns (int256) { unchecked {
        if (x >= 9*W) return W;
        if (x <= -9*W) return 0;
        return (W+erf(x*IS2/W))/2;
    } }
}
```

Bounds behind `unchecked`: normalization shifts a positive input to a mantissa in [1,2) WAD, with exponent between -60 and 195. Log's second adjustment puts the mantissa approximately in [1/sqrt(2),sqrt(2)], so |z| < 0.172 WAD. Log intermediates are below 10^38 in raw units. For exp, the explicit guards give |r| <= ln(2)/2 and k between -61 and 188. Numerator and denominator Horner intermediates are below 2 WAD, the denominator stays above 0.8 WAD, and `p*W` is below 2×10^36. The final shifted result at x=130 is approximately 2.873×10^74 raw, below 2^255. Erf branches before negating, including at `int256.min`; |t| <= WAD and each polynomial's coefficient absolute sum is below 2 WAD, keeping Horner products below 2×10^36. Pdf branches before squaring, bounding the raw square by 10^38. Cdf branches before multiplication, bounding its input to erf by 9/sqrt(2). These bounds also justify the unchecked comparison implementations at their guarded entrypoints.

## Measurement conventions

All following numbers are **local measurements from the verbatim harness below**, not figures taken from a paper. Setup succeeded at amm-challenge commit `6b6921f8d463314aa7c96234e1be25a25245f704`, including wheel checksum and `amm-match validate contracts/src/StarterStrategy.sol`. Toolchain: CPython 3.12.3, py-solc-x 2.1.0, solc 0.8.24, pyrevm 0.3.7. Compilation uses optimizer enabled, 200 runs, Cancun, no via-IR; execution uses pyrevm `spec_id='CANCUN'`.

Each method uses 20,001 points: ln logarithmically spaced on [10^-6,10^6], exp uniformly on [-40,40], erf uniformly on [-6,6], pdf and cdf uniformly on [-8,8]. Each grid is quantized to raw WAD integers before evaluation. Every one of the **200,010** selected/alternative grid evaluations also ran in pyrevm and matched the integer Python model **exactly**. Additional checks passed: 90 successful domain-boundary cases, 228 join/saturation cases, six expected reverts and ten bundled-call cases. Gas samples are ten additional calls per method.

References are Python `math.log`, `math.exp`, `math.erf`, `exp(-x*x/2)/sqrt(2*pi)`, and `erfc(-x/sqrt(2))/2`, respectively, evaluated at the quantized input converted to float. Errors below are in real units; relative error divides by the reference magnitude and excludes a reference exactly zero. Binary64 reference/conversion rounding is included: a measured zero error is not proof of exact arithmetic. Tail relative errors deserve attention even when absolute error is small. The tables give max/mean over the whole stated grid, with an additional central-region relative table later.

Gas is `evm.result.gas_used`. **Execution gas** subtracts exactly `21000 + 4*zero_calldata_bytes + 16*nonzero_calldata_bytes`; it still includes the wrapper's ABI decoding, method dispatch and return. **Transaction gas** is the unadjusted pyrevm result. Calls have a 1,000,000 gas cap, no storage or external calls, and no gas refunds. These are isolated function costs, not `afterSwap` engine statistics. Each gas table uses the ten inputs printed directly below it; negative calldata has more nonzero bytes.

## ln: choice, error and gas

Selected: binary normalization followed by the atanh series
`ln(m) = 2*z*(1 + z²/3 + ... + z²²/23)`, `z=(m-1)/(m+1)`, evaluated in Horner form. With |z| < 0.172, the real-arithmetic omitted tail is below 7×10^-21; fixed-point rounding and the ln(2) constant dominate earlier. The coefficients are integer truncations of `2*W/(2i+1)`.

Alternative: binary normalization plus **64 repeated squarings** to extract fractional log2 bits, then multiply by ln(2). It has similar measured accuracy but much higher gas. Select the short atanh polynomial.

Code: `WadMath.ln` in the complete library above; comparison code: `Alternatives.ln` below.

| Method | Max absolute | Mean absolute | Max relative | Mean relative |
|---|---:|---:|---:|---:|
| Selected | 1.77635684e-15 | 4.01984285e-17 | 3.39021375e-14 | 2.46350141e-17 |
| Alternative | 1.77635684e-15 | 4.06233929e-17 | 3.31173658e-14 | 2.47781721e-17 |

| Method | Execution min / mean / max | Transaction min / mean / max |
|---|---:|---:|
| Selected | 1246 / 1273.4 / 1314 | 22620 / 22664.2 / 22706 |
| Alternative | 6674 / 7691.6 / 8292 | 28078 / 29094.4 / 29708 |

Gas inputs (real units, in call order): `1e-06, 0.0001, 0.01, 0.5, 1.0, 2.0, 10.0, 100.0, 1000.0, 1000000.0`.

Selected grid's observed peak execution gas: **1320**; alternative: **8835**. This grid peak is distinct from the ten-input gas table and is not a proof of a global maximum.

## exp: choice, error and gas

Selected: nearest-ln(2) range reduction, a **[6/6] Padé rational** on [-ln(2)/2,ln(2)/2], then binary scaling. The numerator coefficients are `c_k=(12-k)!*6!/(12!*k!*(6-k)!)`, quantized downward to WAD; the denominator uses alternating signs.

Alternative: identical range reduction, a **16-term Taylor recurrence** plus constant term, then binary scaling. Both have essentially the same measured error, but Padé uses substantially less gas. Select Padé. The large maximum absolute error occurs near the top of this exponential range and must be read with the relative errors; it is not an 800-unit error near exp(0). At the low end, quantization of values around a few raw WAD units dominates relative error.

Code: `WadMath.exp` in the complete library above; comparison code: `Alternatives.exp` below.

| Method | Max absolute | Mean absolute | Max relative | Mean relative |
|---|---:|---:|---:|---:|
| Selected | 8.00000000e+02 | 5.22014483e+00 | 1.97671627e-01 | 1.47681015e-03 |
| Alternative | 8.00000000e+02 | 5.22491259e+00 | 1.97671627e-01 | 1.47681015e-03 |

| Method | Execution min / mean / max | Transaction min / mean / max |
|---|---:|---:|
| Selected | 1183 / 1200.0 / 1226 | 22515 / 22714.4 / 22918 |
| Alternative | 2356 / 2373.0 / 2399 | 23688 / 23887.4 / 24091 |

Gas inputs (real units, in call order): `-40.0, -20.0, -5.0, -1.0, 0.0, 1.0, 5.0, 10.0, 20.0, 40.0`.

Selected grid's observed peak execution gas: **1226**; alternative: **2399**. This grid peak is distinct from the ten-input gas table and is not a proof of a global maximum.

## erf: choice, error and gas

Selected: exploit odd symmetry and approximate **erfc** on six unit intervals [0,1], ..., [5,6]. Each interval uses a degree-18 polynomial in `t=2*(|x|-j)-1`. Nineteen first-kind Chebyshev nodes are interpolated with 95-digit Decimal arithmetic, converted to power coefficients, rounded to nearest WAD, then evaluated with Horner. The harness regenerates and asserts every shipped coefficient before measuring. Approximating erfc directly preserves small tails better than fitting values very close to one; final output is clamped to [-1,1].

Alternative: **Abramowitz–Stegun 7.1.26**, a five-coefficient polynomial in `1/(1+p*|x|)` multiplied by exp(-x²), using the coefficients and formula from [John D. Cook's public-domain implementation](https://www.johndcook.com/blog/cpp_erf/). It uses the selected Padé exp to isolate the erf approximation's cost. The measured error is materially larger and it is slightly slower here; select piecewise interpolation, accepting more bytecode. The underlying erf integral definition is given by [NIST DLMF 7.2.1](https://dlmf.nist.gov/7.2.E1.tex).

Code: `WadMath.erf` in the complete library above; comparison code: `Alternatives.erf` below.

| Method | Max absolute | Mean absolute | Max relative | Mean relative |
|---|---:|---:|---:|---:|
| Selected | 1.11022302e-16 | 4.73888665e-18 | 2.72240627e-15 | 9.15868464e-18 |
| Alternative | 1.39374503e-07 | 4.27681472e-08 | 7.78722220e-06 | 1.05460373e-07 |

| Method | Execution min / mean / max | Transaction min / mean / max |
|---|---:|---:|
| Selected | 471 / 1311.3 / 1710 | 21850 / 22823.3 / 23392 |
| Alternative | 494 / 1402.7 / 1806 | 21873 / 22914.7 / 23488 |

Gas inputs (real units, in call order): `-6.0, -4.0, -2.0, -1.0, 0.0, 0.25, 1.0, 2.0, 4.0, 6.0`.

Selected grid's observed peak execution gas: **1714**; alternative: **1806**. This grid peak is distinct from the ten-input gas table and is not a proof of a global maximum.

## pdf: choice, error and gas

Selected: `exp(-x*x/2)/sqrt(2*pi)`, using Padé exp and a WAD constant for `1/sqrt(2*pi)`.

Alternative: the same density identity using the Taylor exp implementation. This compares two ways of computing the exponential density, with identical squaring and final scaling. Select Padé: no meaningful accuracy advantage from Taylor, and higher gas. In the tails the final WAD multiplication/division truncation dominates relative error.

Code: `WadMath.pdf` in the complete library above; comparison code: `Alternatives.pdf` below.

| Method | Max absolute | Mean absolute | Max relative | Mean relative |
|---|---:|---:|---:|---:|
| Selected | 1.11022302e-16 | 7.92141912e-18 | 2.32433625e-04 | 2.29079950e-06 |
| Alternative | 1.11022302e-16 | 7.97359916e-18 | 2.32433625e-04 | 2.29079950e-06 |

| Method | Execution min / mean / max | Transaction min / mean / max |
|---|---:|---:|
| Selected | 1446 / 1470.4 / 1489 | 22778 / 22982.4 / 23181 |
| Alternative | 2607 / 2631.4 / 2650 | 23939 / 24143.4 / 24342 |

Gas inputs (real units, in call order): `-8.0, -6.0, -3.0, -1.0, 0.0, 0.5, 1.0, 3.0, 6.0, 8.0`.

Selected grid's observed peak execution gas: **1489**; alternative: **2650**. This grid peak is distinct from the ten-input gas table and is not a proof of a global maximum.

## cdf: choice, error and gas

Selected: `(1+erf(x/sqrt(2)))/2`, using piecewise erf. Alternative: the same identity using the A&S rational-in-t/exp erf implementation. These compare the two different probability approximations, not just a changed coefficient count. Select piecewise erf. Python's reference instead uses erfc to avoid binary64 subtraction cancellation in the negative tail. The WAD result remains quantized; a direct higher-precision log-CDF or wider internal scale would be needed for accurate relative probabilities arbitrarily deep in that tail.

Code: `WadMath.cdf` in the complete library above; comparison code: `Alternatives.cdf` below.

| Method | Max absolute | Mean absolute | Max relative | Mean relative |
|---|---:|---:|---:|---:|
| Selected | 1.11022302e-16 | 5.46155064e-18 | 7.55602298e-04 | 3.78047578e-06 |
| Alternative | 6.96874125e-08 | 2.26812438e-08 | 1.00230732e-02 | 1.05332975e-03 |

| Method | Execution min / mean / max | Transaction min / mean / max |
|---|---:|---:|
| Selected | 745 / 1779.0 / 1941 | 22077 / 23291.0 / 23623 |
| Alternative | 768 / 1888.2 / 2033 | 22100 / 23400.2 / 23715 |

Gas inputs (real units, in call order): `-8.0, -6.0, -3.0, -1.0, 0.0, 0.5, 1.0, 3.0, 6.0, 8.0`.

Selected grid's observed peak execution gas: **1941**; alternative: **2033**. This grid peak is distinct from the ten-input gas table and is not a proof of a global maximum.

## Relative error away from the most quantized regions

The following uses subsets of the same grids. The full-range relative errors above remain the appropriate warning for tails. Ln excludes the immediate neighborhood of its zero; erf excludes very small inputs; exp/pdf/cdf subsets reduce deep-tail quantization.

| Routine and region | Selected max / mean relative | Alternative max / mean relative |
|---|---:|---:|
| ln: |ln(x)| >= 0.001 | 3.39021375e-14 / 2.46350141e-17 | 3.31173658e-14 / 2.47781721e-17 |
| exp: x >= -20 | 4.23663826e-10 / 3.88726070e-12 | 4.23663826e-10 / 3.88726094e-12 |
| erf: |x| >= 0.001 | 1.12099122e-15 / 8.88733275e-18 | 6.98647176e-06 / 1.04692120e-07 |
| pdf: |x| <= 6 | 1.90144707e-10 / 3.39328097e-12 | 1.90144707e-10 / 3.39328173e-12 |
| cdf: x >= -6 | 1.27022085e-09 / 1.37881616e-11 | 3.58230714e-03 / 2.63471269e-04 |

## Recommendation

Use all five selected `WadMath` routines. The selected implementation trades several kilobytes of code for cheap Gaussian evaluation, which is a useful trade when the callback needs several Gaussian updates. It embeds 114 erfc coefficients, not the goal's approximately 71,000 policy values. The policy table still needs its own compression.

A separate `SelectedProbe.bundle` calls **each routine five times**, varying its input on every iteration. Its ten base inputs are 0.1, 0.2, 0.5, 1, 2, 3, 4, 5, 6 and 8; within each group the actual input is `base + i/10`, for i=0..4. Every transaction succeeded with a 245,000 gas cap:

| Five calls to each selected routine | Min | Mean | Max |
|---|---:|---:|---:|
| Execution gas | 26230 | 30673.2 | 32401 |
| Transaction gas | 47494 | 51937.2 | 53665 |

The combined selected+alternative `Probe` runtime is **6617 bytes**. The separate selected-only `SelectedProbe`, including both its dispatcher and five-of-each bundle function, is **6002 bytes**. These are measured deployed bytecode sizes, not Solidity source sizes. Both are below 24,576 bytes, but neither includes a strategy or policy table. There is ample measured math-only gas headroom relative to 245,000; the assembled strategy must still establish its own size and worst callback gas on the real engine.

## Comparison contracts and thin wrappers

Append this block to the library when reproducing the measurements. The `run` dispatch IDs are 0/1 ln selected/alternative, 2/3 exp, 4/5 erf, 6/7 pdf, 8/9 cdf. ID 10 returns its input as an optional ABI baseline; no baseline subtraction was used. Helpers are shared where the two approaches intentionally use the same normalization or range reduction.

```solidity
library Alternatives {
    int256 constant W=1e18;
    int256 constant L=693147180559945309;
    function ln(int256 x) internal pure returns (int256) { unchecked {
        (int256 y, int256 k)=WadMath.norm(x);
        uint256 f;
        for (uint256 i=0; i<64; ++i) {
            y=y*y/W;
            if (y>=2*W) {y/=2; f |= uint256(1) << (63-i);}
        }
        return k*L+int256(f)*L/int256(uint256(1)<<64);
    } }
    function exp(int256 x) internal pure returns (int256) { unchecked {
        if (x<=-42*W) return 0;
        if (x>130*W) revert WadMath.ExpDomain();
        (int256 r,int256 k)=WadMath.reduce(x);
        int256 term=W; int256 y=W;
        for (int256 i=1;i<=16;++i) {term=term*r/(W*i); y+=term;}
        return WadMath.scale(y,k);
    } }
    // A&S 7.1.26, coefficients from John D. Cook's public-domain implementation.
    function erf(int256 x) internal pure returns (int256) { unchecked {
        if (x>=6*W) return W;
        if (x<=-6*W) return -W;
        if (x==0) return 0;
        bool neg=x<0; int256 a=neg ? -x : x;
        int256 t=W*W/(W+327591100000000000*a/W);
        int256 p=1061405429000000000;
        p=p*t/W-1453152027000000000;
        p=p*t/W+1421413741000000000;
        p=p*t/W-284496736000000000;
        p=p*t/W+254829592000000000;
        p=p*t/W;
        int256 y=W-p*WadMath.exp(-a*a/W)/W;
        return neg ? -y : y;
    } }
    function pdf(int256 x) internal pure returns (int256) { unchecked {
        if (x>=10*W || x<=-10*W) return 0;
        return exp(-x*x/(2*W))*398942280401432678/W;
    } }
    function cdf(int256 x) internal pure returns (int256) { unchecked {
        if (x>=9*W) return W;
        if (x<=-9*W) return 0;
        return (W+erf(x*707106781186547524/W))/2;
    } }
}
contract Probe {
    function run(uint256 method,int256 x) external pure returns (int256) {
        if(method==0) return WadMath.ln(x);
        if(method==1) return Alternatives.ln(x);
        if(method==2) return WadMath.exp(x);
        if(method==3) return Alternatives.exp(x);
        if(method==4) return WadMath.erf(x);
        if(method==5) return Alternatives.erf(x);
        if(method==6) return WadMath.pdf(x);
        if(method==7) return Alternatives.pdf(x);
        if(method==8) return WadMath.cdf(x);
        if(method==9) return Alternatives.cdf(x);
        return x; // method 10: ABI/dispatch baseline, not subtracted from results
    }
}
contract SelectedProbe {
 function run(uint256 m,int256 x) external pure returns(int256) {
  if(m==0) return WadMath.ln(x);
  if(m==1) return WadMath.exp(x);
  if(m==2) return WadMath.erf(x);
  if(m==3) return WadMath.pdf(x);
  return WadMath.cdf(x);
 }
 function bundle(int256 x) external pure returns(int256 s) {
  for(uint256 i=0;i<5;++i) {
   int256 t=x+int256(i)*1e17;
   s+=WadMath.ln(t);s+=WadMath.exp(t);s+=WadMath.erf(t);
   s+=WadMath.pdf(t);s+=WadMath.cdf(t);
  }
 }
}
```

## Harness: verbatim Python measurement code

Save this Python block as `scratch/measure.py` inside the pinned setup checkout. Save this report anywhere accessible. After the successful setup from the brief, run:

```sh
. scratch/env.sh && python scratch/measure.py /absolute/path/to/result.md
```

The harness reads the first two Solidity blocks from this report, compiles exactly them, regenerates the erfc coefficients, verifies every grid output against the integer model, measures gas and writes `scratch/metrics.json`. It needs only the setup's packages plus Python's standard library and the already-installed solc 0.8.24. It neither fetches a dependency nor installs a compiler. Its synthetic local EVM caller is only a test account; the library takes no owner, key, deployment address or recipient parameter. All source needed by the math library is embedded in this report.

```python
import math,json,re,sys,hashlib,statistics
from pathlib import Path
import solcx,pyrevm
W=10**18; L=693147180559945309
IS2=707106781186547524; IS2PI=398942280401432678
# The delivered report contains the exact two source blocks compiled here.
if len(sys.argv)>1:
    blocks=re.findall(r'```solidity\n(.*?)\n```',Path(sys.argv[1]).read_text(),re.S)
    source='\n'.join(blocks[:2])
else:
    source=Path('scratch/WadMath.sol').read_text()+Path('scratch/Probe.sol').read_text()
# Parse shipped coefficients, rather than regenerating a different approximation.
body=source.split('function erf(int256 x)',1)[1].split('function pdf',1)[0]
rows=[]
for block in re.findall(r'(?:if|else if) \(j == \d\) \{(.*?)\}',body,re.S):
    rows.append(list(reversed([int(x) for x in re.findall(r'p = (-?\d+)',block)])))
# Use exact Solidity literals, including any future edits to the Pade coefficients.
expbody=source.split('function exp(int256 x)',1)[1].split('function erf',1)[0]
pade=list(reversed([int(v) for v in re.findall(r'p = (\d+)',expbody)]))
assert len(rows)==6 and all(len(r)==19 for r in rows) and len(pade)==7

from decimal import Decimal as D, localcontext
with localcontext() as ctx:
 ctx.prec=95
 pi=D('3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117')
 def cos(x):
  a=s=D(1)
  for n in range(1,180):
   a=-a*x*x/((2*n-1)*(2*n));s+=a
   if abs(a)<D('1e-92'):break
  return s
 def erfc(x):
  # Positive series: erf(x) = 2/sqrt(pi) exp(-x*x) sum x(2x*x)^n/(2n+1)!!.
  a=s=x
  for n in range(1,500):
   a=a*2*x*x/(2*n+1);s+=a
   if abs(a)<D('1e-85'):break
  return 1-2/pi.sqrt()*(-x*x).exp()*s
 generated_rows=[]
 for j in range(6):
  n=19; angles=[pi*(D(i)+D('.5'))/n for i in range(n)]
  ys=[erfc(D(j)+(cos(a)+1)/2) for a in angles]
  cheb=[2*sum(y*cos(k*a) for y,a in zip(ys,angles))/n for k in range(n)];cheb[0]/=2
  # T_0=1, T_1=t, T_k=2tT_(k-1)-T_(k-2).
  polys=[[D(1)],[D(0),D(1)]]
  for k in range(2,n):
   p=[D(0)]*(k+1)
   for i,v in enumerate(polys[-1]):p[i+1]+=2*v
   for i,v in enumerate(polys[-2]):p[i]-=v
   polys.append(p)
  power=[sum(cheb[k]*polys[k][i] for k in range(i,n)) for i in range(n)]
  generated_rows.append([int((v*W).to_integral_value()) for v in power])
assert generated_rows == rows, "shipped coefficients differ from 95-digit generator"

def div(a,b):
    return (-1 if (a<0)!=(b<0) else 1)*(abs(a)//abs(b))

def norm(x):
    if x<=0: raise ValueError('ln domain')
    k=x.bit_length()-1-59
    m=x>>k if k>=0 else x<<(-k)
    if m<W: m*=2; k-=1
    return m,k

def ln(x,alt=False):
    m,k=norm(x)
    if alt:
        f=0
        for i in range(64):
            m=m*m//W
            if m>=2*W: m//=2; f|=1<<(63-i)
        return k*L+f*L//(1<<64)
    if m>1414213562373095048: m//=2; k+=1
    z=div((m-W)*W,m+W); q=z*z//W; p=0
    for i in range(11,-1,-1): p=2*W//(2*i+1)+div(p*q,W)
    return k*L+div(z*p,W)

def reduce(x):
    k=div(x,L); r=x-k*L
    if r>L//2: r-=L;k+=1
    if r<-(L//2): r+=L;k-=1
    return r,k

def exp(x,alt=False):
    if x<=-42*W: return 0
    if x>130*W: raise ValueError('exp domain')
    r,k=reduce(x)
    if alt:
        term=y=W
        for i in range(1,17): term=div(term*r,W*i); y+=term
    else:
        p=q=0
        for c in reversed(pade): p=c+div(p*r,W);q=c-div(q*r,W)
        y=p*W//q
    return y<<k if k>=0 else y>>(-k)

def erf(x,alt=False):
    if x>=6*W: return W
    if x<=-6*W: return -W
    if x==0: return 0
    a=abs(x)
    if alt:
        t=W*W//(W+327591100000000000*a//W)
        p=1061405429000000000
        for c in [-1453152027000000000,1421413741000000000,
                  -284496736000000000,254829592000000000]:
            p=div(p*t,W)+c
        p=div(p*t,W)
        y=W-p*exp(div(-a*a,W))//W
    else:
        j=a//W; t=2*(a-j*W)-W; p=0
        for c in reversed(rows[j]): p=c+div(p*t,W)
        y=W-min(W,max(0,p))
    return -y if x<0 else y

def pdf(x,alt=False):
    if x>=10*W or x<=-10*W: return 0
    return exp(div(-x*x,2*W),alt)*IS2PI//W

def cdf(x,alt=False):
    if x>=9*W: return W
    if x<=-9*W: return 0
    return (W+erf(div(x*IS2,W),alt))//2

names=['ln','exp','erf','pdf','cdf']
funcs=[ln,exp,erf,pdf,cdf]
refs=[math.log,math.exp,math.erf,
      lambda x: math.exp(-x*x/2)/math.sqrt(2*math.pi),
      lambda x: math.erfc(-x/math.sqrt(2))/2]
N=20001
# Log-spaced ln inputs; uniform other grids. Endpoints and zero included.
grids=[[int(round(10**(-6+12*i/(N-1))*W)) for i in range(N)]]
for lo,hi in [(-40,40),(-6,6),(-8,8),(-8,8)]:
    grids.append([lo*W+(hi-lo)*W*i//(N-1) for i in range(N)])
inputs=[
 [10**12,10**14,10**16,5*10**17,W,2*W,10*W,100*W,10**21,10**24],
 [-40*W,-20*W,-5*W,-W,0,W,5*W,10*W,20*W,40*W],
 [-6*W,-4*W,-2*W,-W,0,W//4,W,2*W,4*W,6*W],
 [-8*W,-6*W,-3*W,-W,0,W//2,W,3*W,6*W,8*W],
 [-8*W,-6*W,-3*W,-W,0,W//2,W,3*W,6*W,8*W]]
settings={'optimizer':{'enabled':True,'runs':200},'evmVersion':'cancun',
          'outputSelection':{'*':{'*':['evm.bytecode.object',
           'evm.deployedBytecode.object','evm.methodIdentifiers']}}}
builds=solcx.compile_standard({'language':'Solidity',
       'sources':{'Math.sol':{'content':source}},'settings':settings},
       solc_version='0.8.24')['contracts']['Math.sol']
build=builds['Probe']['evm']
evm=pyrevm.EVM(spec_id='CANCUN'); caller='0x'+'11'*20
# A synthetic local caller is used only for pure-function gas measurement.
evm.set_balance(caller,10**30)
addr=evm.deploy(caller,bytes.fromhex(build['bytecode']['object']))
selector=bytes.fromhex(build['methodIdentifiers']['run(uint256,int256)'])
def call(method,x):
    data=selector+method.to_bytes(32,'big')+x.to_bytes(32,'big',signed=True)
    out=evm.message_call(caller,addr,calldata=data,gas=1000000,is_static=True)
    assert evm.result.is_success,(method,x,evm.result.reason)
    intrinsic=21000+sum(4 if b==0 else 16 for b in data)
    return int.from_bytes(out,'big',signed=True),evm.result.gas_used,evm.result.gas_used-intrinsic

def stats(v): return {'min':min(v),'mean':statistics.fmean(v),'max':max(v)}
report={'compiler':'0.8.24','settings':settings,'n_per_method':N,
        'probe_runtime_bytes':len(build['deployedBytecode']['object'])//2,
        'source_sha256':hashlib.sha256(source.encode()).hexdigest(), 'methods':[]}
for i,(name,fn,ref,xs,gasxs) in enumerate(zip(names,funcs,refs,grids,inputs)):
    for alt in [False,True]:
        method=2*i+int(alt); ae=[]; re=[];core=[];worst=None; grid_peak=0
        # Full-grid EVM equality verifies the integer Python model, for BOTH alternatives.
        for x in xs:
            y=fn(x,alt); actual,_,execution_gas=call(method,x)
            grid_peak=max(grid_peak,execution_gas)
            assert actual==y,('mismatch',name,alt,x,actual,y)
            expected=ref(x/W); err=abs(y/W-expected);ae.append(err)
            if expected!=0:
                relative=err/abs(expected); re.append(relative)
                # Regions where relative error is useful before deep-tail quantization.
                is_core=[abs(expected)>=1e-3,x>=-20*W,abs(x)>=W//1000,
                         abs(x)<=6*W,x>=-6*W][i]
                if is_core: core.append(relative)
            if worst is None or err>worst[0]:worst=(err,x/W)
        raw=[]; execution=[]
        for x in gasxs:
            y,g,h=call(method,x); assert y==fn(x,alt)
            raw.append(g);execution.append(h)
        row={'name':name,'alternative':alt,'max_abs':max(ae),
             'mean_abs':statistics.fmean(ae),'max_rel':max(re),
             'mean_rel':statistics.fmean(re),'worst_abs_at':worst[1],
             'core_max_rel':max(core),'core_mean_rel':statistics.fmean(core),
             'grid_peak_execution_gas':grid_peak,
             'tx_gas':stats(raw),'execution_gas':stats(execution),
             'gas_inputs':[x/W for x in gasxs],'tx_samples':raw,
             'execution_samples':execution}
        report['methods'].append(row)
        print(json.dumps(row),flush=True)
# Extremes and saturation/rejection boundaries beyond the error grids.
I=(1<<255)-1; boundary_count=0
for i,fn in enumerate(funcs):
    cases=([1,2,W-1,W,W+1,I] if i==0 else
           [-(1<<255),-42*W-1,-42*W,-42*W+1,130*W-1,130*W] if i==1 else
           [-(1<<255),I,-10*W,-9*W,-6*W,-1,0,1,6*W,9*W,10*W])
    for alt in [False,True]:
        for x in cases:
            assert call(i*2+int(alt),x)[0]==fn(x,alt)
            boundary_count+=1
for method,x in [(0,0),(0,-1),(1,0),(1,-1),(2,130*W+1),(3,130*W+1)]:
    try: call(method,x)
    except (AssertionError,RuntimeError): pass
    else: raise AssertionError(('expected revert',method,x))
report['boundary_success_cases']=boundary_count
report['expected_revert_cases']=6
# Saturation discontinuities and interval boundaries in integer WADs.
join_count=0
for i in [2,4]:
    for alt in [False,True]:
        for j in range(-9,10):
            for delta in [-1,0,1]:
                x=j*W+delta
                assert call(i*2+int(alt),x)[0]==funcs[i](x,alt)
                join_count+=1
report['join_cases']=join_count
# A separate contract measures five calls to each selected routine in one transaction.
b=builds['SelectedProbe']['evm']
report['selected_probe_runtime_bytes']=len(b['deployedBytecode']['object'])//2
bundle_evm=pyrevm.EVM(spec_id='CANCUN');bundle_evm.set_balance(caller,10**30)
bundle_addr=bundle_evm.deploy(caller,bytes.fromhex(b['bytecode']['object']))
bundle_samples=[]
for v in [W//10,W//5,W//2,W,2*W,3*W,4*W,5*W,6*W,8*W]:
    d=bytes.fromhex(b['methodIdentifiers']['bundle(int256)'])+v.to_bytes(32,'big')
    bundle_evm.message_call(caller,bundle_addr,calldata=d,gas=245000,is_static=True)
    assert bundle_evm.result.is_success
    g=bundle_evm.result.gas_used
    bundle_samples.append({'input':v/W,'tx':g,
        'execution':g-21000-sum(4 if z==0 else 16 for z in d)})
report['bundle_samples']=bundle_samples
report['bundle_execution_gas']=stats([v['execution'] for v in bundle_samples])
report['bundle_tx_gas']=stats([v['tx'] for v in bundle_samples])
Path('scratch').mkdir(exist_ok=True)
Path('scratch/metrics.json').write_text(json.dumps(report,indent=2))
print('PASS',json.dumps({k:v for k,v in report.items() if k!='methods'}),flush=True)
```

## Problems and limits

- No setup, compilation or runtime blocker remained. Final measurements include exact Python/EVM equality for both approaches, including signed truncation and saturation branches.
- Error figures are empirical comparisons to binary64 `math`, not certified worst-case bounds. Near a zero reference, relative error is ill-conditioned; near extreme tails, WAD quantization limits relative accuracy. In particular, exp on [-40,40] has roughly 19.8% maximum relative error at tiny positive values even though it is much more accurate in absolute terms there. Both exp approaches share that limitation.
- Piecewise polynomial quantization can cause a few raw-unit artifacts at joins and near zero. Exact-zero and range clamps are included; no formal proof of monotonicity or exhaustive verification over all int256 values is claimed. Pdf and cdf are not specialized log-density/log-probability routines.
- Exp's safe upper input is intentionally restricted to 130; callers must keep it in range. Gaussian routines guard their own inputs. No wider-domain accuracy sweep than those stated was performed.
- The gas figures depend on compiler settings, EVM revision and wrapper dispatch. They do not establish the complete strategy's 245,000-gas guarantee, 0 OutOfGas, bytecode budget, fee parity, or improvement over t03. No real-engine `amm_sim_rs.run_batch`, paired score, `gas_stats()`, replica policy compression test or P/Q/F gate was run for this standalone math item. Those remain integration work for the policy strategy.

External provenance is limited to Cook's A&S comparison implementation and NIST's erf definition, linked above and quoted in `sources.json`. All coefficients for the selected polynomial, all error/gas figures and bytecode sizes are generated locally by the included code. The pinned engine setup and goal reward targets came from the task brief and are not independently asserted research results.