# Q64.64 math for the AMM callback

Use `Q64Math.ln`, `exp`, `erf`, `pdf`, and `cdf` below. The alternatives are included for reproducible comparison. All numerical results in this report are original measurements of the compiled Solidity, called from Python through pyrevm, rather than estimates from a floating-point implementation.

## Format and safe domains

Arguments and results are signed `int128` Q64.64: real value = raw integer / 18446744073709551616. Representable values are [-2^63, 2^63 - 2^-64]; one unit is approximately 5.421010862427522e-20. Intermediates use `int256`. Encode off-chain by truncating x*2^64 toward zero. No storage, dependencies, linked libraries, assembly, or external calls are required by the library.

| Functions | Safe non-reverting input domain (real units) | Boundary behavior |
|---|---|---|
| ln, lnBits | Every positive int128, from 2^-64 to 2^63 - 2^-64 | Reject zero and negative input |
| exp, expTaylor | [-40,40], inclusive | Explicitly reject outside; no overflow inside |
| erf, erfQuad | Every int128 | Saturate to sign(x) for abs(x) >= 6 |
| pdf, pdfTaylor | [-8,8], inclusive | Explicitly reject outside |
| cdf, cdfQuad | Every int128 | Selected cdf saturates when the scaled magnitude reaches 6; cdfQuad also bounds before its int128 conversion |

These are arithmetic safety domains, not claims of uniformly bounded approximation error over all representable inputs. Accuracy domains measured below are ln [1e-6,1e6], exp [-40,40], erf [-6,6], pdf and cdf [-8,8]. lnBits, expTaylor, erfQuad, pdfTaylor and cdfQuad share the respective format and safety rules.

Safety reasoning: ln's normalized mantissa is in [1/sqrt(2),sqrt(2)], and its exponent is between -64 and 62 before the possible normalization increment. Its products fit int256. exp reduces to abs(r) <= ln(2)/2 with abs(k) <= 58; both polynomial and Padé intermediates fit int256 and the largest output is about 2.354e17, below 2^63. The Padé denominator stays positive on that interval. erf promotes before negating, so int128.min is safe, and bounds the argument before squaring; pdf checks its domain before squaring. cdf's largest multiplication has fewer than 192 magnitude bits; tail bounds before squaring. The quadrature integrands have exponent in [-36,0]. `mul` is private and only used with these bounded intermediates, not offered as a general overflow-safe multiply.

Signed right shifts round downward, while division rounds toward zero; the implementation intentionally uses both. Shift overflow has no automatic check: the domain guards and intermediate bounds above are necessary. See [Solidity 0.8.24 types](https://docs.soliditylang.org/en/v0.8.24/types.html).

## Measurement protocol

Compiler: solc 0.8.24, optimizer enabled, runs=200, viaIR=true, EVM Shanghai. Python 3.12.3, pyrevm 0.3.7, py-solc-x 2.1.0, NumPy 2.5.4. The required setup at amm-challenge commit 6b6921f8d463314aa7c96234e1be25a25245f704 completed, including `amm-match validate contracts/src/StarterStrategy.sol`.

Each implementation is evaluated on 10,001 inputs: geometric spacing for ln and linear spacing for all others, including endpoints. Every result comes from actual compiled EVM execution. Inputs are encoded to Q64.64 first; the baseline is evaluated at that decoded input. Baselines are `math.log`, `math.exp`, `math.erf`, `math.exp(-x*x/2)/sqrt(2*pi)`, and `0.5*math.erfc(-x/sqrt(2))`. Python documents these functions and the usefulness of erfc in avoiding subtraction near one in [the math reference](https://docs.python.org/3/library/math.html).

Absolute error is abs(decoded_output - baseline). Relative error divides by abs(baseline); exact-zero baselines are excluded only from relative aggregates. Tables report both all nonzero baselines and a separate restriction abs(baseline) >= 1e-8 (marked R8). The latter avoids presenting near-zero relative errors as ordinary operating accuracy. Values are comparisons in binary64, not correctly-rounded Q64.64 proofs. Numerical maxima are sampled maxima, not certified global bounds.

Gas uses ten inputs for each function, without an RPC or fork. Inner gas is `gasleft()` before dispatch to the internal routine minus `gasleft()` after it; it includes the small function-selection branch and measurement overhead, but excludes ABI decoding, return encoding and transaction intrinsic gas. Total gas is separately read from `evm.result.gas_used` in pyrevm, including that call's intrinsic cost. No access-list or storage costs are involved. Gas depends on compiler settings and call site; these numbers are reproducible wrapper costs, not engine callback measurements.

Ten gas inputs per function family:

- ln: 1e-6, .001, .1, .5, .99, 1, 1.01, 2, 1000, 1e6.
- exp: -40, -32, -10, -1, 0, .1, 1, 10, 32, 40.
- erf: -6, -4, -2, -1, -.01, 0, .01, 1, 3, 6.
- pdf and cdf: -8, -6, -3, -1, -.01, 0, .01, 1, 3, 8.

## ln: alternatives and selection

Select range reduction plus eight terms of 2*(z+z^3/3+...+z^15/15), z=(m-1)/(m+1). The mantissa reduction keeps abs(z) <= about .171573. This is cheaper than the alternative's 48 successive squarings to recover binary log2 digits. lnBits has a lower sampled maximum absolute error, but worse mean absolute error and substantially higher gas. Near ln(1)=0 its finite binary digit count also makes relative error misleading.

| Implementation | Max abs | Mean abs | Max rel, nonzero | Mean rel, nonzero | Max R8 | Mean R8 |
|---|---:|---:|---:|---:|---:|---:|
| ln | 1.24345e-14 | 6.608192e-16 | 3.198353e-14 | 2.369811e-16 | 3.198353e-14 | 2.369061e-16 |
| lnBits | 3.552714e-15 | 1.228004e-15 | 0.2322862 | 2.32263e-05 | 7.28896e-13 | 8.245634e-16 |

| Implementation | Inner min | Inner mean | Inner max | Total min | Total mean | Total max |
|---|---:|---:|---:|---:|---:|---:|
| ln | 1505 | 1601.2 | 1713 | 23115 | 23250.8 | 23395 |
| lnBits | 4824 | 5763 | 6325 | 26446 | 27424.6 | 27971 |

Selected routine's worst absolute-error sample: x=5.395106e-06, output=-12.130018269892629, baseline=-12.130018269892641.

## exp: alternatives and selection

Select a diagonal [5/5] Padé approximation after nearest-ln(2) range reduction, then shift by 2^k. Numerator coefficients, ascending, are 30240,15120,3360,420,30,1; denominator is the same polynomial at -r. They result from matching the exponential power series through degree ten. The alternative uses twelve Taylor terms after the same reduction. Padé sacrifices a small amount of floating-point-scale accuracy for much lower gas. Use expTaylor if absolute errors of tens or hundreds at outputs around 1e17 matter to the application.

| Implementation | Max abs | Mean abs | Max rel, nonzero | Mean rel, nonzero | Max R8 | Mean R8 |
|---|---:|---:|---:|---:|---:|---:|
| exp | 192 | 0.2601508 | 0.01263029 | 8.271384e-05 | 4.729183e-12 | 4.60889e-14 |
| expTaylor | 32 | 0.04101158 | 0.01263029 | 8.271384e-05 | 4.729183e-12 | 4.603901e-14 |

| Implementation | Inner min | Inner mean | Inner max | Total min | Total mean | Total max |
|---|---:|---:|---:|---:|---:|---:|
| exp | 650 | 664.8 | 699 | 22260 | 22403.2 | 22597 |
| expTaylor | 2160 | 2174.4 | 2211 | 23770 | 23912.8 | 24109 |

Selected routine's worst absolute-error sample: x=39.856, output=2.038172186278239e+17, baseline=2.038172186278241e+17.

## erf: alternatives and selection

Select erfc(a) = exp(-a*a)*P(1/(1+a)), a=abs(x), followed by sign restoration. P is a degree-16 least-squares polynomial fitted directly to exp(a*a)*math.erfc(a) on [0,6], not a borrowed coefficient table. Its integer coefficients are adjusted to sum exactly to Q, making erf(0)=0. Clamp the complementary tail to [0,1]. The alternative integrates exp(-t*t) using sixteen Gauss-Legendre nodes and symmetry. Although quadrature has a smaller sampled maximum absolute error, it is much more expensive and was not monotone on the measurement grid near saturation. The selected erf was monotone on that grid; this is not a proof of monotonicity at every raw input.

| Implementation | Max abs | Mean abs | Max rel, nonzero | Mean rel, nonzero | Max R8 | Mean R8 |
|---|---:|---:|---:|---:|---:|---:|
| erf | 6.713962e-11 | 1.212119e-12 | 0.0004130241 | 4.132701e-08 | 1.3415e-08 | 2.873813e-11 |
| erfQuad | 1.56849e-11 | 1.290909e-12 | 7.379493e-05 | 7.380046e-09 | 1.56849e-11 | 1.291042e-12 |

| Implementation | Inner min | Inner mean | Inner max | Total min | Total mean | Total max |
|---|---:|---:|---:|---:|---:|---:|
| erf | 291 | 1154.7 | 1404 | 21913 | 22931.5 | 23349 |
| erfQuad | 306 | 8389.1 | 10684 | 21928 | 30165.9 | 32582 |

Selected routine's worst absolute-error sample: x=-0.0108, output=-0.012186021277425171, baseline=-0.012186021210285549.

## pdf: alternatives and selection

Select exp(-x*x/2)/sqrt(2*pi) using the Padé exp. Compare the same formula with the Taylor exp, so the competing exponential algorithms are isolated. The Taylor version is more accurate at central values by a few binary64 ulps, but costs more. Both lose relative resolution far into the tails because the intermediate and final results are Q64.64.

| Implementation | Max abs | Mean abs | Max rel, nonzero | Mean rel, nonzero | Max R8 | Mean R8 |
|---|---:|---:|---:|---:|---:|---:|
| pdf | 2.775558e-16 | 6.248681e-18 | 1.252073e-05 | 1.225037e-07 | 6.203134e-12 | 1.11948e-13 |
| pdfTaylor | 5.551115e-17 | 3.847469e-18 | 1.252073e-05 | 1.225037e-07 | 6.203134e-12 | 1.119205e-13 |

| Implementation | Inner min | Inner mean | Inner max | Total min | Total mean | Total max |
|---|---:|---:|---:|---:|---:|---:|
| pdf | 963 | 984.9 | 1012 | 22573 | 22761.7 | 22957 |
| pdfTaylor | 2473 | 2495.1 | 2524 | 24083 | 24271.9 | 24467 |

Selected routine's worst absolute-error sample: x=-0.832, output=0.28222503404476229, baseline=0.28222503404476257.

## cdf: alternatives and selection

Select the complementary tail directly: t=erfc(abs(x)/sqrt(2))/2; return t for negative x and 1-t for nonnegative x. This avoids cancellation in the negative tail. Compare numerical quadrature followed by (1+erf(x/sqrt(2)))/2. The latter has a smaller maximum absolute error on this grid, but its absolute quadrature error overwhelms tiny negative-tail probabilities. It also failed sampled monotonicity. Its very large relative error is a reason to reject it despite the lower absolute error.

| Implementation | Max abs | Mean abs | Max rel, nonzero | Mean rel, nonzero | Max R8 | Mean R8 |
|---|---:|---:|---:|---:|---:|---:|
| cdf | 3.363132e-11 | 6.428592e-13 | 7.014483e-05 | 3.796793e-07 | 6.804441e-11 | 6.797257e-12 |
| cdfQuad | 7.842429e-12 | 5.645441e-13 | 11631.68 | 92.98573 | 2.247612e-07 | 1.160978e-09 |

| Implementation | Inner min | Inner mean | Inner max | Total min | Total mean | Total max |
|---|---:|---:|---:|---:|---:|---:|
| cdf | 1459 | 1481.4 | 1509 | 23070 | 23258.2 | 23453 |
| cdfQuad | 10475 | 10737.6 | 11009 | 32085 | 32514.4 | 32907 |

Selected routine's worst absolute-error sample: x=-0.0144, output=0.49425542966162417, baseline=0.49425542969525549.

## Complete code

The following single source compiles as written. All five selected entry points and all five measured alternatives are `internal pure`. The compiler eliminates unused alternatives when the library is inlined into a strategy. Extract this block as Q64Math.sol, import it, and call `Q64Math.ln(x)` etc. The harness appends its own thin wrapper; no `AMMStrategyBase` or policy logic belongs in this math-only deliverable.

```solidity
// SPDX-License-Identifier: MIT
pragma solidity 0.8.24;
// Q64.64 math for the Taylor (gtmcknight) AMM policy port.
// Independently derived numerical routines; no Taylor policy code is included.
library Q64Math {
    int256 internal constant Q = 18446744073709551616;
    int256 internal constant LN2 = 12786308645202655659;
    int256 internal constant SQRT2 = 26087635650665564424;
    function mul(int256 a, int256 b) private pure returns(int256) {
        unchecked { return (a*b)>>64; }
    }
    function msb(uint256 a) private pure returns(uint256 n) {
        unchecked {
            if(a>=1<<64){a>>=64;n+=64;}
            if(a>=1<<32){a>>=32;n+=32;}
            if(a>=1<<16){a>>=16;n+=16;}
            if(a>=1<<8){a>>=8;n+=8;}
            if(a>=1<<4){a>>=4;n+=4;}
            if(a>=1<<2){a>>=2;n+=2;}
            if(a>=1<<1){n+=1;}
        }
    }
    function ln(int128 x) internal pure returns(int128) {
        require(x>0,"ln domain");
        unchecked {
            int256 k=int256(msb(uint256(int256(x))))-64;
            int256 m=k>=0?int256(x)>>uint256(k):int256(x)<<uint256(-k);
            if(m>SQRT2){m>>=1;++k;}
            int256 z=((m-Q)<<64)/(m+Q);
            int256 z2=mul(z,z); int256 term=z; int256 sum=z;
            for(int256 j=3;j<=15;j+=2){term=mul(term,z2);sum+=term/j;}
            return int128(k*LN2+2*sum);
        }
    }
    // Alternative: 48 binary fractional log2 digits by repeated squaring.
    function lnBits(int128 x) internal pure returns(int128) {
        require(x>0,"ln domain");
        unchecked {
            int256 k=int256(msb(uint256(int256(x))))-64;
            int256 m=k>=0?int256(x)>>uint256(k):int256(x)<<uint256(-k);
            int256 f;
            for(uint256 i=1;i<=48;++i){
                m=mul(m,m);
                if(m>=2*Q){m>>=1;f+=Q>>i;}
            }
            return int128(k*LN2+mul(f,LN2));
        }
    }
    function expTaylor(int128 x) internal pure returns(int128) {
        require(x>=-40*Q && x<=40*Q,"exp domain");
        unchecked {
            int256 k=int256(x)/LN2;
            int256 r=int256(x)-k*LN2;
            if(r>LN2/2){r-=LN2;++k;}
            if(r< -LN2/2){r+=LN2;--k;}
            int256 term=Q;int256 sum=Q;
            for(int256 j=1;j<=12;++j){term=mul(term,r)/j;sum+=term;}
            return int128(k>=0?sum<<uint256(k):sum>>uint256(-k));
        }
    }
    // Alternative: diagonal [5/5] Pade on the same reduced interval.
    function exp(int128 x) internal pure returns(int128) {
        require(x>=-40*Q && x<=40*Q,"exp domain");
        unchecked {
            int256 k=int256(x)/LN2;int256 r=int256(x)-k*LN2;
            if(r>LN2/2){r-=LN2;++k;}
            if(r< -LN2/2){r+=LN2;--k;}
            int256 a=Q; int256 b=Q;
            a=mul(a,r)+30*Q; b=mul(b,-r)+30*Q;
            a=mul(a,r)+420*Q; b=mul(b,-r)+420*Q;
            a=mul(a,r)+3360*Q; b=mul(b,-r)+3360*Q;
            a=mul(a,r)+15120*Q; b=mul(b,-r)+15120*Q;
            a=mul(a,r)+30240*Q; b=mul(b,-r)+30240*Q;
            int256 v=(a<<64)/b;
            return int128(k>=0?v<<uint256(k):v>>uint256(-k));
        }
    }
    // Positive-side erfc, also used directly in the negative CDF tail.
    function tail(int256 a) private pure returns(int256) {
        unchecked {
            if(a>=6*Q) return 0;
            int256 t=(Q<<64)/(Q+a);
            int256 p=10577369877819737502;
            p=mul(p,t)+(-94265640308351673045);
            p=mul(p,t)+(378466513249192706555);
            p=mul(p,t)+(-894905970504598549448);
            p=mul(p,t)+(1355237508470430716407);
            p=mul(p,t)+(-1305662266791571766152);
            p=mul(p,t)+(666283532756899319825);
            p=mul(p,t)+(84016965711370517910);
            p=mul(p,t)+(-427921320700291957535);
            p=mul(p,t)+(331518931392328375497);
            p=mul(p,t)+(-117616370305056214833);
            p=mul(p,t)+(16626067415089017687);
            p=mul(p,t)+(-10540927982179093740);
            p=mul(p,t)+(5869848907746932947);
            p=mul(p,t)+(10352386200736936539);
            p=mul(p,t)+(10410177180298953339);
            p=mul(p,t)+(-60496154407839);

            int256 v=mul(p,exp(int128(-mul(a,a))));
            if(v<0)return 0;if(v>Q)return Q;return v;
        }
    }
    function erf(int128 x) internal pure returns(int128) {
        unchecked {
            int256 a=x; bool neg=a<0;if(neg)a=-a;
            int256 v=Q-tail(a);return int128(neg?-v:v);
        }
    }
    // Alternative: 16-point Gauss-Legendre integral of exp(-t*t).
    function erfQuad(int128 x) internal pure returns(int128) {
        unchecked {
            int256 a=x;bool neg=a<0;if(neg)a=-a;
            if(a>=6*Q)return int128(neg?-Q:Q);
            int256 sum;int256 u;int256 v;
            u=mul(a,10099707763242631743);v=mul(a,8347036310466921717);
            sum+=mul(3494746925672694254,int256(exp(int128(-mul(u,u))))+int256(exp(int128(-mul(v,v)))));
            u=mul(a,11820706352591207180);v=mul(a,6626037721118345173);
            sum+=mul(3368438464319070733,int256(exp(int128(-mul(u,u))))+int256(exp(int128(-mul(v,v)))));
            u=mul(a,13447831176308776418);v=mul(a,4998912897400773352);
            sum+=mul(3120387021679099961,int256(exp(int128(-mul(u,u))))+int256(exp(int128(-mul(v,v)))));
            u=mul(a,14922274511714967823);v=mul(a,3524469561994583977);
            sum+=mul(2759558920152906600,int256(exp(int128(-mul(u,u))))+int256(exp(int128(-mul(v,v)))));
            u=mul(a,16190747933393138497);v=mul(a,2255996140316414963);
            sum+=mul(2298998736920541162,int256(exp(int128(-mul(u,u))))+int256(exp(int128(-mul(v,v)))));
            u=mul(a,17207410663187680669);v=mul(a,1239333410521871315);
            sum+=mul(1755364711442041502,int256(exp(int128(-mul(u,u))))+int256(exp(int128(-mul(v,v)))));
            u=mul(a,17935538891379884029);v=mul(a,511205182329667586);
            sum+=mul(1148374823782780663,int256(exp(int128(-mul(u,u))))+int256(exp(int128(-mul(v,v)))));
            u=mul(a,18348984953894730295);v=mul(a,97759119814821892);
            sum+=mul(500874469740415484,int256(exp(int128(-mul(u,u))))+int256(exp(int128(-mul(v,v)))));
            
            int256 out=mul(mul(a,sum),10407460856758233470);
            if(out>Q)out=Q;return int128(neg?-out:out);
        }
    }
    function pdf(int128 x) internal pure returns(int128) {
        require(x>=-8*Q && x<=8*Q,"pdf domain");
        return int128(mul(exp(int128(-mul(x,x)/2)),7359186146747302859));
    }
    function pdfTaylor(int128 x) internal pure returns(int128) {
        require(x>=-8*Q && x<=8*Q,"pdf domain");
        return int128(mul(expTaylor(int128(-mul(x,x)/2)),7359186146747302859));
    }
    function cdf(int128 x) internal pure returns(int128) {
        unchecked {
            int256 a=x;bool neg=a<0;if(neg)a=-a;
            int256 t=tail(mul(a,13043817825332781762))/2;
            return int128(neg?t:Q-t);
        }
    }
    function cdfQuad(int128 x) internal pure returns(int128) {
        // Bound before conversion so every int128 input is safe.
        if(x>=12*Q)return int128(Q);if(x<=-12*Q)return 0;
        return int128((Q+int256(erfQuad(int128(mul(x,13043817825332781762)))))/2);
    }
}
```

## Recommendation

Use the five unsuffixed functions. Their measured inner maxima sum to 6,337 gas for one call to each, or 31,685 for five calls to each, by arithmetic on the measured maxima. That leaves considerable room within a 245,000-gas callback for filtering, tables, and fee calculation. This is a planning estimate; a combined callback, with its actual branch paths, compiler output and data access, must still be measured on the real engine. The ten inputs are representative observations and do not establish a universal gas bound.

The complete benchmark wrapper, including all alternatives, has 2921 bytes of deployed runtime. This is below 24,576 bytes; the enclosing strategy's deployed size must be checked independently. Selected ln, exp, erf and cdf were monotone on their respective measurement grids. Endpoints and rejection checks in the harness also passed, including both int128 extremes for saturating functions and raw=1 / int128.max for ln.

Choose the Taylor exponential only when its extra accuracy is useful. Do not use quadrature for this callback: it spends more gas and its tail error is less useful than the fitted complementary form. Do not infer a strategy's fee accuracy or score from the math errors alone: amplification by divisions, filtering and interpolation still matters.

## Harness (verbatim)

Save the following Python block as `measure.py` under the setup checkout's scratch directory. With the required venv and solc already installed, run:

```sh
. scratch/env.sh && python scratch/measure.py /path/to/result.md 0 5 > scratch/first.json
. scratch/env.sh && python scratch/measure.py /path/to/result.md 5 10 > scratch/second.json
```

Each invocation compiles locally and measures five functions, well below the eight-minute command limit in this run. Omitting the final two arguments measures all ten. It performs no installation, network access, or reads of machine configuration. The source SHA-256 measured here is `549ea0f2a21c08083c14b186b7cacc225542a7177e655822bfacbed8a9186fea` (contents inside the Solidity block, excluding its Markdown fences). The wrapper's test-only caller address is local EVM scaffolding and has no production meaning.

```python
"""Run inside the requested venv; python measure.py result.md [first last].
Without .md, read the Solidity source path directly. Outputs JSON to stdout.
Dependencies: solcx 2.1.0, pyrevm 0.3.7, numpy (grids only).
"""
import sys, re, json, math, statistics, hashlib
from pathlib import Path
import numpy as np
import solcx
from pyrevm import EVM
path=Path(sys.argv[1]); text=path.read_text()
source=re.search(r'```solidity\n(.*?)```',text,re.S)[1] if path.suffix=='.md' else text
names=['ln','lnBits','exp','expTaylor','erf','erfQuad','pdf','pdfTaylor','cdf','cdfQuad']
wrapper='''\ncontract Bench {
function run(uint8 f,int128 x) external view returns(int128 y,uint256 used){
uint256 start=gasleft();
'''
for i,n in enumerate(names):
 wrapper+=('if' if i==0 else 'else if')+f'(f=={i}) y=Q64Math.{n}(x);\n'
wrapper+='else revert(); used=start-gasleft(); }}'
compiled=solcx.compile_standard({'language':'Solidity','sources':{'M.sol':{'content':source+wrapper}},'settings':{'optimizer':{'enabled':True,'runs':200},'viaIR':True,'evmVersion':'shanghai','outputSelection':{'*':{'*':['evm.bytecode.object','evm.deployedBytecode.object','evm.methodIdentifiers']}}}},solc_version='0.8.24')
c=compiled['contracts']['M.sol']['Bench']['evm']; selector=bytes.fromhex(c['methodIdentifiers']['run(uint8,int128)'])
e=EVM(spec_id='SHANGHAI'); caller='0x'+'10'*20
address=e.deploy(deployer=caller,code=bytes.fromhex(c['bytecode']['object']),gas=10000000)
assert e.result.is_success
Q=1<<64

def call(i,raw):
 data=selector+i.to_bytes(32,'big')+int(raw).to_bytes(32,'big',signed=True)
 out=e.message_call(caller=caller,to=address,calldata=data,gas=1000000)
 assert e.result.is_success,(i,raw,e.result.reason)
 assert len(out)==64
 return int.from_bytes(out[:32],'big',signed=True),int.from_bytes(out[32:],'big'),e.result.gas_used

def target(name,x):
 if name.startswith('ln'):return math.log(x)
 if name.startswith('exp'):return math.exp(x)
 if name.startswith('erf'):return math.erf(x)
 if name.startswith('pdf'):return math.exp(-x*x/2)/math.sqrt(2*math.pi)
 return .5*math.erfc(-x/math.sqrt(2))

def stats(a):return {'min':min(a),'mean':statistics.fmean(a),'max':max(a)}

result={'solc':'0.8.24 optimizer=200 viaIR=true shanghai','source_sha256':hashlib.sha256(source.encode()).hexdigest(),'runtime_bytes':len(c['deployedBytecode']['object'])//2,'functions':{}}
first=int(sys.argv[2]) if len(sys.argv)>2 else 0
last=int(sys.argv[3]) if len(sys.argv)>3 else 10
for i in range(first,last):
 name=names[i]
 if name.startswith('ln'):grid=np.geomspace(1e-6,1e6,10001); reps=[1e-6,.001,.1,.5,.99,1,1.01,2,1000,1e6]
 elif name.startswith('exp'):grid=np.linspace(-40,40,10001);reps=[-40,-32,-10,-1,0,.1,1,10,32,40]
 elif name.startswith('erf'):grid=np.linspace(-6,6,10001);reps=[-6,-4,-2,-1,-.01,0,.01,1,3,6]
 else:grid=np.linspace(-8,8,10001);reps=[-8,-6,-3,-1,-.01,0,.01,1,3,8]
 errors=[];relative=[]; meaningful=[]; worst=None; previous=None; monotone=True
 for requested in grid:
  raw=int(float(requested)*Q);x=raw/Q; y,_,_=call(i,raw); y/=Q; exact=target(name,x)
  err=abs(y-exact);errors.append(err)
  if exact!=0:relative.append(err/abs(exact))
  if abs(exact)>=1e-8:meaningful.append(err/abs(exact))
  if not name.startswith('pdf'):
   if previous is not None and y<previous:monotone=False
   previous=y
  if worst is None or err>worst['abs']:worst={'x':x,'actual':y,'expected':exact,'abs':err}
 gas=[];inner=[];rawgas=[]
 for x in reps:
  _,g,total=call(i,int(x*Q));inner.append(g);gas.append(total);rawgas.append([x,g,total])
 result['functions'][name]={'n':len(errors),'max_abs':max(errors),'mean_abs':statistics.fmean(errors),'max_rel_nonzero':max(relative),'mean_rel_nonzero':statistics.fmean(relative),'max_rel_target_ge_1e-8':max(meaningful),'mean_rel_target_ge_1e-8':statistics.fmean(meaningful),'monotone_grid':monotone,'worst_abs':worst,'inner_gas':stats(inner),'transaction_gas':stats(gas),'gas_inputs':rawgas}
 print(name,result['functions'][name],file=sys.stderr)
# Boundary and rejection checks, independently of the measurement grid.
for i in [0,1]:
 call(i,1);call(i,(1<<127)-1)
for i in [4,5,8,9]:
 call(i,-(1<<127));call(i,(1<<127)-1)
for i,raw in [(0,0),(1,-1),(2,40*Q+1),(3,-40*Q-1),(6,8*Q+1),(7,-8*Q-1)]:
 try:call(i,raw)
 except (AssertionError,RuntimeError):pass
 else:raise AssertionError(('expected rejection',i,raw))
print(json.dumps(result,indent=2))
```

## Coefficient construction

The fitted coefficients are frozen in the complete Solidity source, so reproducing the benchmark does not refit anything. For transparency, this is the fitting procedure used with NumPy 2.5.4. `ci` is in ascending polynomial order; emit it in reverse order for Horner evaluation. The exact-Q sum adjustment was applied after the floating-point endpoint correction.

```python
import math
from decimal import Decimal, getcontext
import numpy as np
from numpy.polynomial import Polynomial
getcontext().prec = 80
Q = 1 << 64
C = lambda a: int(Decimal(str(a)) * Q)
x = np.linspace(0, 6, 20001)
t = 1 / (1 + x)
y = np.array([math.erfc(float(v))*math.exp(float(v*v)) for v in x])
p = Polynomial.fit(t, y, 16).convert()
c = p.coef.copy()
c[-1] += 1 - sum(c)
ci = [C(a) for a in c]
ci[0] += Q - sum(ci)
print(ci)
# Quadrature constants: positive nodes/weights, with pair symmetry.
nodes, weights = np.polynomial.legendre.leggauss(16)
for n, w in zip(nodes[8:], weights[8:]):
    print(C((1+n)/2), C((1-n)/2), C(w))
```

LN2 and SQRT2 were calculated with Decimal at precision 80 and truncated to Q64.64. Inverse square-root constants were calculated with Python math and then truncated; their binary64 rounding is included in the measured error. The uncorrected floating-point polynomial fit had maximum absolute residual 5.3164583846410096e-11 on its 20,001-point fitting grid. That residual is for scaled erfc, not a substitute for the separately measured fixed-point erf errors.

Validation: extracting the frozen Solidity and harness from this delivered Markdown and rerunning all ten functions reproduced all 100,010 accuracy calls and 100 gas calls exactly, including the runtime byte count and boundary/rejection checks. Both sources.json quotations were also checked as exact substrings of fresh HTTP responses from their listed URLs.

## Problems and limits

- At exp(-40), the true output is only about 78 Q64.64 units. Flooring a fraction of one unit can therefore produce roughly one-percent relative error. The measured maximum relative error over the exp grid is about 1.263%, for both approaches. This is primarily a format-resolution limitation; the central approximation's ulp-scale error is much smaller. The absolute-error table is dominated by very large positive outputs.
- At normal-tail probabilities near x=-8, Q64.64 rounding of exp and the subsequent multiplication/division dominates relative error. Selected pdf and cdf have sampled maximum relative errors about 1.252e-5 and 7.014e-5 respectively. A policy depending on relative accuracy of tinier probabilities should carry more fractional bits internally.
- Grids produced by NumPy can include numbers very close to, but not exactly, zero. This explains the large nonzero-relative log/erf figures compared with their absolute error and R8 columns. Exact erf(0)=0, ln(1)=0 and cdf(0)=1/2 are preserved by the selected code.
- Saturation deliberately discards tiny tails. At erf magnitude 6 the omitted complementary tail is about 2.15e-17. No broad uniform-error certification, exact rounding claim, or proof of monotonicity is provided. Refit results can differ with numerical-library versions; benchmark the frozen integer constants.
- The library does not implement or score a complete policy, compress tables, or execute gates P/Q/F, paired seed batches, engine gas_stats, or a 50x2,000-callback fee comparison. Those are integration gates for the other strategy work, not results established by this bounded math deliverable. No improvement over t03 is claimed.
- Scope was limited to the requested math library and measurements. Setup and measurement scaffolding stayed under test/scratch; only result.md and sources.json are delivered. No requester input was missing and no setup step failed.