# r1-port / port-1: Taylor's Gaussian + plain policy (real.py) ported to Solidity

**Bottom line.** The port matches the Python reference to 6.6e-5 bps at most over 100,050 callbacks (Gate F passes). It is **not deployable yet**: on the real engine **290,102 of 536,430 afterSwap calls (54%) ran out of gas**, and it scores **−101.7 (se 7.7) against t03** on the 64 assigned seeds. One finding matters for every copy of this port: **an `afterSwap` call that runs out of gas is reverted, so the strategy's state stops moving**. After the first OutOfGas, every later call starts from the same stale state and usually runs out of gas again. Gas is the work still to do (see Problems). Policy, filter and tables are by Taylor (gtmcknight), github.com/gtmcknight/amm-sim @ 4812a6382311624e1506aca131ca35c378f930d9, research/solved. The contract header credits him.

## Setup (per step)

| step | result |
|---|---|
| clone gtmcknight/amm-sim, checkout 4812a638… | ok; `research/solved` has README.md, real.py, pol_mid.json, gates.py, replica2.py, evalpol.py, … |
| Gate P: `python3 -I gates.py` (before any change) | **reproduced to the digit.** gaussian+plain column: 10000 **443.455598**, 10001 **571.463409**, 10002 **525.168028**, the same as the README's gates table. The other columns matched too (perfect-info 448.261189 / 576.779199 / 531.886497). 8.8 s |
| clone benedictbrady/amm-challenge, checkout 6b6921f8… | ok |
| `sh scratch/setup.sh` (verbatim from the brief) | ok. The wheel's sha256 check passed; uv venv py3.12; installed amm_sim_rs wheel, numpy, numba, cma, py-solc-x 2.1.0, pyrevm 0.3.7, and `-e .`; solc 0.8.24; `amm-match validate contracts/src/StarterStrategy.sol` → "validated successfully" |
| t03 | fetched api.imd.fun/artifacts/2333ecf2…; the file's sha256 is `2333ecf28d4356bd232b85c1e55842e94e1ada05e1cd067830b5625dfdbc5336`; `code_zb64` → base64 → zlib → 27,712-byte `contract Strategy`. The decoded source's sha256 is cb6d945511e1e10af78baa64c2a6cb25806f0998403ce6440654bfc3c70092d9, which matches the brief's cb6d9455 prefix. |
| `amm-match validate Strategy.sol` (the port) | **"validated successfully"**, but only with `ulimit -s 16384`. At the default 8 MB stack, solc 0.8.24 (viaIR, optimizer 200, paris; the settings in `compiler.py`) **segfaults (return code -11)** on this source, probably because of the 13.7 KB `hex"…"` literal. See Problems. |

Sandbox note: everything ran under `test/scratch/`, which is not delivered. The only deliverables are this file and sources.json, so the contract is given in full below rather than as `contracts/src/Strategy.sol`.

## Design

- **Filter.** A line-by-line port of `real.py` `RealPolicy.on_trade` with `{'soft': 1}` (the "plain:soft=1" gate variant) and `hard` = 0.995. All quantities are in natural-log space as signed WAD (1e18) values: m = mean of ln p; v = variance; ln(own spot) and ln(rival-spot estimate) are kept as logs. The fees are stored as L = −ln(1 − fee), so the no-arb interval is `lo = ln sp − Lb, hi = ln sp + La`, and the arb reading is `ln s − Lb` (isBuy) or `ln s + La`, with no extra ln or exp. The 10% fee clamp becomes `L ≤ −ln 0.9`, which is exactly the same.
- **Steps.** Same as README steps 1–6: the silent-step loop is capped at 12 (`min(dt−1,12)`), each step adding v += σ² and then truncating. Then the arb-vs-retail likelihood, the soft two-point mixture for 0.005 < π < 0.995, the band clip of the rival estimate against ln G = ln 0.997, the router update when not an arb, and phases 0, 1, 2.
- **Routines** (all fixed point, no assembly, since the validator rejects inline assembly):
  - `_ln`: msb range reduction to [1/√2, √2), then an atanh series of 7 terms.
  - `_exp`: k·ln2 reduction, then a 10-term Taylor series. Returns 0 below −42.
  - `_sqrt`: an msb seed, then 5 Newton steps.
  - Normal pdf: done in the log domain, `−x²/2 − ln√(2π)`.
  - Normal tail / cdf: `_lnQ(x)` = ln Q(x) through the rational-Chebyshev **erfc of Numerical Recipes (`erfcc`, fractional error < 1.2e-7 everywhere)**, kept in the log domain.
  - `_trunc`: the truncated-normal moments with the mass computed as ln Z. Z = Q(a) − Q(b) when a ≥ 0, Q(−b) − Q(−a) when b ≤ 0, else 1 − Q(−a) − Q(b). This avoids the tail cancellation that WAD precision would otherwise suffer. The Z < 1e-12 snap rule is `ln Z < ln 1e-12`, with snap to the nearer edge and v/4, as in real.py.
  - The arb/retail posterior is π = 1/(1 + exp(ln L_ret − ln L_arb)), so neither likelihood can underflow.
- **Errors vs the reference:** erfc 1.2e-7 relative (real.py uses `math.erf`), and ln/exp ~1e-13 relative. The result is a fee error of at most 6.6e-5 bps (Gate F).
- **Table access.** All reads go through one internal function, `_tv(blob, k)`: value k of a single `bytes private constant BLOB`, 12 bits per value, unit 0.25 bps, stored as u = 4·bps + 256. The header holds the knot lists as int16 bps, read by `_knot`. To swap in another table format, replace `_tv`, `_knot` and the layout constants that `gen.py` emits.
- **Interpolation.** `_seg` (binary search over the knots) and `_bil` do bilinear interpolation in (w, w2), clamped at the ends. The original lattice semantics are kept: rows are w in bps; phases 1 and 2 use the W2L columns ×5 bps (W2L is in 5 bps lattice units); phase 0 uses ±30 bps.
- **Storage.** `slots[0..7]`. Every value that can be zero is stored +1, and the timestamp, count and synced flag are packed into one word, so a slot never returns to zero. Without this, each zero→nonzero write cost about 22k gas.

## Table

- **Coarse table "A"**, built from `pol_mid.json` (plain, not qm). The knots are a subset of the original lattice, so values at the knots are exact before rounding.
  - Rows (own displacement w, bps; 57 knots): 0, ±5…±50 step 5, ±60…±100 step 10, ±120…±200 step 20, ±250…±400 step 50, ±500…±800 step 100.
  - Phase 1/2 columns (rival displacement, bps; 33 knots): 0, ±5, ±10, ±15, ±20, ±25, ±30, ±35, ±40, ±50, ±60, ±70, ±100, ±150, ±250, ±500, ±800 (W2L entries 1–8, 10, 12, 14, 20, 30, 50, 100, 160, ×5).
  - Phase 0 columns: all 13 (−30…30 step 5).
- **Size and coding.** 6 tables, 9,006 values, rounded to 0.25 bps, packed 12 bits each. The blob is 13,717 bytes (sha256 6c59b8b5ad5146e3d2d5a5d31abde17d2028c3d54384e61458c1bc18be6f5346).
- **Replica cost of the coarse table (indicative, not the formal Gate Q).** Coarse vs full, paired, real.py soft=1, replica seeds 60000..60255 (choose set):
  - table A at 0.1 bps (before the 0.25 bps rounding): **−0.103 se 0.019**, worst −1.00.
  - table B (85 rows × 35 columns, 28 KB at int16): −0.034 se 0.013, but too big for the bytecode.
  - Not run: the 12-bit A table exactly as shipped, and the holdout seeds 10000..12047.
- A first attempt read W2L as bps instead of 5 bps lattice units and lost 41 on the replica. It was fixed before anything shipped.

## Gate F

The Python reference is `ref.py`: real.py's `RealPolicy` unchanged except `_fees`, which is replaced by a bilinear lookup on exactly the shipped coarse table (same knots, same 0.25 bps values).

How it was recorded: replica runs (replica2.py) on seeds 70000..70049 were driven by the reference. Its first 2,000 `on_trade` callbacks per run were logged as (isBuy, amountX, amountY, timestamp, reserveX, reserveY) together with the reference's clamped (bid, ask) fees. Each run then deployed the compiled contract fresh in pyrevm, called `afterInitialize(100e18, 10000e18)` and replayed the 2,000 callbacks (WAD-encoded) through `afterSwap`. The contract's own state evolves; it is not reset to the reference.

| | value |
|---|---|
| callbacks compared | 100,050 (50 × 2,000 afterSwap + 50 afterInitialize) |
| mean abs fee diff (max of bid/ask) | **2.77e-7 bps** |
| median / p99.9 | 1.05e-7 / 1.19e-5 bps |
| **max** | **6.57e-5 bps** (seed 70033, callback 1319) |
| callbacks over 0.01 bps | 0 |

This replay used a 50M gas limit so the fee comparison could run; under the engine's limit the same replay fails (see Gas). A 1-run check before the gas trims gave a max of 2.7e-5 bps.

## Size

Deployed bytecode is **22,645 bytes** (limit 24,576), measured from `SolidityCompiler.compile(...).deployed_bytecode`. Creation and runtime pass the opcode scan.

## Gas

Engine `amm_sim_rs.gas_stats()`, 64 seeds 9510000..9510063, Submission role. gas_stats is cumulative within a process, so the t03 figures are its row minus the port's.

| | calls | ok | **OutOfGas** | max ok gas | mean ok gas |
|---|---|---|---|---|---|
| port afterSwap | 536,430 | 246,328 | **290,102** | 247,188 | 97,032 |
| port afterInitialize | 64 | 64 | 0 | 198,194 | |
| t03 afterSwap | 689,546 | 689,546 | 0 | (≤ 247,188; not separable) | 135,651 |

**Gas fails:** OutOfGas is not 0, and the 245,000 peak is not met.

- **What runs out.** With the engine limit of 250,000 in pyrevm, the first failing call in replay seed 70000 is callback 707, at dt = 10, which means 9 silent-step truncations. Each `_trunc` costs roughly 15–20k gas after the series trims.
- **Why it keeps failing.** The failed call is reverted, so the stored timestamp does not advance. Every later call then sees a larger dt, runs the full 12-step loop, and runs out of gas again. That accounts for the 54% failure rate.
- **History.** The first build had 23,329 failures in 24,254 calls. The +1 bias on the slots and shorter series cut that, but did not remove it.

## Paired vs t03

Each strategy ran alone against the 30 bps VanillaStrategy with `amm_sim_rs.run_batch(X, Vanilla, configs, 8)`. The configs are the `MatchRunner._build_configs` body (BASELINE_SETTINGS / BASELINE_VARIANCE) with i taken from the assigned seeds **9510000..9510063** (fresh, used once).

| | mean edge |
|---|---|
| port | 421.542 |
| t03 | 523.268 |
| **paired port − t03** | **n 64, mean_diff −101.726, se 7.696, worst5 −258.57, −205.97, −191.24, −182.92, −181.98** |

This measures a contract that runs out of gas on more than half its calls and then quotes stale fees. It says nothing about the policy itself, which on the replica is HA +2.83 against about +1.45 for t03, according to the brief.

## Problems

1. **Gas (blocking).** Ways to fit under 245k:
   - Collapse k silent steps analytically, or precompute: the cut prior after k silent steps depends only on (m−lo, hi−lo, v) through a short table.
   - Use cheaper erfc and exp (fewer terms; the error budget is about 1e-6, against today's ~1e-13).
   - Drop the per-call ~13.7 KB memory copy of the blob: it is a `bytes constant` and assembly is banned, so this is not measured yet.
   - Make the silent-step loop safe against stale timestamps. The cap stays 12, but an OutOfGas must never leave the timestamp unwritten.

   The cap of 12 is kept as the brief requires, so this has to be done by making each step cheaper.
2. **solc stack.** Compiling the port segfaults (rc −11) at an 8 MB stack and works at 16 MB. If the verifier compiles at the default stack, the contract will not build. A possible fix is to split the hex literal into several constants (untested).
3. Not done for lack of time: the formal Gate Q (paired on choose 60000..61023 and report 10000..12047 with the exact 12-bit table), a qm-table version, and a finer-row table (B scored −0.034 but needs a smaller encoding, for example 8-bit deltas).

## Code

The hex blob is elided from the Strategy.sol listing (27,434 hex characters) to fit the 64 KB limit. `python3 -I gen.py Strategy.sol` regenerates it byte for byte (sha256 of the full file: 63a8f38f1a126bee4753ba5bb4a3b0c2bc1899b1664e4737455d492ffb022ae4). gen.py imports ref.py, which reads `pol_mid.json` from the amm-sim checkout.

### Strategy.sol (blob elided)

```solidity
// SPDX-License-Identifier: MIT
pragma solidity ^0.8.24;

import {AMMStrategyBase} from "./AMMStrategyBase.sol";
import {TradeInfo} from "./IAMMStrategy.sol";

/// @title Gaussian-filter port of Taylor McKnight's solved AMM policy
/// @notice Policy, filter and tables are by Taylor (gtmcknight): github.com/gtmcknight/amm-sim,
///         commit 4812a6382311624e1506aca131ca35c378f930d9, research/solved (real.py, soft=1, plain tables pol_mid.json).
///         This file is a fixed-point port of real.py with a coarse copy of pol_mid.json:
///         rows (own displacement, bps) -800 -700 -600 -500 -400 -350 -300 -250 -200 -180 -160 -140 -120 -100 -90 -80 -70 -60 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 60 70 80 90 100 120 140 160 180 200 250 300 350 400 500 600 700 800;
///         phase 1/2 columns (rival displacement, bps) -800 -500 -250 -150 -100 -70 -60 -50 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 50 60 70 100 150 250 500 800; phase 0 columns -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30;
///         offsets rounded to 0.25 bps, packed 12 bits per value.
contract Strategy is AMMStrategyBase {
    int256 private constant W = 1e18;
    int256 private constant LN2 = 693147180559945309;
    int256 private constant LN_SQ2PI = 918938533204672742;
    int256 private constant SIG2 = 893025000000;      // (9.45 bps)^2
    int256 private constant LN_G = -3004509020298722;      // ln(0.997)
    int256 private constant LN_RHO = -1347073647966609323;  // ln(0.26)
    int256 private constant LD_M = -6801595083153746419;      // ln(11.12e-4)
    int256 private constant LD_S = 1176000000000000000;      // 1.176
    int256 private constant LN_LDS = 162118849476435178;
    int256 private constant LCAP = 105360515657826301;      // -ln(0.9): fee clamp at 10%
    int256 private constant LN_ZMIN = -27631021115928548208; // ln(1e-12)
    int256 private constant HARD = 995000000000000000;
    int256 private constant SQRT2 = 1414213562373095049;
    int256 private constant NEG = -1e30;

    uint256 private constant NR = 57;
    uint256 private constant NC1 = 33;
    uint256 private constant NC0 = 13;
    uint256 private constant OFF_C1 = 114;
    uint256 private constant OFF_C0 = 180;
    uint256 private constant T0 = 206;
    uint256 private constant N0 = 741;
    uint256 private constant N1 = 1881;

    // slots (stored +1 where a value can be 0, so a slot never returns to zero):
    //   0 m = mean of ln p, 1 v + 1, 2 Lb + 1 with Lb = -ln(1 - bid fee), 3 La + 1, 4 ln(own spot after last call),
    //   5 (last timestamp << 16) | (own retail count this step << 8) | (synced << 1) | 1, 7 ln(rival spot estimate)
    bytes private constant BLOB = hex"<13,717-byte blob elided: regenerate with gen.py>";

    function afterInitialize(uint256 initialX, uint256 initialY) external override returns (uint256, uint256) {
        int256 m = _ln(int256(initialY * 1e18 / initialX));
        slots[0] = uint256(m);
        slots[1] = 1;
        slots[4] = uint256(m);
        slots[5] = 1;
        slots[7] = uint256(m);
        return _quote(BLOB, 0, 0, 0);
    }

    function afterSwap(TradeInfo calldata trade) external override returns (uint256, uint256) {
        int256 ls = _ln(int256(trade.reserveY * 1e18 / trade.reserveX));
        int256 m = int256(slots[0]);
        int256 v = int256(slots[1]) - 1;
        int256 lb = int256(slots[2]) - 1;
        int256 la = int256(slots[3]) - 1;
        int256 lsp = int256(slots[4]);
        uint256 t = trade.timestamp;
        uint256 pk = slots[5];
        uint256 tl = pk >> 16;
        uint256 cnt = (pk >> 8) & 0xFF;
        int256 ls2 = int256(slots[7]);
        uint256 synced = (pk >> 1) & 1;
        bool isArb = false;
        if (t >= tl + 1) {
            cnt = 0;
            int256 lo = lsp - lb;
            int256 hi = lsp + la;
            uint256 n = t - tl - 1;
            if (n > 12) n = 12;
            int256 lz;
            for (uint256 k = 0; k < n; ++k) {
                v += SIG2;
                (lz, m, v) = _trunc(m, v, lo, hi);
            }
            v += SIG2;
            int256 lp = trade.isBuy ? ls - lb : ls + la;
            int256 sd = int256(_sqrt(uint256(v) * 1e18));
            int256 za = (lp - m) * W / sd;
            int256 lnArb = -(za * za / (2 * W)) - _ln(sd) - LN_SQ2PI;
            int256 mt;
            int256 vt;
            (lz, mt, vt) = _trunc(m, v, lo, hi);
            int256 d = ls - lsp;
            if (d < 0) d = -d;
            if (d < 1e9) d = 1e9;
            int256 lnd = _ln(d);
            int256 u = (lnd - LD_M) * W / LD_S;
            int256 lnRet = lz + LN_RHO - (u * u / (2 * W)) - lnd - LN_LDS - LN_SQ2PI;
            int256 dl = lnRet - lnArb;
            int256 pi;
            if (dl > 100 * W) pi = 0;
            else if (dl < -100 * W) pi = W;
            else pi = W * W / (W + _exp(dl));
            isArb = pi > W / 2;
            if (pi < HARD && pi > W - HARD) {
                int256 m2 = (pi * lp + (W - pi) * mt) / W;
                int256 e1 = lp - m2;
                int256 e2 = mt - m2;
                v = (pi * (e1 * e1 / W) + (W - pi) * (vt + e2 * e2 / W)) / W;
                m = m2;
            } else if (isArb) {
                m = lp;
                v = 0;
            } else {
                m = mt;
                v = vt;
            }
            if (ls2 < m + LN_G) ls2 = m + LN_G;
            else if (ls2 > m - LN_G) ls2 = m - LN_G;
            synced = 0;
        }
        uint256 phs = 0;
        if (!isArb) {
            if (trade.isBuy) {
                int256 cand = ls - lb - LN_G;
                if (cand < ls2) { ls2 = cand; synced = 1; }
            } else {
                int256 cand = ls + LN_G + la;
                if (cand > ls2) { ls2 = cand; synced = 1; }
            }
            if (cnt < 255) cnt += 1;
            phs = cnt <= 1 ? 1 : 2;
        }
        slots[0] = uint256(m);
        slots[1] = uint256(v + 1);
        slots[4] = uint256(ls);
        slots[5] = (t << 16) | (cnt << 8) | (synced << 1) | 1;
        slots[7] = uint256(ls2);
        return _quote(BLOB, phs, (ls - m) * 1e4, (ls2 - m) * 1e4);
    }

    /// fees from the table at own displacement w and rival displacement w2 (bps, WAD); stores La, Lb
    function _quote(bytes memory blob, uint256 phs, int256 w, int256 w2) internal returns (uint256 bidFee, uint256 askFee) {
        (uint256 i, int256 f) = _seg(blob, 0, NR, w);
        uint256 j;
        int256 g;
        uint256 nc;
        uint256 base;
        if (phs == 0) {
            (j, g) = _seg(blob, OFF_C0, NC0, w2);
            nc = NC0;
            base = 0;
        } else {
            (j, g) = _seg(blob, OFF_C1, NC1, w2);
            nc = NC1;
            base = 2 * N0 + (phs - 1) * 2 * N1;
        }
        uint256 tabSize = phs == 0 ? N0 : N1;
        int256 a = _bil(blob, base + i * nc + j, nc, f, g);
        int256 b = _bil(blob, base + tabSize + i * nc + j, nc, f, g);
        int256 la = 0;
        int256 lb = 0;
        if (a > w) { la = (a - w) / 1e4; if (la > LCAP) la = LCAP; }
        if (b > -w) { lb = (b + w) / 1e4; if (lb > LCAP) lb = LCAP; }
        slots[2] = uint256(lb + 1);
        slots[3] = uint256(la + 1);
        bidFee = uint256(W - _exp(-lb));
        askFee = uint256(W - _exp(-la));
    }

    function _bil(bytes memory blob, uint256 k, uint256 nc, int256 f, int256 g) internal pure returns (int256) {
        unchecked {
            int256 v00 = _tv(blob, k);
            int256 v10 = _tv(blob, k + nc);
            int256 v01 = _tv(blob, k + 1);
            int256 v11 = _tv(blob, k + nc + 1);
            int256 c0 = v00 * (W - f) + v10 * f;
            int256 c1 = v01 * (W - f) + v11 * f;
            return (c0 * (W - g) + c1 * g) / W / 4;   // 0.25 bps units -> bps WAD
        }
    }

    /// the one table read: value k of the packed table, in 0.25 bps units
    function _tv(bytes memory blob, uint256 k) internal pure returns (int256 r) {
        unchecked {
            uint256 off = T0 + (k * 3) / 2;
            uint256 x = (uint256(uint8(blob[off])) << 8) | uint256(uint8(blob[off + 1]));
            if (k & 1 == 0) r = int256(x >> 4);
            else r = int256(x & 0xFFF);
            r -= 256;
        }
    }

    function _knot(bytes memory blob, uint256 off, uint256 i) internal pure returns (int256 r) {
        unchecked {
            uint256 p = off + 2 * i;
            uint256 x = (uint256(uint8(blob[p])) << 8) | uint256(uint8(blob[p + 1]));
            r = int256(int16(uint16(x))) * W;
        }
    }

    /// segment and fraction of x (bps, WAD) on a knot list; clamps at the ends
    function _seg(bytes memory blob, uint256 off, uint256 n, int256 x) internal pure returns (uint256, int256) {
        unchecked {
            if (x <= _knot(blob, off, 0)) return (0, 0);
            if (x >= _knot(blob, off, n - 1)) return (n - 2, W);
            uint256 lo = 0;
            uint256 hi = n - 1;
            while (hi - lo > 1) {
                uint256 mid = (lo + hi) >> 1;
                if (_knot(blob, off, mid) <= x) lo = mid;
                else hi = mid;
            }
            int256 kl = _knot(blob, off, lo);
            return (lo, (x - kl) * W / (_knot(blob, off, hi) - kl));
        }
    }

    /// moments of N(m, v) restricted to (lo, hi): (ln mass, mean, var)
    function _trunc(int256 m, int256 v, int256 lo, int256 hi) internal pure returns (int256 lz, int256 mean, int256 var_) {
        unchecked {
            if (v <= 0) return ((lo <= m && m <= hi) ? int256(0) : NEG, m, 0);
            int256 sd = int256(_sqrt(uint256(v) * 1e18));
            int256 a = (lo - m) * W / sd;
            int256 b = (hi - m) * W / sd;
            if (a >= 0) {
                lz = _lnDiff(_lnQ(a), _lnQ(b));
            } else if (b <= 0) {
                lz = _lnDiff(_lnQ(-b), _lnQ(-a));
            } else {
                int256 z = W - _exp(_lnQ(-a)) - _exp(_lnQ(b));
                lz = z > 0 ? _ln(z) : NEG;
            }
            if (lz < LN_ZMIN) return (lz, m < lo ? lo : hi, v / 4);
            int256 pa = _exp(-(a * a / (2 * W)) - LN_SQ2PI - lz);
            int256 pb = _exp(-(b * b / (2 * W)) - LN_SQ2PI - lz);
            int256 r = pa - pb;
            mean = m + sd * r / W;
            var_ = v * (W + (a * pa - b * pb) / W - r * r / W) / W;
            if (var_ < 0) var_ = 0;
        }
    }

    /// ln(e^x - e^y) for x >= y
    function _lnDiff(int256 x, int256 y) internal pure returns (int256) {
        unchecked {
            int256 e = _exp(y - x);
            if (e >= W) return NEG;
            return x + _ln(W - e);
        }
    }

    /// ln of the upper normal tail Q(x), x >= 0, via the erfc of Numerical Recipes (fractional error < 1.2e-7)
    function _lnQ(int256 x) internal pure returns (int256) {
        unchecked {
            int256 z = x * W / SQRT2;
            int256 t = W * W / (W + z / 2);
            int256 p = 170872770000000000;
            p = -822152230000000000 + p * t / W;
            p = 1488515870000000000 + p * t / W;
            p = -1135203980000000000 + p * t / W;
            p = 278868070000000000 + p * t / W;
            p = -186288060000000000 + p * t / W;
            p = 96784180000000000 + p * t / W;
            p = 374091960000000000 + p * t / W;
            p = 1000023680000000000 + p * t / W;
            p = -1265512230000000000 + p * t / W;
            return -LN2 + _ln(t) - z * z / W + p;
        }
    }

    function _ln(int256 x) internal pure returns (int256) {
        unchecked {
            uint256 ux = uint256(x);
            int256 e = int256(_msb(ux)) - 59;
            uint256 y = e >= 0 ? ux >> uint256(e) : ux << uint256(-e);
            if (y < 707106781186547524) { y <<= 1; e -= 1; }
            int256 s = (int256(y) - W) * W / (int256(y) + W);
            int256 s2 = s * s / W;
            int256 term = s;
            int256 sum = s;
            for (int256 i = 3; i <= 15; i += 2) {
                term = term * s2 / W;
                sum += term / i;
            }
            return e * LN2 + 2 * sum;
        }
    }

    function _exp(int256 x) internal pure returns (int256) {
        unchecked {
            if (x < -42 * W) return 0;
            require(x < 130 * W, "exp");
            int256 k = x >= 0 ? (x + LN2 / 2) / LN2 : (x - LN2 / 2) / LN2;
            int256 r = x - k * LN2;
            int256 term = W;
            int256 sum = W;
            for (int256 i = 1; i <= 10; ++i) {
                term = term * r / (W * i);
                sum += term;
            }
            if (k >= 0) return sum << uint256(k);
            return sum >> uint256(-k);
        }
    }

    function _msb(uint256 x) internal pure returns (uint256 r) {
        if (x >= 1 << 128) { x >>= 128; r += 128; }
        if (x >= 1 << 64) { x >>= 64; r += 64; }
        if (x >= 1 << 32) { x >>= 32; r += 32; }
        if (x >= 1 << 16) { x >>= 16; r += 16; }
        if (x >= 1 << 8) { x >>= 8; r += 8; }
        if (x >= 1 << 4) { x >>= 4; r += 4; }
        if (x >= 1 << 2) { x >>= 2; r += 2; }
        if (x >= 1 << 1) { r += 1; }
    }

    function _sqrt(uint256 x) internal pure returns (uint256 y) {
        unchecked {
            if (x == 0) return 0;
            y = 1 << ((_msb(x) + 1) / 2);
            for (uint256 i = 0; i < 5; ++i) {
                y = (y + x / y) >> 1;
            }
            uint256 z = x / y;
            if (z < y) y = z;
        }
    }

    function getName() external pure override returns (string memory) {
        return "Taylor solved policy (Gaussian+plain coarse)";
    }
}

```

### gen.py (table packer + contract template)

```python
"""Build Strategy.sol: coarse table (ref.py knots) packed 12 bits per value into one bytes constant.
usage: python3 -I gen.py OUT.sol"""
import sys, os
from decimal import Decimal, getcontext
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import ref

getcontext().prec = 60
D = Decimal
W = D(10) ** 18
def wad(x): return int((D(x) * W).to_integral_value())
def ln(x): return D(x).ln()

T = ref.T_full(); C = ref.coarse(T)
RK, CK1, CK0 = ref.RK, ref.CK1, ref.CK0
NR, NC1, NC0 = len(RK), len(CK1), len(CK0)
blob = bytearray()
for k in RK + CK1 + CK0:
    blob += (k & 0xFFFF).to_bytes(2, 'big')
T0 = len(blob)
vals = []
for name in ref.NAMES:
    for row in C[name]:
        for v in row:
            u = v + 256                       # v in 0.25 bps; u = 4*(bps + 64)
            assert 0 <= u < 4096, (name, v)
            vals.append(u)
if len(vals) % 2: vals.append(0)
for k in range(0, len(vals), 2):
    x = (vals[k] << 12) | vals[k + 1]
    blob += x.to_bytes(3, 'big')
blob += b'\x00\x00'                           # padding so a 2-byte read never runs past the end
N0 = NR * NC0; N1 = NR * NC1

const = dict(
    LN2=wad(ln(2)), LN_SQ2PI=wad(ln(D(2) * D('3.14159265358979323846264338327950288419716939937510')) / 2),
    SIG2=wad(D('0.000945') ** 2), LN_G=wad(ln(D('0.997'))), LN_RHO=wad(ln(D('0.26'))),
    LD_M=wad(ln(D('11.12e-4'))), LD_S=wad('1.176'), LN_LDS=wad(ln(D('1.176'))),
    LCAP=wad(-ln(D('0.9'))), LN_ZMIN=wad(ln(D('1e-12'))), HARD=wad('0.995'),
    SQRT2=wad(D(2).sqrt()),
)

SRC = r'''// SPDX-License-Identifier: MIT
pragma solidity ^0.8.24;

import {AMMStrategyBase} from "./AMMStrategyBase.sol";
import {TradeInfo} from "./IAMMStrategy.sol";

/// @title Gaussian-filter port of Taylor McKnight's solved AMM policy
/// @notice Policy, filter and tables are by Taylor (gtmcknight): github.com/gtmcknight/amm-sim,
///         commit 4812a6382311624e1506aca131ca35c378f930d9, research/solved (real.py, soft=1, plain tables pol_mid.json).
///         This file is a fixed-point port of real.py with a coarse copy of pol_mid.json:
///         rows (own displacement, bps) %%RK%%;
///         phase 1/2 columns (rival displacement, bps) %%CK1%%; phase 0 columns %%CK0%%;
///         offsets rounded to 0.25 bps, packed 12 bits per value.
contract Strategy is AMMStrategyBase {
    int256 private constant W = 1e18;
    int256 private constant LN2 = %%LN2%%;
    int256 private constant LN_SQ2PI = %%LN_SQ2PI%%;
    int256 private constant SIG2 = %%SIG2%%;      // (9.45 bps)^2
    int256 private constant LN_G = %%LN_G%%;      // ln(0.997)
    int256 private constant LN_RHO = %%LN_RHO%%;  // ln(0.26)
    int256 private constant LD_M = %%LD_M%%;      // ln(11.12e-4)
    int256 private constant LD_S = %%LD_S%%;      // 1.176
    int256 private constant LN_LDS = %%LN_LDS%%;
    int256 private constant LCAP = %%LCAP%%;      // -ln(0.9): fee clamp at 10%
    int256 private constant LN_ZMIN = %%LN_ZMIN%%; // ln(1e-12)
    int256 private constant HARD = %%HARD%%;
    int256 private constant SQRT2 = %%SQRT2%%;
    int256 private constant NEG = -1e30;

    uint256 private constant NR = %%NR%%;
    uint256 private constant NC1 = %%NC1%%;
    uint256 private constant NC0 = %%NC0%%;
    uint256 private constant OFF_C1 = %%OFF_C1%%;
    uint256 private constant OFF_C0 = %%OFF_C0%%;
    uint256 private constant T0 = %%T0%%;
    uint256 private constant N0 = %%N0%%;
    uint256 private constant N1 = %%N1%%;

    // slots (stored +1 where a value can be 0, so a slot never returns to zero):
    //   0 m = mean of ln p, 1 v + 1, 2 Lb + 1 with Lb = -ln(1 - bid fee), 3 La + 1, 4 ln(own spot after last call),
    //   5 (last timestamp << 16) | (own retail count this step << 8) | (synced << 1) | 1, 7 ln(rival spot estimate)
    bytes private constant BLOB = hex"%%BLOB%%";

    function afterInitialize(uint256 initialX, uint256 initialY) external override returns (uint256, uint256) {
        int256 m = _ln(int256(initialY * 1e18 / initialX));
        slots[0] = uint256(m);
        slots[1] = 1;
        slots[4] = uint256(m);
        slots[5] = 1;
        slots[7] = uint256(m);
        return _quote(BLOB, 0, 0, 0);
    }

    function afterSwap(TradeInfo calldata trade) external override returns (uint256, uint256) {
        int256 ls = _ln(int256(trade.reserveY * 1e18 / trade.reserveX));
        int256 m = int256(slots[0]);
        int256 v = int256(slots[1]) - 1;
        int256 lb = int256(slots[2]) - 1;
        int256 la = int256(slots[3]) - 1;
        int256 lsp = int256(slots[4]);
        uint256 t = trade.timestamp;
        uint256 pk = slots[5];
        uint256 tl = pk >> 16;
        uint256 cnt = (pk >> 8) & 0xFF;
        int256 ls2 = int256(slots[7]);
        uint256 synced = (pk >> 1) & 1;
        bool isArb = false;
        if (t >= tl + 1) {
            cnt = 0;
            int256 lo = lsp - lb;
            int256 hi = lsp + la;
            uint256 n = t - tl - 1;
            if (n > 12) n = 12;
            int256 lz;
            for (uint256 k = 0; k < n; ++k) {
                v += SIG2;
                (lz, m, v) = _trunc(m, v, lo, hi);
            }
            v += SIG2;
            int256 lp = trade.isBuy ? ls - lb : ls + la;
            int256 sd = int256(_sqrt(uint256(v) * 1e18));
            int256 za = (lp - m) * W / sd;
            int256 lnArb = -(za * za / (2 * W)) - _ln(sd) - LN_SQ2PI;
            int256 mt;
            int256 vt;
            (lz, mt, vt) = _trunc(m, v, lo, hi);
            int256 d = ls - lsp;
            if (d < 0) d = -d;
            if (d < 1e9) d = 1e9;
            int256 lnd = _ln(d);
            int256 u = (lnd - LD_M) * W / LD_S;
            int256 lnRet = lz + LN_RHO - (u * u / (2 * W)) - lnd - LN_LDS - LN_SQ2PI;
            int256 dl = lnRet - lnArb;
            int256 pi;
            if (dl > 100 * W) pi = 0;
            else if (dl < -100 * W) pi = W;
            else pi = W * W / (W + _exp(dl));
            isArb = pi > W / 2;
            if (pi < HARD && pi > W - HARD) {
                int256 m2 = (pi * lp + (W - pi) * mt) / W;
                int256 e1 = lp - m2;
                int256 e2 = mt - m2;
                v = (pi * (e1 * e1 / W) + (W - pi) * (vt + e2 * e2 / W)) / W;
                m = m2;
            } else if (isArb) {
                m = lp;
                v = 0;
            } else {
                m = mt;
                v = vt;
            }
            if (ls2 < m + LN_G) ls2 = m + LN_G;
            else if (ls2 > m - LN_G) ls2 = m - LN_G;
            synced = 0;
        }
        uint256 phs = 0;
        if (!isArb) {
            if (trade.isBuy) {
                int256 cand = ls - lb - LN_G;
                if (cand < ls2) { ls2 = cand; synced = 1; }
            } else {
                int256 cand = ls + LN_G + la;
                if (cand > ls2) { ls2 = cand; synced = 1; }
            }
            if (cnt < 255) cnt += 1;
            phs = cnt <= 1 ? 1 : 2;
        }
        slots[0] = uint256(m);
        slots[1] = uint256(v + 1);
        slots[4] = uint256(ls);
        slots[5] = (t << 16) | (cnt << 8) | (synced << 1) | 1;
        slots[7] = uint256(ls2);
        return _quote(BLOB, phs, (ls - m) * 1e4, (ls2 - m) * 1e4);
    }

    /// fees from the table at own displacement w and rival displacement w2 (bps, WAD); stores La, Lb
    function _quote(bytes memory blob, uint256 phs, int256 w, int256 w2) internal returns (uint256 bidFee, uint256 askFee) {
        (uint256 i, int256 f) = _seg(blob, 0, NR, w);
        uint256 j;
        int256 g;
        uint256 nc;
        uint256 base;
        if (phs == 0) {
            (j, g) = _seg(blob, OFF_C0, NC0, w2);
            nc = NC0;
            base = 0;
        } else {
            (j, g) = _seg(blob, OFF_C1, NC1, w2);
            nc = NC1;
            base = 2 * N0 + (phs - 1) * 2 * N1;
        }
        uint256 tabSize = phs == 0 ? N0 : N1;
        int256 a = _bil(blob, base + i * nc + j, nc, f, g);
        int256 b = _bil(blob, base + tabSize + i * nc + j, nc, f, g);
        int256 la = 0;
        int256 lb = 0;
        if (a > w) { la = (a - w) / 1e4; if (la > LCAP) la = LCAP; }
        if (b > -w) { lb = (b + w) / 1e4; if (lb > LCAP) lb = LCAP; }
        slots[2] = uint256(lb + 1);
        slots[3] = uint256(la + 1);
        bidFee = uint256(W - _exp(-lb));
        askFee = uint256(W - _exp(-la));
    }

    function _bil(bytes memory blob, uint256 k, uint256 nc, int256 f, int256 g) internal pure returns (int256) {
        unchecked {
            int256 v00 = _tv(blob, k);
            int256 v10 = _tv(blob, k + nc);
            int256 v01 = _tv(blob, k + 1);
            int256 v11 = _tv(blob, k + nc + 1);
            int256 c0 = v00 * (W - f) + v10 * f;
            int256 c1 = v01 * (W - f) + v11 * f;
            return (c0 * (W - g) + c1 * g) / W / 4;   // 0.25 bps units -> bps WAD
        }
    }

    /// the one table read: value k of the packed table, in 0.25 bps units
    function _tv(bytes memory blob, uint256 k) internal pure returns (int256 r) {
        unchecked {
            uint256 off = T0 + (k * 3) / 2;
            uint256 x = (uint256(uint8(blob[off])) << 8) | uint256(uint8(blob[off + 1]));
            if (k & 1 == 0) r = int256(x >> 4);
            else r = int256(x & 0xFFF);
            r -= 256;
        }
    }

    function _knot(bytes memory blob, uint256 off, uint256 i) internal pure returns (int256 r) {
        unchecked {
            uint256 p = off + 2 * i;
            uint256 x = (uint256(uint8(blob[p])) << 8) | uint256(uint8(blob[p + 1]));
            r = int256(int16(uint16(x))) * W;
        }
    }

    /// segment and fraction of x (bps, WAD) on a knot list; clamps at the ends
    function _seg(bytes memory blob, uint256 off, uint256 n, int256 x) internal pure returns (uint256, int256) {
        unchecked {
            if (x <= _knot(blob, off, 0)) return (0, 0);
            if (x >= _knot(blob, off, n - 1)) return (n - 2, W);
            uint256 lo = 0;
            uint256 hi = n - 1;
            while (hi - lo > 1) {
                uint256 mid = (lo + hi) >> 1;
                if (_knot(blob, off, mid) <= x) lo = mid;
                else hi = mid;
            }
            int256 kl = _knot(blob, off, lo);
            return (lo, (x - kl) * W / (_knot(blob, off, hi) - kl));
        }
    }

    /// moments of N(m, v) restricted to (lo, hi): (ln mass, mean, var)
    function _trunc(int256 m, int256 v, int256 lo, int256 hi) internal pure returns (int256 lz, int256 mean, int256 var_) {
        unchecked {
            if (v <= 0) return ((lo <= m && m <= hi) ? int256(0) : NEG, m, 0);
            int256 sd = int256(_sqrt(uint256(v) * 1e18));
            int256 a = (lo - m) * W / sd;
            int256 b = (hi - m) * W / sd;
            if (a >= 0) {
                lz = _lnDiff(_lnQ(a), _lnQ(b));
            } else if (b <= 0) {
                lz = _lnDiff(_lnQ(-b), _lnQ(-a));
            } else {
                int256 z = W - _exp(_lnQ(-a)) - _exp(_lnQ(b));
                lz = z > 0 ? _ln(z) : NEG;
            }
            if (lz < LN_ZMIN) return (lz, m < lo ? lo : hi, v / 4);
            int256 pa = _exp(-(a * a / (2 * W)) - LN_SQ2PI - lz);
            int256 pb = _exp(-(b * b / (2 * W)) - LN_SQ2PI - lz);
            int256 r = pa - pb;
            mean = m + sd * r / W;
            var_ = v * (W + (a * pa - b * pb) / W - r * r / W) / W;
            if (var_ < 0) var_ = 0;
        }
    }

    /// ln(e^x - e^y) for x >= y
    function _lnDiff(int256 x, int256 y) internal pure returns (int256) {
        unchecked {
            int256 e = _exp(y - x);
            if (e >= W) return NEG;
            return x + _ln(W - e);
        }
    }

    /// ln of the upper normal tail Q(x), x >= 0, via the erfc of Numerical Recipes (fractional error < 1.2e-7)
    function _lnQ(int256 x) internal pure returns (int256) {
        unchecked {
            int256 z = x * W / SQRT2;
            int256 t = W * W / (W + z / 2);
            int256 p = 170872770000000000;
            p = -822152230000000000 + p * t / W;
            p = 1488515870000000000 + p * t / W;
            p = -1135203980000000000 + p * t / W;
            p = 278868070000000000 + p * t / W;
            p = -186288060000000000 + p * t / W;
            p = 96784180000000000 + p * t / W;
            p = 374091960000000000 + p * t / W;
            p = 1000023680000000000 + p * t / W;
            p = -1265512230000000000 + p * t / W;
            return -LN2 + _ln(t) - z * z / W + p;
        }
    }

    function _ln(int256 x) internal pure returns (int256) {
        unchecked {
            uint256 ux = uint256(x);
            int256 e = int256(_msb(ux)) - 59;
            uint256 y = e >= 0 ? ux >> uint256(e) : ux << uint256(-e);
            if (y < 707106781186547524) { y <<= 1; e -= 1; }
            int256 s = (int256(y) - W) * W / (int256(y) + W);
            int256 s2 = s * s / W;
            int256 term = s;
            int256 sum = s;
            for (int256 i = 3; i <= 15; i += 2) {
                term = term * s2 / W;
                sum += term / i;
            }
            return e * LN2 + 2 * sum;
        }
    }

    function _exp(int256 x) internal pure returns (int256) {
        unchecked {
            if (x < -42 * W) return 0;
            require(x < 130 * W, "exp");
            int256 k = x >= 0 ? (x + LN2 / 2) / LN2 : (x - LN2 / 2) / LN2;
            int256 r = x - k * LN2;
            int256 term = W;
            int256 sum = W;
            for (int256 i = 1; i <= 10; ++i) {
                term = term * r / (W * i);
                sum += term;
            }
            if (k >= 0) return sum << uint256(k);
            return sum >> uint256(-k);
        }
    }

    function _msb(uint256 x) internal pure returns (uint256 r) {
        if (x >= 1 << 128) { x >>= 128; r += 128; }
        if (x >= 1 << 64) { x >>= 64; r += 64; }
        if (x >= 1 << 32) { x >>= 32; r += 32; }
        if (x >= 1 << 16) { x >>= 16; r += 16; }
        if (x >= 1 << 8) { x >>= 8; r += 8; }
        if (x >= 1 << 4) { x >>= 4; r += 4; }
        if (x >= 1 << 2) { x >>= 2; r += 2; }
        if (x >= 1 << 1) { r += 1; }
    }

    function _sqrt(uint256 x) internal pure returns (uint256 y) {
        unchecked {
            if (x == 0) return 0;
            y = 1 << ((_msb(x) + 1) / 2);
            for (uint256 i = 0; i < 5; ++i) {
                y = (y + x / y) >> 1;
            }
            uint256 z = x / y;
            if (z < y) y = z;
        }
    }

    function getName() external pure override returns (string memory) {
        return "Taylor solved policy (Gaussian+plain coarse)";
    }
}
'''

rep = dict(const)
rep.update(NR=NR, NC1=NC1, NC0=NC0, OFF_C1=2 * NR, OFF_C0=2 * (NR + NC1), T0=T0, N0=N0, N1=N1,
           BLOB=blob.hex(), RK=' '.join(map(str, RK)), CK1=' '.join(map(str, CK1)), CK0=' '.join(map(str, CK0)))
out = SRC
for k, v in rep.items():
    out = out.replace('%%' + k + '%%', str(v))
assert '%%' not in out
open(sys.argv[1], 'w').write(out)
print('blob bytes', len(blob), 'values', len(vals), 'source bytes', len(out))

```

### ref.py (Python reference)

```python
"""Python reference: Taylor's real.py (Gaussian filter, soft=1) with the coarse table that the contract ships.
Only the table lookup (_fees) is replaced; the filter is real.py unchanged."""
import sys, os, math, json
HERE = os.path.dirname(os.path.abspath(__file__))
SOL = os.path.join(HERE, 'amm-sim', 'research', 'solved')
sys.path.insert(0, SOL)
import replica2 as R, real
from evalpol import load

# coarse knots (bps). Rows/cols are subsets of the original lattice, so knot values are exact (then rounded to 0.1 bps).
_V = os.environ.get('KNOTS', 'A')
if _V == 'A':
    _side = list(range(5, 51, 5)) + list(range(60, 101, 10)) + list(range(120, 201, 20)) + list(range(250, 401, 50)) + list(range(500, 801, 100))
    _c = [1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 20, 30, 50, 100, 160]
else:
    _side = list(range(5, 101, 5)) + list(range(110, 201, 10)) + list(range(225, 401, 25)) + list(range(500, 801, 100))
    _c = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 20, 30, 50, 100, 160]
RK = [-x for x in reversed(_side)] + [0] + _side
CK1 = [-5 * x for x in reversed(_c)] + [0] + [5 * x for x in _c]          # bps; W2L is in 5 bps lattice units
CK0 = [-30, -25, -20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30]
NAMES = ['A0', 'B0', 'A1', 'B1', 'A2', 'B2']


def coarse(T):
    W2L = T['W2L']; out = {}
    for k in NAMES:
        ck = CK0 if k.endswith('0') else CK1
        cols = [(c // 5 + T['FB']) if k.endswith('0') else W2L.index(float(c // 5)) for c in ck]
        rows = [r // 5 + T['W0'] for r in RK]
        out[k] = [[int(round(T[k][r][c] * 4.0)) for c in cols] for r in rows]   # 0.25 bps units (12-bit packed as u = v + 256)
    return out


def _seg(K, v):
    n = len(K)
    if v <= K[0]: return 0, 0.0
    if v >= K[n - 1]: return n - 2, 1.0
    lo, hi = 0, n - 1
    while hi - lo > 1:
        mid = (lo + hi) // 2
        if K[mid] <= v: lo = mid
        else: hi = mid
    return lo, (v - K[lo]) / (K[hi] - K[lo])


class RefPolicy(real.RealPolicy):
    def __init__(self, T, opt=None, C=None):
        self.C = C if C is not None else coarse(T)
        real.RealPolicy.__init__(self, T, opt if opt is not None else {'soft': 1.0})
    def _fees(self, w, w2, ph):
        C = self.C
        i, f = _seg(RK, w)
        if ph == 0:
            j, g = _seg(CK0, w2)
        else:
            j, g = _seg(CK1, w2)
        A, B = C['A%d' % ph], C['B%d' % ph]
        a = ((A[i][j] * (1 - f) + A[i + 1][j] * f) * (1 - g) + (A[i][j + 1] * (1 - f) + A[i + 1][j + 1] * f) * g) * 0.25
        b = ((B[i][j] * (1 - f) + B[i + 1][j] * f) * (1 - g) + (B[i][j + 1] * (1 - f) + B[i + 1][j + 1] * f) * g) * 0.25
        fa = 1.0 - math.exp(-(a - w) * 1e-4) if a > w else 0.0
        fb = 1.0 - math.exp(-(b + w) * 1e-4) if b > -w else 0.0
        return fb, fa


def T_full():
    return load(os.path.join(SOL, 'pol_mid.json'))

```

### gatef.py (Gate F harness)

```python
"""Gate F: record the first 2,000 callbacks of replica runs driven by the Python reference (ref.py), replay the same
callbacks through the compiled contract in pyrevm, compare fees.
usage: python gatef.py Strategy.sol SEED0 NRUNS"""
import sys, os
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import ref
from pyrevm import EVM
from amm_competition.evm.compiler import SolidityCompiler
from amm_competition.evm.executor import EVMStrategyExecutor as X

NCB = 2000
GAS = int(os.environ.get('GAS', X.GAS_LIMIT_TRADE))


class Stop(Exception):
    pass


class Rec(ref.RefPolicy):
    def __init__(self, T, C):
        ref.RefPolicy.__init__(self, T, {'soft': 1.0}, C)
        self.log = []
    def init(self, *a):
        fb, fa = ref.RefPolicy.init(self, *a)
        self.f0 = (self.fb, self.fa)
        return fb, fa
    def on_trade(self, is_buy, ax, ay, t, x, y, arb, p, s2, j):
        fb, fa = ref.RefPolicy.on_trade(self, is_buy, ax, ay, t, x, y, arb, p, s2, j)
        self.log.append((is_buy, ax, ay, t, x, y, self.fb, self.fa))
        if len(self.log) >= NCB:
            raise Stop()
        return fb, fa


def w32(v): return int(v).to_bytes(32, 'big')


def main():
    src = open(sys.argv[1]).read(); s0 = int(sys.argv[2]); n = int(sys.argv[3])
    res = SolidityCompiler().compile(src)
    assert res.success, res.errors
    T = ref.T_full(); C = ref.coarse(T)
    diffs = []; worst = (0, None)
    for seed in range(s0, s0 + n):
        pol = Rec(T, C)
        try:
            ref.R.run(seed, pol)
        except Stop:
            pass
        evm = EVM(); addr = evm.deploy(deployer=X.CALLER_ADDRESS, code=res.bytecode, value=0, gas=X.GAS_LIMIT_DEPLOY)
        out = evm.message_call(caller=X.CALLER_ADDRESS, to=addr, calldata=X.SELECTOR_AFTER_INITIALIZE + w32(100 * 10**18) + w32(10000 * 10**18), value=0, gas=X.GAS_LIMIT_INIT)
        fb0 = int.from_bytes(out[:32], 'big') / 1e18; fa0 = int.from_bytes(out[32:64], 'big') / 1e18
        diffs.append(max(abs(fb0 - pol.f0[0]), abs(fa0 - pol.f0[1])) * 1e4)
        for k, (b, ax, ay, t, x, y, rfb, rfa) in enumerate(pol.log):
            cd = X.SELECTOR_AFTER_SWAP + w32(1 if b else 0) + w32(ax * 1e18) + w32(ay * 1e18) + w32(t) + w32(x * 1e18) + w32(y * 1e18)
            try:
                out = evm.message_call(caller=X.CALLER_ADDRESS, to=addr, calldata=cd, value=0, gas=GAS)
            except Exception as e:
                print('FAIL', seed, k, (b, ax, ay, t, x, y), pol.log[k-1][3] if k else None, e); raise
            cfb = int.from_bytes(out[:32], 'big') / 1e18; cfa = int.from_bytes(out[32:64], 'big') / 1e18
            d = max(abs(cfb - rfb), abs(cfa - rfa)) * 1e4
            diffs.append(d)
            if d > worst[0]: worst = (d, (seed, k, b, t, rfb * 1e4, rfa * 1e4, cfb * 1e4, cfa * 1e4))
    diffs.sort()
    m = sum(diffs) / len(diffs)
    print(f'runs {n} seeds {s0}..{s0+n-1} callbacks {len(diffs)} (incl. {n} afterInitialize)')
    print(f'fee diff bps: mean {m:.3e} median {diffs[len(diffs)//2]:.3e} p99.9 {diffs[int(.999*(len(diffs)-1))]:.3e} max {diffs[-1]:.3e}; >0.01 bps: {sum(1 for d in diffs if d > 0.01)}')
    print('worst (seed, idx, isBuy, t, ref bid, ref ask, sol bid, sol ask):', worst)


main()

```

### engine.py (engine harness)

```python
"""Engine: each strategy alone vs the 30 bps VanillaStrategy on fresh seeds, configs per MatchRunner._build_configs.
usage: python engine.py OUT.json SEED0 N file.sol [file.sol ...]"""
import sys, os, json, time
import amm_sim_rs
from amm_competition.competition.match import MatchRunner
from amm_competition.competition.config import BASELINE_SETTINGS, BASELINE_VARIANCE
from amm_competition.evm.compiler import SolidityCompiler
from amm_competition.evm.baseline import get_vanilla_bytecode_and_abi


def configs(seeds):
    """MatchRunner._build_configs body with seed i taken from `seeds` (baseline settings and variance)."""
    import numpy as np
    B, V = BASELINE_SETTINGS, BASELINE_VARIANCE
    out = []
    for i in seeds:
        rng = np.random.default_rng(seed=i)
        size = rng.uniform(V.retail_mean_size_min, V.retail_mean_size_max)
        rate = rng.uniform(V.retail_arrival_rate_min, V.retail_arrival_rate_max)
        sig = rng.uniform(V.gbm_sigma_min, V.gbm_sigma_max)
        out.append(amm_sim_rs.SimulationConfig(n_steps=B.n_steps, initial_price=B.initial_price, initial_x=B.initial_x,
                   initial_y=B.initial_y, gbm_mu=B.gbm_mu, gbm_sigma=sig, gbm_dt=B.gbm_dt, retail_arrival_rate=rate,
                   retail_mean_size=size, retail_size_sigma=B.retail_size_sigma, retail_buy_prob=B.retail_buy_prob, seed=i))
    return out


def main():
    outp, s0, n = sys.argv[1], int(sys.argv[2]), int(sys.argv[3])
    seeds = list(range(s0, s0 + n)); cfgs = configs(seeds)
    van, _ = get_vanilla_bytecode_and_abi()
    res = {}
    for path in sys.argv[4:]:
        c = SolidityCompiler().compile(open(path).read())
        assert c.success, c.errors
        t0 = time.time()
        b = amm_sim_rs.run_batch(list(c.bytecode), list(van), cfgs, 8)
        edges = [r.edges.get('submission', 0.0) for r in b.results]
        gs = amm_sim_rs.gas_stats()
        res[os.path.basename(path)] = dict(edges=edges, seeds=seeds, deployed=len(c.deployed_bytecode), gas=str(gs), secs=time.time() - t0)
        print(os.path.basename(path), 'mean edge', sum(edges) / n, 'deployed', len(c.deployed_bytecode), 'gas', gs, f'{time.time()-t0:.0f}s', flush=True)
    json.dump(res, open(outp, 'w'))


main()

```

### qcheck.py (replica coarse-vs-full)

```python
"""Paired replica check: real.py full tables vs the same filter with the coarse table (ref.py).
usage: python3 -I qcheck.py START N"""
import sys, os, math
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import ref
from multiprocessing import Pool
R = ref.R


def one(seed):
    T = ref.T_full()
    a = R.run(seed, ref.real.RealPolicy(T, {'soft': 1.0}))['edge']
    b = R.run(seed, ref.RefPolicy(T))['edge']
    return seed, a, b


if __name__ == '__main__':
    s0, n = int(sys.argv[1]), int(sys.argv[2])
    with Pool(os.cpu_count()) as p:
        res = p.map(one, range(s0, s0 + n), chunksize=1)
    d = [b - a for _, a, b in res]
    m = sum(d) / n; se = (sum((x - m) ** 2 for x in d) / (n - 1)) ** .5 / n ** .5
    print(f'seeds {s0}..{s0+n-1}: full {sum(r[1] for r in res)/n:.3f} coarse {sum(r[2] for r in res)/n:.3f} coarse-full {m:+.4f} se {se:.4f} worst {min(d):+.3f}')

```