diff --git a/.github/workflows/.pre-commit-config.yaml b/.github/workflows/.pre-commit-config.yaml
new file mode 100644
index 00000000..3c39b1a2
--- /dev/null
+++ b/.github/workflows/.pre-commit-config.yaml
@@ -0,0 +1,11 @@
+name: pre-commit
+
+on: [push]
+
+jobs:
+  pre-commit:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - uses: actions/setup-python@v4
+      - uses: pre-commit/action@v3.0.0
diff --git a/.github/workflows/unit-tests.yaml b/.github/workflows/unit-tests.yaml
new file mode 100644
index 00000000..d02b9553
--- /dev/null
+++ b/.github/workflows/unit-tests.yaml
@@ -0,0 +1,30 @@
+name: unit-tests-boa
+
+on: [push]
+
+env:
+  GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+
+jobs:
+  boa-tests:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+
+      - name: Cache Compiler Installations
+        uses: actions/cache@v2
+        with:
+          path: |
+            ~/.vvm
+          key: compiler-cache
+
+      - name: Setup Python 3.10.4
+        uses: actions/setup-python@v2
+        with:
+          python-version: 3.10.4
+
+      - name: Install Requirements
+        run: pip install -r requirements.txt
+
+      - name: Run Tests
+        run: python -m pytest tests/ -n auto
diff --git a/ape-config.yaml b/ape-config.yaml
deleted file mode 100644
index 33b4e7ae..00000000
--- a/ape-config.yaml
+++ /dev/null
@@ -1,59 +0,0 @@
-name: curve-tricrypto-ng
-contracts_folder: contracts/main/
-
-plugins:
-  - name: vyper
-  - name: alchemy
-  - name: hardhat
-  - name: ledger
-  - name: etherscan
-  - name: arbitrum
-  - name: optimism
-  - name: polygon
-
-default_ecosystem: ethereum
-
-# vyper:
-#   evm_version: paris  # enable for non PUSH0 evm networks
-
-hardhat:
-  port: auto
-  fork:
-    ethereum:
-      mainnet:
-        upstream_provider: alchemy
-      sepolia:
-        upstream_provider: alchemy
-    arbitrum:
-      mainnet:
-        upstream_provider: geth
-
-ethereum:
-  default_network: mainnet-fork
-  mainnet_fork:
-    transaction_acceptance_timeout: 99999999
-    default_provider: hardhat
-  mainnet:
-    transaction_acceptance_timeout: 99999999
-  sepolia:
-    transaction_acceptance_timeout: 99999999
-
-arbitrum:
-  default_network: mainnet-fork
-  mainnet_fork:
-    transaction_acceptance_timeout: 99999999
-    default_provider: hardhat
-  mainnet:
-    transaction_acceptance_timeout: 99999999
-
-geth:
-  ethereum:
-    mainnet:
-      uri: http://localhost:9090
-  arbitrum:
-    mainnet:
-      uri: https://arb-mainnet.g.alchemy.com/v2/{some_key}
-
-test:
-  mnemonic: test test test test test test test test test test test junk
-  number_of_accounts: 5
diff --git a/contracts/experimental/n=2.vy b/contracts/experimental/n=2.vy
new file mode 100644
index 00000000..d7ddbc70
--- /dev/null
+++ b/contracts/experimental/n=2.vy
@@ -0,0 +1,195 @@
+# @version 0.3.10
+
+
+N_COINS: constant(uint256) = 2
+PRECISION: constant(uint256) = 10 ** 18  # The precision to convert to
+A_MULTIPLIER: constant(uint256) = 10000
+
+MIN_GAMMA: constant(uint256) = 10**10
+MAX_GAMMA: constant(uint256) = 2 * 10**15
+
+MIN_A: constant(uint256) = N_COINS**N_COINS * A_MULTIPLIER / 10
+MAX_A: constant(uint256) = N_COINS**N_COINS * A_MULTIPLIER * 1000
+
+
+@internal
+@pure
+def _newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i: uint256) -> uint256:
+    """
+    Calculating x[i] given other balances x[0..N_COINS-1] and invariant D
+    ANN = A * N**N
+    """
+    # Safety checks
+    assert ANN > MIN_A - 1 and ANN < MAX_A + 1  # dev: unsafe values A
+    assert gamma > MIN_GAMMA - 1 and gamma < MAX_GAMMA + 1  # dev: unsafe values gamma
+    assert D > 10**17 - 1 and D < 10**15 * 10**18 + 1 # dev: unsafe values D
+
+    x_j: uint256 = x[1 - i]
+    y: uint256 = D**2 / (x_j * N_COINS**2)
+    K0_i: uint256 = (10**18 * N_COINS) * x_j / D
+    # S_i = x_j
+
+    # frac = x_j * 1e18 / D => frac = K0_i / N_COINS
+    # assert (K0_i > 10**16*N_COINS - 1) and (K0_i < 10**20*N_COINS + 1)  # dev: unsafe values x[i]
+
+    # x_sorted: uint256[N_COINS] = x
+    # x_sorted[i] = 0
+    # x_sorted = self.sort(x_sorted)  # From high to low
+    # x[not i] instead of x_sorted since x_soted has only 1 element
+
+    convergence_limit: uint256 = max(max(x_j / 10**14, D / 10**14), 100)
+
+    for j in range(255):
+        y_prev: uint256 = y
+
+        K0: uint256 = K0_i * y * N_COINS / D
+        S: uint256 = x_j + y
+
+        _g1k0: uint256 = gamma + 10**18
+        if _g1k0 > K0:
+            _g1k0 = _g1k0 - K0 + 1
+        else:
+            _g1k0 = K0 - _g1k0 + 1
+
+        # D / (A * N**N) * _g1k0**2 / gamma**2
+        mul1: uint256 = 10**18 * D / gamma * _g1k0 / gamma * _g1k0 * A_MULTIPLIER / ANN
+
+        # 2*K0 / _g1k0
+        mul2: uint256 = 10**18 + (2 * 10**18) * K0 / _g1k0
+
+        yfprime: uint256 = 10**18 * y + S * mul2 + mul1
+        _dyfprime: uint256 = D * mul2
+        if yfprime < _dyfprime:
+            y = y_prev / 2
+            continue
+        else:
+            yfprime -= _dyfprime
+        fprime: uint256 = yfprime / y
+
+        # y -= f / f_prime;  y = (y * fprime - f) / fprime
+        # y = (yfprime + 10**18 * D - 10**18 * S) // fprime + mul1 // fprime * (10**18 - K0) // K0
+        y_minus: uint256 = mul1 / fprime
+        y_plus: uint256 = (yfprime + 10**18 * D) / fprime + y_minus * 10**18 / K0
+        y_minus += 10**18 * S / fprime
+
+        if y_plus < y_minus:
+            y = y_prev / 2
+        else:
+            y = y_plus - y_minus
+
+        diff: uint256 = 0
+        if y > y_prev:
+            diff = y - y_prev
+        else:
+            diff = y_prev - y
+        if diff < max(convergence_limit, y / 10**14):
+            frac: uint256 = y * 10**18 / D
+            assert (frac > 10**16 - 1) and (frac < 10**20 + 1)  # dev: unsafe value for y
+            return y
+
+    raise "Did not converge"
+
+@external
+@pure
+def newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i: uint256) -> uint256:
+    return self._newton_y(ANN, gamma, x, D, i)
+
+
+@internal
+@pure
+def geometric_mean(unsorted_x: uint256[N_COINS], sort: bool) -> uint256:
+    """
+    (x[0] * x[1] * ...) ** (1/N)
+    """
+    x: uint256[N_COINS] = unsorted_x
+    if sort and x[0] < x[1]:
+        x = [unsorted_x[1], unsorted_x[0]]
+    D: uint256 = x[0]
+    diff: uint256 = 0
+    for i in range(255):
+        D_prev: uint256 = D
+        # tmp: uint256 = 10**18
+        # for _x in x:
+        #     tmp = tmp * _x / D
+        # D = D * ((N_COINS - 1) * 10**18 + tmp) / (N_COINS * 10**18)
+        # line below makes it for 2 coins
+        D = unsafe_div(D + x[0] * x[1] / D, N_COINS)
+        if D > D_prev:
+            diff = unsafe_sub(D, D_prev)
+        else:
+            diff = unsafe_sub(D_prev, D)
+        if diff <= 1 or diff * 10**18 < D:
+            return D
+    raise "Did not converge"
+
+
+@external
+@view
+def newton_D(ANN: uint256, gamma: uint256, x_unsorted: uint256[N_COINS]) -> uint256:
+    """
+    Finding the invariant using Newton method.
+    ANN is higher by the factor A_MULTIPLIER
+    ANN is already A * N**N
+
+    Currently uses 60k gas
+    """
+
+    # Initial value of invariant D is that for constant-product invariant
+    x: uint256[N_COINS] = x_unsorted
+    if x[0] < x[1]:
+        x = [x_unsorted[1], x_unsorted[0]]
+
+    D: uint256 = N_COINS * self.geometric_mean(x, False)
+    S: uint256 = x[0] + x[1]
+    __g1k0: uint256 = gamma + 10**18
+
+    for i in range(255):
+        D_prev: uint256 = D
+        assert D > 0
+        # Unsafe ivision by D is now safe
+
+        # K0: uint256 = 10**18
+        # for _x in x:
+        #     K0 = K0 * _x * N_COINS / D
+        # collapsed for 2 coins
+        K0: uint256 = unsafe_div(unsafe_div((10**18 * N_COINS**2) * x[0], D) * x[1], D)
+
+        _g1k0: uint256 = __g1k0
+        if _g1k0 > K0:
+            _g1k0 = unsafe_sub(_g1k0, K0) + 1  # > 0
+        else:
+            _g1k0 = unsafe_sub(K0, _g1k0) + 1  # > 0
+
+        # D / (A * N**N) * _g1k0**2 / gamma**2
+        mul1: uint256 = unsafe_div(unsafe_div(unsafe_div(10**18 * D, gamma) * _g1k0, gamma) * _g1k0 * A_MULTIPLIER, ANN)
+
+        # 2*N*K0 / _g1k0
+        mul2: uint256 = unsafe_div(((2 * 10**18) * N_COINS) * K0, _g1k0)
+
+        neg_fprime: uint256 = (S + unsafe_div(S * mul2, 10**18)) + mul1 * N_COINS / K0 - unsafe_div(mul2 * D, 10**18)
+
+        # D -= f / fprime
+        D_plus: uint256 = D * (neg_fprime + S) / neg_fprime
+        D_minus: uint256 = D*D / neg_fprime
+        if 10**18 > K0:
+            D_minus += unsafe_div(D * (mul1 / neg_fprime), 10**18) * unsafe_sub(10**18, K0) / K0
+        else:
+            D_minus -= unsafe_div(D * (mul1 / neg_fprime), 10**18) * unsafe_sub(K0, 10**18) / K0
+
+        if D_plus > D_minus:
+            D = unsafe_sub(D_plus, D_minus)
+        else:
+            D = unsafe_div(unsafe_sub(D_minus, D_plus), 2)
+
+        diff: uint256 = 0
+        if D > D_prev:
+            diff = unsafe_sub(D, D_prev)
+        else:
+            diff = unsafe_sub(D_prev, D)
+        if diff * 10**14 < max(10**16, D):  # Could reduce precision for gas efficiency here
+            # Test that we are safe with the next newton_y
+            for _x in x:
+                frac: uint256 = _x * 10**18 / D
+            return D
+
+    raise "Did not converge"
diff --git a/contracts/experimental/n=2_optimized.vy b/contracts/experimental/n=2_optimized.vy
new file mode 100644
index 00000000..69a8ef13
--- /dev/null
+++ b/contracts/experimental/n=2_optimized.vy
@@ -0,0 +1,368 @@
+# @version 0.3.10
+
+
+N_COINS: constant(uint256) = 2
+PRECISION: constant(uint256) = 10 ** 18  # The precision to convert to
+A_MULTIPLIER: constant(uint256) = 10000
+
+
+@internal
+@pure
+def _cbrt(x: uint256) -> uint256:
+
+    # we artificially set a cap to the values for which we can compute the
+    # cube roots safely. This is not to say that there are no values above
+    # max(uint256) // 10**36 for which we cannot get good cube root estimates.
+    # However, beyond this point, accuracy is not guaranteed since overflows
+    # start to occur.
+    # assert x < 115792089237316195423570985008687907853269, "inaccurate cbrt"  # TODO: check limits again
+
+    # we increase precision of input `x` by multiplying 10 ** 36.
+    # in such cases: cbrt(10**18) = 10**18, cbrt(1) = 10**12
+    xx: uint256 = 0
+    if x >= 115792089237316195423570985008687907853269 * 10**18:
+        xx = x
+    elif x >= 115792089237316195423570985008687907853269:
+        xx = unsafe_mul(x, 10**18)
+    else:
+        xx = unsafe_mul(x, 10**36)
+
+    # get log2(x) for approximating initial value
+    # logic is: cbrt(a) = cbrt(2**(log2(a))) = 2**(log2(a) / 3) ≈ 2**|log2(a)/3|
+    # from: https://github.com/transmissions11/solmate/blob/b9d69da49bbbfd090f1a73a4dba28aa2d5ee199f/src/utils/FixedPointMathLib.sol#L352
+
+    a_pow: int256 = 0
+    if xx > 340282366920938463463374607431768211455:
+        a_pow = 128
+    if unsafe_div(xx, shift(2, a_pow)) > 18446744073709551615:
+        a_pow = a_pow | 64
+    if unsafe_div(xx, shift(2, a_pow)) > 4294967295:
+        a_pow = a_pow | 32
+    if unsafe_div(xx, shift(2, a_pow)) > 65535:
+        a_pow = a_pow | 16
+    if unsafe_div(xx, shift(2, a_pow)) > 255:
+        a_pow = a_pow | 8
+    if unsafe_div(xx, shift(2, a_pow)) > 15:
+        a_pow = a_pow | 4
+    if unsafe_div(xx, shift(2, a_pow)) > 3:
+        a_pow = a_pow | 2
+    if unsafe_div(xx, shift(2, a_pow)) > 1:
+        a_pow = a_pow | 1
+
+    # initial value: 2**|log2(a)/3|
+    # which is: 2 ** (n / 3) * 1260 ** (n % 3) / 1000 ** (n % 3)
+    a_pow_mod: uint256 = convert(a_pow, uint256) % 3
+    a: uint256 = unsafe_div(
+        unsafe_mul(
+            pow_mod256(
+                2,
+                unsafe_div(
+                    convert(a_pow, uint256), 3
+                )
+            ),
+            pow_mod256(1260, a_pow_mod)
+        ),
+        pow_mod256(1000, a_pow_mod)
+    )
+
+    # 7 newton raphson iterations:
+    a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
+    a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
+    a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
+    a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
+    a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
+    a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
+    a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
+
+    if x >= 115792089237316195423570985008687907853269 * 10**18:
+        return a*10**12
+    elif x >= 115792089237316195423570985008687907853269:
+        return a*10**6
+    else:
+        return a
+
+
+@internal
+@pure
+def _newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i: uint256) -> uint256:
+    """
+    Calculating x[i] given other balances x[0..N_COINS-1] and invariant D
+    ANN = A * N**N
+    """
+    # Safety checks
+    # assert ANN > MIN_A - 1 and ANN < MAX_A + 1  # dev: unsafe values A
+    # assert gamma > MIN_GAMMA - 1 and gamma < MAX_GAMMA + 1  # dev: unsafe values gamma
+    # assert D > 10**17 - 1 and D < 10**15 * 10**18 + 1 # dev: unsafe values D
+
+    x_j: uint256 = x[1 - i]
+    y: uint256 = D**2 / (x_j * N_COINS**2)
+    K0_i: uint256 = (10**18 * N_COINS) * x_j / D
+    # S_i = x_j
+
+    # frac = x_j * 1e18 / D => frac = K0_i / N_COINS
+    # assert (K0_i > 10**16*N_COINS - 1) and (K0_i < 10**20*N_COINS + 1)  # dev: unsafe values x[i]
+
+    # x_sorted: uint256[N_COINS] = x
+    # x_sorted[i] = 0
+    # x_sorted = self.sort(x_sorted)  # From high to low
+    # x[not i] instead of x_sorted since x_soted has only 1 element
+
+    convergence_limit: uint256 = max(max(x_j / 10**14, D / 10**14), 100)
+
+    for j in range(255):
+        y_prev: uint256 = y
+
+        K0: uint256 = K0_i * y * N_COINS / D
+        S: uint256 = x_j + y
+
+        _g1k0: uint256 = gamma + 10**18
+        if _g1k0 > K0:
+            _g1k0 = _g1k0 - K0 + 1
+        else:
+            _g1k0 = K0 - _g1k0 + 1
+
+        # D / (A * N**N) * _g1k0**2 / gamma**2
+        mul1: uint256 = 10**18 * D / gamma * _g1k0 / gamma * _g1k0 * A_MULTIPLIER / ANN
+
+        # 2*K0 / _g1k0
+        mul2: uint256 = 10**18 + (2 * 10**18) * K0 / _g1k0
+
+        yfprime: uint256 = 10**18 * y + S * mul2 + mul1
+        _dyfprime: uint256 = D * mul2
+        if yfprime < _dyfprime:
+            y = y_prev / 2
+            continue
+        else:
+            yfprime -= _dyfprime
+        fprime: uint256 = yfprime / y
+
+        # y -= f / f_prime;  y = (y * fprime - f) / fprime
+        # y = (yfprime + 10**18 * D - 10**18 * S) // fprime + mul1 // fprime * (10**18 - K0) // K0
+        y_minus: uint256 = mul1 / fprime
+        y_plus: uint256 = (yfprime + 10**18 * D) / fprime + y_minus * 10**18 / K0
+        y_minus += 10**18 * S / fprime
+
+        if y_plus < y_minus:
+            y = y_prev / 2
+        else:
+            y = y_plus - y_minus
+
+        diff: uint256 = 0
+        if y > y_prev:
+            diff = y - y_prev
+        else:
+            diff = y_prev - y
+        if diff < max(convergence_limit, y / 10**14):
+            frac: uint256 = y * 10**18 / D
+            # assert (frac > 10**16 - 1) and (frac < 10**20 + 1)  # dev: unsafe value for y
+            return y
+
+    raise "Did not converge"
+
+@external
+@pure
+def newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i: uint256) -> uint256:
+    return self._newton_y(ANN, gamma, x, D, i)
+
+
+@external
+@pure
+def get_y(_ANN: uint256, _gamma: uint256, _x: uint256[N_COINS], _D: uint256, i: uint256) -> uint256[2]:
+
+    # Safety checks go here
+
+    j: uint256 = 0
+    if i == 0:
+        j = 1
+    elif i == 1:
+        j = 0
+
+    ANN: int256 = convert(_ANN, int256)
+    gamma: int256 = convert(_gamma, int256)
+    D: int256 = convert(_D, int256)
+    x_j: int256 = convert(_x[j], int256)
+    gamma2: int256 = gamma**2
+
+    a: int256 = 10**32
+    b: int256 = ANN*D*gamma2/4/10000/x_j/10**4 - 10**32*3 - 2*gamma*10**14
+    c: int256 = 10**32*3 + 4*gamma*10**14 + gamma2/10**4 + 4*ANN*gamma2*x_j/D/10000/4/10**4 - 4*ANN*gamma2/10000/4/10**4
+    d: int256 = -(10**18+gamma)**2 / 10**4
+
+    delta0: int256 = 3*a*c/b - b
+    delta1: int256 = 9*a*c/b - 2*b - 27*a**2/b*d/b
+
+    divider: int256 = 0
+    threshold: int256 = min(min(abs(delta0), abs(delta1)), a)
+    if threshold > 10**48:
+        divider = 10**30
+    elif threshold > 10**46:
+        divider = 10**28
+    elif threshold > 10**44:
+        divider = 10**26
+    elif threshold > 10**42:
+        divider = 10**24
+    elif threshold > 10**40:
+        divider = 10**22
+    elif threshold > 10**38:
+        divider = 10**20
+    elif threshold > 10**36:
+        divider = 10**18
+    elif threshold > 10**34:
+        divider = 10**16
+    elif threshold > 10**32:
+        divider = 10**14
+    elif threshold > 10**30:
+        divider = 10**12
+    elif threshold > 10**28:
+        divider = 10**10
+    elif threshold > 10**26:
+        divider = 10**8
+    elif threshold > 10**24:
+        divider = 10**6
+    elif threshold > 10**20:
+        divider = 10**2
+    else:
+        divider = 1
+
+    a /= divider
+    b /= divider
+    c /= divider
+    d /= divider
+
+    delta0 = 3*a*c/b - b
+    delta1 = 9*a*c/b - 2*b - 27*a**2/b*d/b
+
+    sqrt_arg: int256 = delta1**2 + 4*delta0**2/b*delta0
+    sqrt_val: int256 = 0
+    if sqrt_arg > 0:
+        sqrt_val = convert(isqrt(convert(sqrt_arg, uint256)), int256)
+    else:
+        return [self._newton_y(_ANN, _gamma, _x, _D, i), 0]
+
+    b_cbrt: int256 = 0
+    if b > 0:
+        b_cbrt = convert(self._cbrt(convert(b, uint256)), int256)
+    else:
+        b_cbrt = -convert(self._cbrt(convert(-b, uint256)), int256)
+
+    second_cbrt: int256 = 0
+    if delta1 > 0:
+        second_cbrt = convert(self._cbrt(convert((delta1 + sqrt_val), uint256) / 2), int256)
+    else:
+        second_cbrt = -convert(self._cbrt(convert(-(delta1 - sqrt_val), uint256) / 2), int256)
+
+    C1: int256 = b_cbrt*b_cbrt/10**18*second_cbrt/10**18
+
+    root: int256 = -(10**18*b - 10**18*C1 + 10**18*b*delta0/C1)/(3*a)
+
+    return [convert(D**2/x_j*root/4/10**18, uint256), convert(root, uint256)]
+
+
+@internal
+@pure
+def geometric_mean(unsorted_x: uint256[N_COINS], sort: bool) -> uint256:
+    """
+    (x[0] * x[1] * ...) ** (1/N)
+    """
+    x: uint256[N_COINS] = unsorted_x
+    if sort and x[0] < x[1]:
+        x = [unsorted_x[1], unsorted_x[0]]
+    D: uint256 = x[0]
+    diff: uint256 = 0
+    for i in range(255):
+        D_prev: uint256 = D
+        # tmp: uint256 = 10**18
+        # for _x in x:
+        #     tmp = tmp * _x / D
+        # D = D * ((N_COINS - 1) * 10**18 + tmp) / (N_COINS * 10**18)
+        # line below makes it for 2 coins
+        D = unsafe_div(D + x[0] * x[1] / D, N_COINS)
+        if D > D_prev:
+            diff = unsafe_sub(D, D_prev)
+        else:
+            diff = unsafe_sub(D_prev, D)
+        if diff <= 1 or diff * 10**18 < D:
+            return D
+    raise "Did not converge"
+
+
+@external
+@view
+def newton_D(ANN: uint256, gamma: uint256, x_unsorted: uint256[N_COINS], K0_prev: uint256 = 0) -> uint256:
+    """
+    Finding the invariant using Newton method.
+    ANN is higher by the factor A_MULTIPLIER
+    ANN is already A * N**N
+
+    Currently uses 60k gas
+    """
+
+    # Initial value of invariant D is that for constant-product invariant
+    x: uint256[N_COINS] = x_unsorted
+    if x[0] < x[1]:
+        x = [x_unsorted[1], x_unsorted[0]]
+
+    S: uint256 = x[0] + x[1]
+
+    D: uint256 = 0
+    if K0_prev == 0:
+        # Geometric mean of 3 numbers cannot be larger than the largest number
+        # so the following is safe to do:
+        D = N_COINS * self.geometric_mean(x, False)
+    else:
+        D = isqrt(x[0] * x[1] * 4 / K0_prev * 10**18)
+        if S < D:
+            D = S
+
+    __g1k0: uint256 = gamma + 10**18
+
+    for i in range(255):
+        D_prev: uint256 = D
+        assert D > 0
+        # Unsafe ivision by D is now safe
+
+        # K0: uint256 = 10**18
+        # for _x in x:
+        #     K0 = K0 * _x * N_COINS / D
+        # collapsed for 2 coins
+        K0: uint256 = unsafe_div(unsafe_div((10**18 * N_COINS**2) * x[0], D) * x[1], D)
+
+        _g1k0: uint256 = __g1k0
+        if _g1k0 > K0:
+            _g1k0 = unsafe_sub(_g1k0, K0) + 1  # > 0
+        else:
+            _g1k0 = unsafe_sub(K0, _g1k0) + 1  # > 0
+
+        # D / (A * N**N) * _g1k0**2 / gamma**2
+        mul1: uint256 = unsafe_div(unsafe_div(unsafe_div(10**18 * D, gamma) * _g1k0, gamma) * _g1k0 * A_MULTIPLIER, ANN)
+
+        # 2*N*K0 / _g1k0
+        mul2: uint256 = unsafe_div(((2 * 10**18) * N_COINS) * K0, _g1k0)
+
+        neg_fprime: uint256 = (S + unsafe_div(S * mul2, 10**18)) + mul1 * N_COINS / K0 - unsafe_div(mul2 * D, 10**18)
+
+        # D -= f / fprime
+        D_plus: uint256 = D * (neg_fprime + S) / neg_fprime
+        D_minus: uint256 = D*D / neg_fprime
+        if 10**18 > K0:
+            D_minus += unsafe_div(D * (mul1 / neg_fprime), 10**18) * unsafe_sub(10**18, K0) / K0
+        else:
+            D_minus -= unsafe_div(D * (mul1 / neg_fprime), 10**18) * unsafe_sub(K0, 10**18) / K0
+
+        if D_plus > D_minus:
+            D = unsafe_sub(D_plus, D_minus)
+        else:
+            D = unsafe_div(unsafe_sub(D_minus, D_plus), 2)
+
+        diff: uint256 = 0
+        if D > D_prev:
+            diff = unsafe_sub(D, D_prev)
+        else:
+            diff = unsafe_sub(D_prev, D)
+        if diff * 10**14 < max(10**16, D):  # Could reduce precision for gas efficiency here
+            # Test that we are safe with the next newton_y
+            for _x in x:
+                frac: uint256 = _x * 10**18 / D
+            return D
+
+    raise "Did not converge"
diff --git a/contracts/main/CurveCryptoMathOptimized2.vy b/contracts/main/CurveCryptoMathOptimized2.vy
index fd399394..9b155388 100644
--- a/contracts/main/CurveCryptoMathOptimized2.vy
+++ b/contracts/main/CurveCryptoMathOptimized2.vy
@@ -1,5 +1,5 @@
-# @version 0.3.10
-#pragma optimize gas
+# pragma version 0.3.10
+# pragma optimize gas
 # (c) Curve.Fi, 2020-2023
 # AMM Math for 2-coin Curve Cryptoswap Pools
 #
@@ -7,7 +7,7 @@
 # Swiss Stake GmbH are allowed to call this contract.
 
 """
-@title CurveTricryptoMathOptimized
+@title CurveTwocryptoMathOptimized
 @author Curve.Fi
 @license Copyright (c) Curve.Fi, 2020-2023 - all rights reserved
 @notice Curve AMM Math for 2 unpegged assets (e.g. ETH <> USD).
@@ -17,10 +17,10 @@ N_COINS: constant(uint256) = 2
 A_MULTIPLIER: constant(uint256) = 10000
 
 MIN_GAMMA: constant(uint256) = 10**10
-MAX_GAMMA: constant(uint256) = 5 * 10**16
+MAX_GAMMA: constant(uint256) = 2 * 10**15
 
 MIN_A: constant(uint256) = N_COINS**N_COINS * A_MULTIPLIER / 10
-MAX_A: constant(uint256) = N_COINS**N_COINS * A_MULTIPLIER * 100000
+MAX_A: constant(uint256) = N_COINS**N_COINS * A_MULTIPLIER * 1000
 
 version: public(constant(String[8])) = "v2.0.0"
 
@@ -30,17 +30,60 @@ version: public(constant(String[8])) = "v2.0.0"
 
 @internal
 @pure
-def _cbrt(x: uint256) -> uint256:
+def _snekmate_log_2(x: uint256, roundup: bool) -> uint256:
+    """
+    @notice An `internal` helper function that returns the log in base 2
+         of `x`, following the selected rounding direction.
+    @dev This implementation is derived from Snekmate, which is authored
+         by pcaversaccio (Snekmate), distributed under the AGPL-3.0 license.
+         https://github.com/pcaversaccio/snekmate
+    @dev Note that it returns 0 if given 0. The implementation is
+         inspired by OpenZeppelin's implementation here:
+         https://github.com/OpenZeppelin/openzeppelin-contracts/blob/master/contracts/utils/math/Math.sol.
+    @param x The 32-byte variable.
+    @param roundup The Boolean variable that specifies whether
+           to round up or not. The default `False` is round down.
+    @return uint256 The 32-byte calculation result.
+    """
+    value: uint256 = x
+    result: uint256 = empty(uint256)
+
+    # The following lines cannot overflow because we have the well-known
+    # decay behaviour of `log_2(max_value(uint256)) < max_value(uint256)`.
+    if x >> 128 != empty(uint256):
+        value = x >> 128
+        result = 128
+    if value >> 64 != empty(uint256):
+        value = value >> 64
+        result = unsafe_add(result, 64)
+    if value >> 32 != empty(uint256):
+        value = value >> 32
+        result = unsafe_add(result, 32)
+    if value >> 16 != empty(uint256):
+        value = value >> 16
+        result = unsafe_add(result, 16)
+    if value >> 8 != empty(uint256):
+        value = value >> 8
+        result = unsafe_add(result, 8)
+    if value >> 4 != empty(uint256):
+        value = value >> 4
+        result = unsafe_add(result, 4)
+    if value >> 2 != empty(uint256):
+        value = value >> 2
+        result = unsafe_add(result, 2)
+    if value >> 1 != empty(uint256):
+        result = unsafe_add(result, 1)
+
+    if (roundup and (1 << result) < x):
+        result = unsafe_add(result, 1)
+
+    return result
+
 
-    # we artificially set a cap to the values for which we can compute the
-    # cube roots safely. This is not to say that there are no values above
-    # max(uint256) // 10**36 for which we cannot get good cube root estimates.
-    # However, beyond this point, accuracy is not guaranteed since overflows
-    # start to occur.
-    # assert x < 115792089237316195423570985008687907853269, "inaccurate cbrt"  # TODO: check limits again
+@internal
+@pure
+def _cbrt(x: uint256) -> uint256:
 
-    # we increase precision of input `x` by multiplying 10 ** 36.
-    # in such cases: cbrt(10**18) = 10**18, cbrt(1) = 10**12
     xx: uint256 = 0
     if x >= 115792089237316195423570985008687907853269 * 10**18:
         xx = x
@@ -49,45 +92,34 @@ def _cbrt(x: uint256) -> uint256:
     else:
         xx = unsafe_mul(x, 10**36)
 
-    # get log2(x) for approximating initial value
-    # logic is: cbrt(a) = cbrt(2**(log2(a))) = 2**(log2(a) / 3) ≈ 2**|log2(a)/3|
-    # from: https://github.com/transmissions11/solmate/blob/b9d69da49bbbfd090f1a73a4dba28aa2d5ee199f/src/utils/FixedPointMathLib.sol#L352
-
-    a_pow: int256 = 0
-    if xx > 340282366920938463463374607431768211455:
-        a_pow = 128
-    if unsafe_div(xx, shift(2, a_pow)) > 18446744073709551615:
-        a_pow = a_pow | 64
-    if unsafe_div(xx, shift(2, a_pow)) > 4294967295:
-        a_pow = a_pow | 32
-    if unsafe_div(xx, shift(2, a_pow)) > 65535:
-        a_pow = a_pow | 16
-    if unsafe_div(xx, shift(2, a_pow)) > 255:
-        a_pow = a_pow | 8
-    if unsafe_div(xx, shift(2, a_pow)) > 15:
-        a_pow = a_pow | 4
-    if unsafe_div(xx, shift(2, a_pow)) > 3:
-        a_pow = a_pow | 2
-    if unsafe_div(xx, shift(2, a_pow)) > 1:
-        a_pow = a_pow | 1
-
-    # initial value: 2**|log2(a)/3|
-    # which is: 2 ** (n / 3) * 1260 ** (n % 3) / 1000 ** (n % 3)
-    a_pow_mod: uint256 = convert(a_pow, uint256) % 3
+    log2x: int256 = convert(self._snekmate_log_2(xx, False), int256)
+
+    # When we divide log2x by 3, the remainder is (log2x % 3).
+    # So if we just multiply 2**(log2x/3) and discard the remainder to calculate our
+    # guess, the newton method will need more iterations to converge to a solution,
+    # since it is missing that precision. It's a few more calculations now to do less
+    # calculations later:
+    # pow = log2(x) // 3
+    # remainder = log2(x) % 3
+    # initial_guess = 2 ** pow * cbrt(2) ** remainder
+    # substituting -> 2 = 1.26 ≈ 1260 / 1000, we get:
+    #
+    # initial_guess = 2 ** pow * 1260 ** remainder // 1000 ** remainder
+
+    remainder: uint256 = convert(log2x, uint256) % 3
     a: uint256 = unsafe_div(
         unsafe_mul(
-            pow_mod256(
-                2,
-                unsafe_div(
-                    convert(a_pow, uint256), 3
-                )
-            ),
-            pow_mod256(1260, a_pow_mod)
+            pow_mod256(2, unsafe_div(convert(log2x, uint256), 3)),  # <- pow
+            pow_mod256(1260, remainder),
         ),
-        pow_mod256(1000, a_pow_mod)
+        pow_mod256(1000, remainder),
     )
 
-    # 7 newton raphson iterations:
+    # Because we chose good initial values for cube roots, 7 newton raphson iterations
+    # are just about sufficient. 6 iterations would result in non-convergences, and 8
+    # would be one too many iterations. Without initial values, the iteration count
+    # can go up to 20 or greater. The iterations are unrolled. This reduces gas costs
+    # but takes up more bytecode:
     a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
     a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
     a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
@@ -97,11 +129,11 @@ def _cbrt(x: uint256) -> uint256:
     a = unsafe_div(unsafe_add(unsafe_mul(2, a), unsafe_div(xx, unsafe_mul(a, a))), 3)
 
     if x >= 115792089237316195423570985008687907853269 * 10**18:
-        return a*10**12
+        a = unsafe_mul(a, 10**12)
     elif x >= 115792089237316195423570985008687907853269:
-        return a*10**6
-    else:
-        return a
+        a = unsafe_mul(a, 10**6)
+
+    return a
 
 
 @internal
@@ -110,24 +142,14 @@ def _newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i:
     """
     Calculating x[i] given other balances x[0..N_COINS-1] and invariant D
     ANN = A * N**N
+    This is computationally expensive.
     """
-    # Safety checks
-    # assert ANN > MIN_A - 1 and ANN < MAX_A + 1  # dev: unsafe values A
-    # assert gamma > MIN_GAMMA - 1 and gamma < MAX_GAMMA + 1  # dev: unsafe values gamma
-    # assert D > 10**17 - 1 and D < 10**15 * 10**18 + 1 # dev: unsafe values D
 
     x_j: uint256 = x[1 - i]
     y: uint256 = D**2 / (x_j * N_COINS**2)
     K0_i: uint256 = (10**18 * N_COINS) * x_j / D
-    # S_i = x_j
 
-    # frac = x_j * 1e18 / D => frac = K0_i / N_COINS
-    # assert (K0_i > 10**16*N_COINS - 1) and (K0_i < 10**20*N_COINS + 1)  # dev: unsafe values x[i]
-
-    # x_sorted: uint256[N_COINS] = x
-    # x_sorted[i] = 0
-    # x_sorted = self.sort(x_sorted)  # From high to low
-    # x[not i] instead of x_sorted since x_soted has only 1 element
+    assert (K0_i > 10**16*N_COINS - 1) and (K0_i < 10**20*N_COINS + 1)  # dev: unsafe values x[i]
 
     convergence_limit: uint256 = max(max(x_j / 10**14, D / 10**14), 100)
 
@@ -174,48 +196,88 @@ def _newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i:
             diff = y - y_prev
         else:
             diff = y_prev - y
+
         if diff < max(convergence_limit, y / 10**14):
-            frac: uint256 = y * 10**18 / D
-            # assert (frac > 10**16 - 1) and (frac < 10**20 + 1)  # dev: unsafe value for y
             return y
 
     raise "Did not converge"
 
+
 @external
 @pure
 def newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i: uint256) -> uint256:
-    return self._newton_y(ANN, gamma, x, D, i)
+
+    # Safety checks
+    assert ANN > MIN_A - 1 and ANN < MAX_A + 1  # dev: unsafe values A
+    assert gamma > MIN_GAMMA - 1 and gamma < MAX_GAMMA + 1  # dev: unsafe values gamma
+    assert D > 10**17 - 1 and D < 10**15 * 10**18 + 1 # dev: unsafe values D
+
+    y: uint256 = self._newton_y(ANN, gamma, x, D, i)
+    frac: uint256 = y * 10**18 / D
+    assert (frac >= 10**16 - 1) and (frac < 10**20 + 1)  # dev: unsafe value for y
+
+    return y
 
 
 @external
 @pure
-def get_y(_ANN: uint256, _gamma: uint256, _x: uint256[N_COINS], _D: uint256, i: uint256) -> uint256[2]:
+def get_y(
+    _ANN: uint256,
+    _gamma: uint256,
+    _x: uint256[N_COINS],
+    _D: uint256,
+    i: uint256
+) -> uint256[2]:
 
-    # Safety checks go here
-
-    j: uint256 = 0
-    if i == 0:
-        j = 1
-    elif i == 1:
-        j = 0
+    # Safety checks
+    assert _ANN > MIN_A - 1 and _ANN < MAX_A + 1  # dev: unsafe values A
+    assert _gamma > MIN_GAMMA - 1 and _gamma < MAX_GAMMA + 1  # dev: unsafe values gamma
+    assert _D > 10**17 - 1 and _D < 10**15 * 10**18 + 1 # dev: unsafe values D
 
     ANN: int256 = convert(_ANN, int256)
     gamma: int256 = convert(_gamma, int256)
     D: int256 = convert(_D, int256)
-    x_j: int256 = convert(_x[j], int256)
-    gamma2: int256 = gamma**2
+    x_j: int256 = convert(_x[1 - i], int256)
+    gamma2: int256 = unsafe_mul(gamma, gamma)
 
+    # savediv by x_j done here:
+    y: int256 = D**2 / (x_j * N_COINS**2)
+
+    # K0_i: int256 = (10**18 * N_COINS) * x_j / D
+    K0_i: int256 = unsafe_div(10**18 * N_COINS * x_j, D)
+    assert (K0_i > 10**16 * N_COINS - 1) and (K0_i < 10**20 * N_COINS + 1)  # dev: unsafe values x[i]
+
+    ann_gamma2: int256 = ANN * gamma2
+
+    # a = 10**36 / N_COINS**2
     a: int256 = 10**32
-    b: int256 = ANN*D*gamma2/4/10000/x_j/10**4 - 10**32*3 - 2*gamma*10**14
-    c: int256 = 10**32*3 + 4*gamma*10**14 + gamma2/10**4 + 4*ANN*gamma2*x_j/D/10000/4/10**4 - 4*ANN*gamma2/10000/4/10**4
-    d: int256 = -(10**18+gamma)**2 / 10**4
+
+    # b = ANN*D*gamma2/4/10000/x_j/10**4 - 10**32*3 - 2*gamma*10**14
+    b: int256 = (
+        D*ann_gamma2/400000000/x_j
+        - convert(unsafe_mul(10**32, 3), int256)
+        - unsafe_mul(unsafe_mul(2, gamma), 10**14)
+    )
+
+    # c = 10**32*3 + 4*gamma*10**14 + gamma2/10**4 + 4*ANN*gamma2*x_j/D/10000/4/10**4 - 4*ANN*gamma2/10000/4/10**4
+    c: int256 = (
+        unsafe_mul(10**32, convert(3, int256))
+        + unsafe_mul(unsafe_mul(4, gamma), 10**14)
+        + unsafe_div(gamma2, 10**4)
+        + unsafe_div(unsafe_div(unsafe_mul(4, ann_gamma2), 400000000) * x_j, D)
+        - unsafe_div(unsafe_mul(4, ann_gamma2), 400000000)
+    )
+
+    # d = -(10**18+gamma)**2 / 10**4
+    d: int256 = -unsafe_div(unsafe_add(10**18, gamma) ** 2, 10**4)
 
     # delta0: int256 = 3*a*c/b - b
-    delta0: int256 = unsafe_sub(unsafe_div(unsafe_mul(unsafe_mul(3, a), c), b), b)
+    delta0: int256 = 3 * a * c / b - b  # safediv by b
+
     # delta1: int256 = 9*a*c/b - 2*b - 27*a**2/b*d/b
-    delta1: int256 = unsafe_div(unsafe_mul(unsafe_mul(9, a), c), b) - unsafe_mul(2, b) - unsafe_div(unsafe_mul(unsafe_div(unsafe_mul(27, a**2), b), d), b)
+    delta1: int256 = 3 * delta0 + b - 27*a**2/b*d/b
 
-    divider: int256 = 0
+    divider: int256 = 1
     threshold: int256 = min(min(abs(delta0), abs(delta1)), a)
     if threshold > 10**48:
         divider = 10**30
@@ -245,18 +307,17 @@ def get_y(_ANN: uint256, _gamma: uint256, _x: uint256[N_COINS], _D: uint256, i:
         divider = 10**6
     elif threshold > 10**20:
         divider = 10**2
-    else:
-        divider = 1
 
     a = unsafe_div(a, divider)
     b = unsafe_div(b, divider)
     c = unsafe_div(c, divider)
     d = unsafe_div(d, divider)
 
-    # delta0 = 3*a*c/b - b
-    delta0 = unsafe_sub(unsafe_div(unsafe_mul(unsafe_mul(3, a), c), b), b)
+    # delta0 = 3*a*c/b - b: here we can do more unsafe ops now:
+    delta0 = unsafe_div(unsafe_mul(unsafe_mul(3, a), c), b) - b
+
     # delta1 = 9*a*c/b - 2*b - 27*a**2/b*d/b
-    delta1 = unsafe_div(unsafe_mul(unsafe_mul(9, a), c), b) - unsafe_mul(2, b) - unsafe_div(unsafe_mul(unsafe_div(unsafe_mul(27, a**2), b), d), b)
+    delta1 = 3 * delta0 + b - unsafe_div(unsafe_mul(unsafe_div(unsafe_mul(27, a**2), b), d), b)
 
     # sqrt_arg: int256 = delta1**2 + 4*delta0**2/b*delta0
     sqrt_arg: int256 = delta1**2 + unsafe_mul(unsafe_div(4*delta0**2, b), delta0)
@@ -264,7 +325,10 @@ def get_y(_ANN: uint256, _gamma: uint256, _x: uint256[N_COINS], _D: uint256, i:
     if sqrt_arg > 0:
         sqrt_val = convert(isqrt(convert(sqrt_arg, uint256)), int256)
     else:
-        return [self._newton_y(_ANN, _gamma, _x, _D, i), 0]
+        return [
+            self._newton_y(_ANN, _gamma, _x, _D, i),
+            0
+        ]
 
     b_cbrt: int256 = 0
     if b > 0:
@@ -284,38 +348,18 @@ def get_y(_ANN: uint256, _gamma: uint256, _x: uint256[N_COINS], _D: uint256, i:
     C1: int256 = unsafe_div(unsafe_mul(unsafe_div(b_cbrt**2, 10**18), second_cbrt), 10**18)
 
     # root: int256 = (10**18*C1 - 10**18*b - 10**18*b*delta0/C1)/(3*a), keep 2 safe ops here.
-    root: int256 = (unsafe_mul(10**18, C1) - unsafe_mul(10**18, b) - unsafe_div(unsafe_mul(10**18, b)*delta0, C1))/unsafe_mul(3, a)
+    root: int256 = (unsafe_mul(10**18, C1) - unsafe_mul(10**18, b) - unsafe_mul(10**18, b)/C1*delta0)/unsafe_mul(3, a)
 
-    # return [convert(D**2/x_j*root/4/10**18, uint256), convert(root, uint256)]
-    return [convert(unsafe_div(unsafe_div(unsafe_mul(unsafe_div(D**2, x_j), root), 4), 10**18), uint256), convert(root, uint256)]
+    # y_out: uint256[2] =  [
+    #     convert(D**2/x_j*root/4/10**18, uint256),   # <--- y
+    #     convert(root, uint256)  # <----------------------- K0Prev
+    # ]
+    y_out: uint256[2] = [convert(unsafe_div(unsafe_div(unsafe_mul(unsafe_div(D**2, x_j), root), 4), 10**18), uint256), convert(root, uint256)]
 
+    frac: uint256 = unsafe_div(y_out[0] * 10**18, _D)
+    assert (frac >= 10**16 - 1) and (frac < 10**20 + 1)  # dev: unsafe value for y
 
-@internal
-@pure
-def geometric_mean(unsorted_x: uint256[N_COINS], sort: bool) -> uint256:
-    """
-    (x[0] * x[1] * ...) ** (1/N)
-    """
-    x: uint256[N_COINS] = unsorted_x
-    if sort and x[0] < x[1]:
-        x = [unsorted_x[1], unsorted_x[0]]
-    D: uint256 = x[0]
-    diff: uint256 = 0
-    for i in range(255):
-        D_prev: uint256 = D
-        # tmp: uint256 = 10**18
-        # for _x in x:
-        #     tmp = tmp * _x / D
-        # D = D * ((N_COINS - 1) * 10**18 + tmp) / (N_COINS * 10**18)
-        # line below makes it for 2 coins
-        D = unsafe_div(D + x[0] * x[1] / D, N_COINS)
-        if D > D_prev:
-            diff = unsafe_sub(D, D_prev)
-        else:
-            diff = unsafe_sub(D_prev, D)
-        if diff <= 1 or diff * 10**18 < D:
-            return D
-    raise "Did not converge"
+    return y_out
 
 
 @external
@@ -325,16 +369,21 @@ def newton_D(ANN: uint256, gamma: uint256, x_unsorted: uint256[N_COINS], K0_prev
     Finding the invariant using Newton method.
     ANN is higher by the factor A_MULTIPLIER
     ANN is already A * N**N
-
-    Currently uses 60k gas
     """
 
+    # Safety checks
+    assert ANN > MIN_A - 1 and ANN < MAX_A + 1  # dev: unsafe values A
+    assert gamma > MIN_GAMMA - 1 and gamma < MAX_GAMMA + 1  # dev: unsafe values gamma
+
     # Initial value of invariant D is that for constant-product invariant
     x: uint256[N_COINS] = x_unsorted
     if x[0] < x[1]:
         x = [x_unsorted[1], x_unsorted[0]]
 
-    S: uint256 = x[0] + x[1]
+    assert x[0] > 10**9 - 1 and x[0] < 10**15 * 10**18 + 1  # dev: unsafe values x[0]
+    assert unsafe_div(x[1] * 10**18, x[0]) > 10**14 - 1  # dev: unsafe values x[i] (input)
+
+    S: uint256 = unsafe_add(x[0], x[1])  # can unsafe add here because we checked x[0] bounds
 
     D: uint256 = 0
     if K0_prev == 0:
@@ -346,11 +395,12 @@ def newton_D(ANN: uint256, gamma: uint256, x_unsorted: uint256[N_COINS], K0_prev
             D = S
 
     __g1k0: uint256 = gamma + 10**18
+    diff: uint256 = 0
 
     for i in range(255):
         D_prev: uint256 = D
         assert D > 0
-        # Unsafe ivision by D is now safe
+        # Unsafe division by D and D_prev is now safe
 
         # K0: uint256 = 10**18
         # for _x in x:
@@ -360,9 +410,9 @@ def newton_D(ANN: uint256, gamma: uint256, x_unsorted: uint256[N_COINS], K0_prev
 
         _g1k0: uint256 = __g1k0
         if _g1k0 > K0:
-            _g1k0 = unsafe_sub(_g1k0, K0) + 1  # > 0
+            _g1k0 = unsafe_add(unsafe_sub(_g1k0, K0), 1)  # > 0
         else:
-            _g1k0 = unsafe_sub(K0, _g1k0) + 1  # > 0
+            _g1k0 = unsafe_add(unsafe_sub(K0, _g1k0), 1)  # > 0
 
         # D / (A * N**N) * _g1k0**2 / gamma**2
         mul1: uint256 = unsafe_div(unsafe_div(unsafe_div(10**18 * D, gamma) * _g1k0, gamma) * _g1k0 * A_MULTIPLIER, ANN)
@@ -370,30 +420,32 @@ def newton_D(ANN: uint256, gamma: uint256, x_unsorted: uint256[N_COINS], K0_prev
         # 2*N*K0 / _g1k0
         mul2: uint256 = unsafe_div(((2 * 10**18) * N_COINS) * K0, _g1k0)
 
+        # calculate neg_fprime. here K0 > 0 is being validated (safediv).
         neg_fprime: uint256 = (S + unsafe_div(S * mul2, 10**18)) + mul1 * N_COINS / K0 - unsafe_div(mul2 * D, 10**18)
 
-        # D -= f / fprime
+        # D -= f / fprime; neg_fprime safediv being validated
         D_plus: uint256 = D * (neg_fprime + S) / neg_fprime
-        D_minus: uint256 = D*D / neg_fprime
+        D_minus: uint256 = unsafe_div(D * D,  neg_fprime)
         if 10**18 > K0:
-            D_minus += unsafe_div(D * (mul1 / neg_fprime), 10**18) * unsafe_sub(10**18, K0) / K0
+            D_minus += unsafe_div(unsafe_div(D * unsafe_div(mul1, neg_fprime), 10**18) * unsafe_sub(10**18, K0), K0)
         else:
-            D_minus -= unsafe_div(D * (mul1 / neg_fprime), 10**18) * unsafe_sub(K0, 10**18) / K0
+            D_minus -= unsafe_div(unsafe_div(D * unsafe_div(mul1, neg_fprime), 10**18) * unsafe_sub(K0, 10**18), K0)
 
         if D_plus > D_minus:
             D = unsafe_sub(D_plus, D_minus)
         else:
             D = unsafe_div(unsafe_sub(D_minus, D_plus), 2)
 
-        diff: uint256 = 0
         if D > D_prev:
             diff = unsafe_sub(D, D_prev)
         else:
             diff = unsafe_sub(D_prev, D)
+
         if diff * 10**14 < max(10**16, D):  # Could reduce precision for gas efficiency here
-            # Test that we are safe with the next newton_y
+
             for _x in x:
                 frac: uint256 = _x * 10**18 / D
+                assert (frac >= 10**16 - 1) and (frac < 10**20 + 1)  # dev: unsafe values x[i]
             return D
 
     raise "Did not converge"
@@ -444,9 +496,76 @@ def get_p(
     denominator: uint256 = (GK0 + unsafe_div(unsafe_div(NNAG2 * _xp[0], _D) * K0, 10**36) )
 
     # p_xy = x * (GK0 + NNAG2 * y / D * K0 / 10**36) / y * 10**18 / denominator
-    # p_xz = x * (GK0 + NNAG2 * z / D * K0 / 10**36) / z * 10**18 / denominator
     # p is in 10**18 precision.
     return unsafe_div(
-               _xp[0] * ( GK0 + unsafe_div(unsafe_div(NNAG2 * _xp[1], _D) * K0, 10**36) ) / _xp[1] * 10**18,
-               denominator
-           )
+        _xp[0] * ( GK0 + unsafe_div(unsafe_div(NNAG2 * _xp[1], _D) * K0, 10**36) ) / _xp[1] * 10**18,
+        denominator
+    )
+
+
+@external
+@pure
+def wad_exp(x: int256) -> int256:
+    """
+    @dev Calculates the natural exponential function of a signed integer with
+         a precision of 1e18.
+    @notice Note that this function consumes about 810 gas units. The implementation
+            is inspired by Remco Bloemen's implementation under the MIT license here:
+            https://xn--2-umb.com/22/exp-ln.
+    @param x The 32-byte variable.
+    @return int256 The 32-byte calculation result.
+    """
+    value: int256 = x
+
+    # If the result is `< 0.5`, we return zero. This happens when we have the following:
+    # "x <= floor(log(0.5e18) * 1e18) ~ -42e18".
+    if (x <= -42_139_678_854_452_767_551):
+        return empty(int256)
+
+    # When the result is "> (2 ** 255 - 1) / 1e18" we cannot represent it as a signed integer.
+    # This happens when "x >= floor(log((2 ** 255 - 1) / 1e18) * 1e18) ~ 135".
+    assert x < 135_305_999_368_893_231_589, "Math: wad_exp overflow"
+
+    # `x` is now in the range "(-42, 136) * 1e18". Convert to "(-42, 136) * 2 ** 96" for higher
+    # intermediate precision and a binary base. This base conversion is a multiplication with
+    # "1e18 / 2 ** 96 = 5 ** 18 / 2 ** 78".
+    value = unsafe_div(x << 78, 5 ** 18)
+
+    # Reduce the range of `x` to "(-½ ln 2, ½ ln 2) * 2 ** 96" by factoring out powers of two
+    # so that "exp(x) = exp(x') * 2 ** k", where `k` is a signer integer. Solving this gives
+    # "k = round(x / log(2))" and "x' = x - k * log(2)". Thus, `k` is in the range "[-61, 195]".
+    k: int256 = unsafe_add(unsafe_div(value << 96, 54_916_777_467_707_473_351_141_471_128), 2 ** 95) >> 96
+    value = unsafe_sub(value, unsafe_mul(k, 54_916_777_467_707_473_351_141_471_128))
+
+    # Evaluate using a "(6, 7)"-term rational approximation. Since `p` is monic,
+    # we will multiply by a scaling factor later.
+    y: int256 = unsafe_add(unsafe_mul(unsafe_add(value, 1_346_386_616_545_796_478_920_950_773_328), value) >> 96, 57_155_421_227_552_351_082_224_309_758_442)
+    p: int256 = unsafe_add(unsafe_mul(unsafe_add(unsafe_mul(unsafe_sub(unsafe_add(y, value), 94_201_549_194_550_492_254_356_042_504_812), y) >> 96,\
+                           28_719_021_644_029_726_153_956_944_680_412_240), value), 4_385_272_521_454_847_904_659_076_985_693_276 << 96)
+
+    # We leave `p` in the "2 ** 192" base so that we do not have to scale it up
+    # again for the division.
+    q: int256 = unsafe_add(unsafe_mul(unsafe_sub(value, 2_855_989_394_907_223_263_936_484_059_900), value) >> 96, 50_020_603_652_535_783_019_961_831_881_945)
+    q = unsafe_sub(unsafe_mul(q, value) >> 96, 533_845_033_583_426_703_283_633_433_725_380)
+    q = unsafe_add(unsafe_mul(q, value) >> 96, 3_604_857_256_930_695_427_073_651_918_091_429)
+    q = unsafe_sub(unsafe_mul(q, value) >> 96, 14_423_608_567_350_463_180_887_372_962_807_573)
+    q = unsafe_add(unsafe_mul(q, value) >> 96, 26_449_188_498_355_588_339_934_803_723_976_023)
+
+    # The polynomial `q` has no zeros in the range because all its roots are complex.
+    # No scaling is required, as `p` is already "2 ** 96" too large. Also,
+    # `r` is in the range "(0.09, 0.25) * 2**96" after the division.
+    r: int256 = unsafe_div(p, q)
+
+    # To finalise the calculation, we have to multiply `r` by:
+    #   - the scale factor "s = ~6.031367120",
+    #   - the factor "2 ** k" from the range reduction, and
+    #   - the factor "1e18 / 2 ** 96" for the base conversion.
+    # We do this all at once, with an intermediate result in "2**213" base,
+    # so that the final right shift always gives a positive value.
+
+    # Note that to circumvent Vyper's safecast feature for the potentially
+    # negative parameter value `r`, we first convert `r` to `bytes32` and
+    # subsequently to `uint256`. Remember that the EVM default behaviour is
+    # to use two's complement representation to handle signed integers.
+    return convert(unsafe_mul(convert(convert(r, bytes32), uint256), 3_822_833_074_963_236_453_042_738_258_902_158_003_155_416_615_667) >>\
+           convert(unsafe_sub(195, k), uint256), int256)
diff --git a/contracts/main/CurveCryptoViews2Optimized.vy b/contracts/main/CurveCryptoViews2Optimized.vy
index 739a3f5d..90c62640 100644
--- a/contracts/main/CurveCryptoViews2Optimized.vy
+++ b/contracts/main/CurveCryptoViews2Optimized.vy
@@ -1,4 +1,4 @@
-# @version 0.3.10
+# pragma version 0.3.10
 """
 @title CurveCryptoViews2Optimized
 @author Curve.Fi
@@ -43,8 +43,12 @@ interface Math:
         D: uint256,
         i: uint256,
     ) -> uint256[2]: view
-    def reduction_coefficient(
-        x: uint256[N_COINS], fee_gamma: uint256
+    def newton_y(
+        ANN: uint256,
+        gamma: uint256,
+        x: uint256[N_COINS],
+        D: uint256,
+        i: uint256,
     ) -> uint256: view
 
 
@@ -162,14 +166,11 @@ def _calc_D_ramp(
 ) -> uint256:
 
     math: Math = Curve(swap).MATH()
-
     D: uint256 = Curve(swap).D()
     if Curve(swap).future_A_gamma_time() > block.timestamp:
         _xp: uint256[N_COINS] = xp
         _xp[0] *= precisions[0]
-        _xp[1] = (
-            _xp[1] * price_scale * precisions[1] / PRECISION
-        )
+        _xp[1] = _xp[1] * price_scale * precisions[1] / PRECISION
         D = math.newton_D(A, gamma, _xp, 0)
 
     return D
@@ -201,16 +202,15 @@ def _get_dx_fee(
     # adjust xp with output dy. dy contains fee element, which we handle later
     # (hence this internal method is called _get_dx_fee)
     xp[j] -= dy
-    xp = [xp[0] * precisions[0], xp[1] * price_scale / PRECISION]
+    xp = [xp[0] * precisions[0], xp[1] * price_scale * precisions[1] / PRECISION]
 
     x_out: uint256[2] = math.get_y(A, gamma, xp, D, i)
     dx: uint256 = x_out[0] - xp[i]
     xp[i] = x_out[0]
 
     if i > 0:
-        dx = dy * PRECISION / price_scale
-    else:
-        dx /= precisions[0]
+        dx = dx * PRECISION / price_scale
+    dx /= precisions[i]
 
     return dx, xp
 
@@ -238,15 +238,18 @@ def _get_dy_nofee(
 
     # adjust xp with input dx
     xp[i] += dx
-    xp = [xp[0] * precisions[0], xp[1] * price_scale / PRECISION]
+    xp = [
+        xp[0] * precisions[0],
+        xp[1] * price_scale * precisions[1] / PRECISION
+    ]
 
     y_out: uint256[2] = math.get_y(A, gamma, xp, D, j)
+
     dy: uint256 = xp[j] - y_out[0] - 1
     xp[j] = y_out[0]
     if j > 0:
         dy = dy * PRECISION / price_scale
-    else:
-        dy /= precisions[0]
+    dy /= precisions[j]
 
     return dy, xp
 
@@ -277,15 +280,13 @@ def _calc_dtoken_nofee(
         for k in range(N_COINS):
             xp[k] -= amounts[k]
 
-    amountsp[0] *= precisions[0]
-
     xp = [
         xp[0] * precisions[0],
-        xp[1] * price_scale / PRECISION
+        xp[1] * price_scale * precisions[1] / PRECISION
     ]
     amountsp = [
-        amountsp[0] * precisions[0],
-        amountsp[1] * price_scale / PRECISION
+        amountsp[0]* precisions[0],
+        amountsp[1] * price_scale * precisions[1] / PRECISION
     ]
 
     D: uint256 = math.newton_D(A, gamma, xp, 0)
@@ -314,7 +315,6 @@ def _calc_withdraw_one_coin(
     math: Math = Curve(swap).MATH()
 
     xx: uint256[N_COINS] = empty(uint256[N_COINS])
-    price_scale: uint256 = Curve(swap).price_scale()
     for k in range(N_COINS):
         xx[k] = Curve(swap).balances(k)
 
@@ -324,7 +324,7 @@ def _calc_withdraw_one_coin(
     D0: uint256 = 0
     p: uint256 = 0
 
-    price_scale_i: uint256 = PRECISION * precisions[0]
+    price_scale_i: uint256 = Curve(swap).price_scale() * precisions[1]
     xp: uint256[N_COINS] = [
         xx[0] * precisions[0],
         unsafe_div(xx[1] * price_scale_i, PRECISION)
@@ -357,10 +357,15 @@ def _calc_withdraw_one_coin(
 @internal
 @view
 def _fee(xp: uint256[N_COINS], swap: address) -> uint256:
-    math: Math = Curve(swap).MATH()
+
     packed_fee_params: uint256 = Curve(swap).packed_fee_params()
-    fee_params: uint256[3] = self._unpack(packed_fee_params)
-    f: uint256 = math.reduction_coefficient(xp, fee_params[2])
+    fee_params: uint256[3] = self._unpack_3(packed_fee_params)
+    f: uint256 = xp[0] + xp[1]
+    f = fee_params[2] * 10**18 / (
+        fee_params[2] + 10**18 -
+        (10**18 * N_COINS**N_COINS) * xp[0] / f * xp[1] / f
+    )
+
     return (fee_params[0] * f + fee_params[1] * (10**18 - f)) / 10**18
 
 
@@ -395,7 +400,7 @@ def _prep_calc(swap: address) -> (
 
 @internal
 @view
-def _unpack(_packed: uint256) -> uint256[3]:
+def _unpack_3(_packed: uint256) -> uint256[3]:
     """
     @notice Unpacks a uint256 into 3 integers (values must be <= 10**18)
     @param val The uint256 to unpack
diff --git a/contracts/main/CurveTwocryptoFactory.vy b/contracts/main/CurveTwocryptoFactory.vy
index d55e16a5..e125673c 100644
--- a/contracts/main/CurveTwocryptoFactory.vy
+++ b/contracts/main/CurveTwocryptoFactory.vy
@@ -1,4 +1,5 @@
-# @version 0.3.10
+# pragma version 0.3.10
+# pragma optimize codesize
 """
 @title CurveTwocryptoFactory
 @author Curve.Fi
@@ -71,12 +72,6 @@ A_MULTIPLIER: constant(uint256) = 10000
 # Limits
 MAX_FEE: constant(uint256) = 10 * 10 ** 9
 
-MIN_GAMMA: constant(uint256) = 10 ** 10
-MAX_GAMMA: constant(uint256) = 5 * 10**16
-
-MIN_A: constant(uint256) = N_COINS ** N_COINS * A_MULTIPLIER / 100
-MAX_A: constant(uint256) = 1000 * A_MULTIPLIER * N_COINS**N_COINS
-
 admin: public(address)
 future_admin: public(address)
 
@@ -91,12 +86,9 @@ math_implementation: public(address)
 # mapping of coins -> pools for trading
 # a mapping key is generated for each pair of addresses via
 # `bitwise_xor(convert(a, uint256), convert(b, uint256))`
-markets: HashMap[uint256, address[4294967296]]
-market_counts: HashMap[uint256, uint256]
-
-pool_count: public(uint256)              # actual length of pool_list
+markets: HashMap[uint256, DynArray[address, 4294967296]]
 pool_data: HashMap[address, PoolArray]
-pool_list: public(address[4294967296])   # master list of pools
+pool_list: public(DynArray[address, 4294967296])   # master list of pools
 
 
 @external
@@ -110,8 +102,8 @@ def __init__(_fee_receiver: address, _admin: address):
 
 
 @internal
-@view
-def _pack(x: uint256[3]) -> uint256:
+@pure
+def _pack_3(x: uint256[3]) -> uint256:
     """
     @notice Packs 3 integers with values <= 10**18 into a uint256
     @param x The uint256[3] to pack
@@ -120,6 +112,11 @@ def _pack(x: uint256[3]) -> uint256:
     return (x[0] << 128) | (x[1] << 64) | x[2]
 
 
+@pure
+@internal
+def _pack_2(p1: uint256, p2: uint256) -> uint256:
+    return p1 | (p2 << 128)
+
 
 # <--- Pool Deployers --->
 
@@ -147,14 +144,9 @@ def deploy_pool(
     @return Address of the deployed pool
     """
     pool_implementation: address = self.pool_implementations[implementation_id]
+    _math_implementation: address = self.math_implementation
     assert pool_implementation != empty(address), "Pool implementation not set"
-
-    # Validate parameters
-    assert A > MIN_A-1
-    assert A < MAX_A+1
-
-    assert gamma > MIN_GAMMA-1
-    assert gamma < MAX_GAMMA+1
+    assert _math_implementation != empty(address), "Math implementation not set"
 
     assert mid_fee < MAX_FEE-1  # mid_fee can be zero
     assert out_fee >= mid_fee
@@ -180,44 +172,44 @@ def deploy_pool(
         d: uint256 = ERC20(_coins[i]).decimals()
         assert d < 19, "Max 18 decimals for coins"
         decimals[i] = d
-        precisions[i] = 10** (18 - d)
+        precisions[i] = 10 ** (18 - d)
+
+    # pack precision
+    packed_precisions: uint256 = self._pack_2(precisions[0], precisions[1])
 
     # pack fees
-    packed_fee_params: uint256 = self._pack(
+    packed_fee_params: uint256 = self._pack_3(
         [mid_fee, out_fee, fee_gamma]
     )
 
     # pack liquidity rebalancing params
-    packed_rebalancing_params: uint256 = self._pack(
+    packed_rebalancing_params: uint256 = self._pack_3(
         [allowed_extra_profit, adjustment_step, ma_exp_time]
     )
 
-    # pack A_gamma
-    packed_A_gamma: uint256 = A << 128
-    packed_A_gamma = packed_A_gamma | gamma
+    # pack gamma and A
+    packed_gamma_A: uint256 = self._pack_2(gamma, A)
 
     # pool is an ERC20 implementation
     _salt: bytes32 = block.prevhash
-    _math_implementation: address = self.math_implementation
     pool: address = create_from_blueprint(
-        pool_implementation,
-        _name,
-        _symbol,
-        _coins,
-        _math_implementation,
-        _salt,
-        precisions,
-        packed_A_gamma,
-        packed_fee_params,
-        packed_rebalancing_params,
-        initial_price,
+        pool_implementation,  # blueprint: address
+        _name,  # String[64]
+        _symbol,  # String[32]
+        _coins,  # address[N_COINS]
+        _math_implementation,  # address
+        _salt,  # bytes32
+        packed_precisions,  # uint256
+        packed_gamma_A,  # uint256
+        packed_fee_params,  # uint256
+        packed_rebalancing_params,  # uint256
+        initial_price,  # uint256
         code_offset=3
     )
 
     # populate pool data
-    length: uint256 = self.pool_count
-    self.pool_list[length] = pool
-    self.pool_count = length + 1
+    self.pool_list.append(pool)
+
     self.pool_data[pool].decimals = decimals
     self.pool_data[pool].coins = _coins
     self.pool_data[pool].implementation = pool_implementation
@@ -233,7 +225,7 @@ def deploy_pool(
         _math_implementation,
         _salt,
         precisions,
-        packed_A_gamma,
+        packed_gamma_A,
         packed_fee_params,
         packed_rebalancing_params,
         initial_price,
@@ -249,10 +241,7 @@ def _add_coins_to_market(coin_a: address, coin_b: address, pool: address):
     key: uint256 = (
         convert(coin_a, uint256) ^ convert(coin_b, uint256)
     )
-
-    length: uint256 = self.market_counts[key]
-    self.markets[key][length] = pool
-    self.market_counts[key] = length + 1
+    self.markets[key].append(pool)
 
 
 @external
@@ -390,6 +379,16 @@ def find_pool_for_coins(_from: address, _to: address, i: uint256 = 0) -> address
 # <--- Pool Getters --->
 
 
+@view
+@external
+def pool_count() -> uint256:
+    """
+    @notice Get number of pools deployed from the factory
+    @return Number of pools deployed from factory
+    """
+    return len(self.pool_list)
+
+
 @view
 @external
 def get_coins(_pool: address) -> address[N_COINS]:
@@ -472,4 +471,4 @@ def get_market_counts(coin_a: address, coin_b: address) -> uint256:
         convert(coin_a, uint256) ^ convert(coin_b, uint256)
     )
 
-    return self.market_counts[key]
+    return len(self.markets[key])
diff --git a/contracts/main/CurveTwocryptoFactoryHandler.vy b/contracts/main/CurveTwocryptoFactoryHandler.vy
index f2e0fd05..ffb200cc 100644
--- a/contracts/main/CurveTwocryptoFactoryHandler.vy
+++ b/contracts/main/CurveTwocryptoFactoryHandler.vy
@@ -1,4 +1,4 @@
-# @version 0.3.10
+# pragma version 0.3.10
 """
 @title CurveTwocryptoFactoryHandler
 @author Curve.Fi
@@ -20,7 +20,7 @@ interface BaseRegistry:
 
 interface CurvePool:
     def adjustment_step() -> uint256: view
-    def admin_fee() -> uint256: view
+    def ADMIN_FEE() -> uint256: view
     def allowed_extra_profit() -> uint256: view
     def A() -> uint256: view
     def balances(i: uint256) -> uint256: view
@@ -29,7 +29,7 @@ interface CurvePool:
     def fee_gamma() -> uint256: view
     def gamma() -> uint256: view
     def get_virtual_price() -> uint256: view
-    def ma_half_time() -> uint256: view
+    def ma_time() -> uint256: view
     def mid_fee() -> uint256: view
     def out_fee() -> uint256: view
     def virtual_price() -> uint256: view
@@ -119,7 +119,7 @@ def _get_gauge_type(_gauge: address) -> int128:
     success, response = raw_call(
         GAUGE_CONTROLLER,
         concat(
-            method_id("gauge_type(address)"),
+            method_id("gauge_types(address)"),
             convert(_gauge, bytes32),
         ),
         max_outsize=32,
@@ -170,7 +170,7 @@ def get_admin_balances(_pool: address) -> uint256[MAX_METAREGISTRY_COINS]:
 
     xcp_profit: uint256 = CurvePool(_pool).xcp_profit()
     xcp_profit_a: uint256 = CurvePool(_pool).xcp_profit_a()
-    admin_fee: uint256 = CurvePool(_pool).admin_fee()
+    admin_fee: uint256 = CurvePool(_pool).ADMIN_FEE()
     admin_balances: uint256[MAX_METAREGISTRY_COINS] = empty(uint256[MAX_METAREGISTRY_COINS])
 
     # admin balances are non zero if pool has made more than allowed profits:
@@ -275,7 +275,7 @@ def get_fees(_pool: address) -> uint256[10]:
     fees: uint256[10] = empty(uint256[10])
     pool_fees: uint256[4] = [
         CurvePool(_pool).fee(),
-        CurvePool(_pool).admin_fee(),
+        CurvePool(_pool).ADMIN_FEE(),
         CurvePool(_pool).mid_fee(),
         CurvePool(_pool).out_fee()
     ]
@@ -390,7 +390,7 @@ def get_pool_params(_pool: address) -> uint256[20]:
     pool_params[3] = CurvePool(_pool).allowed_extra_profit()
     pool_params[4] = CurvePool(_pool).fee_gamma()
     pool_params[5] = CurvePool(_pool).adjustment_step()
-    pool_params[6] = CurvePool(_pool).ma_half_time()
+    pool_params[6] = CurvePool(_pool).ma_time()
     return pool_params
 
 
diff --git a/contracts/main/CurveTwocryptoOptimized.vy b/contracts/main/CurveTwocryptoOptimized.vy
index 7cb52444..d036988d 100644
--- a/contracts/main/CurveTwocryptoOptimized.vy
+++ b/contracts/main/CurveTwocryptoOptimized.vy
@@ -1,5 +1,5 @@
-# @version 0.3.10
-#pragma optimize gas
+# pragma version 0.3.10
+# pragma optimize gas
 """
 @title CurveTwocryptoOptimized
 @author Curve.Fi
@@ -123,7 +123,6 @@ event ClaimAdminFee:
 
 N_COINS: constant(uint256) = 2
 PRECISION: constant(uint256) = 10**18  # <------- The precision to convert to.
-A_MULTIPLIER: constant(uint256) = 10000
 PRECISIONS: immutable(uint256[N_COINS])
 
 MATH: public(immutable(Math))
@@ -134,8 +133,8 @@ cached_price_scale: uint256  # <------------------------ Internal price scale.
 cached_price_oracle: uint256  # <------- Price target given by moving average.
 cached_xcp_oracle: uint256  # <----------- EMA of totalSupply * virtual_price.
 
-last_prices: uint256
-last_prices_timestamp: public(uint256)
+last_prices: public(uint256)
+last_timestamp: public(uint256)    # idx 0 is for prices, idx 1 is for xcp.
 last_xcp: public(uint256)
 xcp_ma_time: public(uint256)
 
@@ -173,12 +172,14 @@ NOISE_FEE: constant(uint256) = 10**5  # <---------------------------- 0.1 BPS.
 # ----------------------- Admin params ---------------------------------------
 
 last_admin_fee_claim_timestamp: uint256
+admin_lp_virtual_balance: uint256
 
 MIN_RAMP_TIME: constant(uint256) = 86400
 MIN_ADMIN_FEE_CLAIM_INTERVAL: constant(uint256) = 86400
 
+A_MULTIPLIER: constant(uint256) = 10000
 MIN_A: constant(uint256) = N_COINS**N_COINS * A_MULTIPLIER / 10
-MAX_A: constant(uint256) = N_COINS**N_COINS * A_MULTIPLIER * 100000
+MAX_A: constant(uint256) = N_COINS**N_COINS * A_MULTIPLIER * 1000
 MAX_A_CHANGE: constant(uint256) = 10
 MIN_GAMMA: constant(uint256) = 10**10
 MAX_GAMMA: constant(uint256) = 5 * 10**16
@@ -217,8 +218,8 @@ def __init__(
     _coins: address[N_COINS],
     _math: address,
     _salt: bytes32,
-    precisions: uint256[N_COINS],
-    packed_A_gamma: uint256,
+    packed_precisions: uint256,
+    packed_gamma_A: uint256,
     packed_fee_params: uint256,
     packed_rebalancing_params: uint256,
     initial_price: uint256,
@@ -231,10 +232,20 @@ def __init__(
     symbol = _symbol
     coins = _coins
 
-    PRECISIONS = precisions  # <------------------------ Precisions of coins.
+    PRECISIONS = self._unpack_2(packed_precisions)  # <-- Precisions of coins.
+
+    # --------------- Validate A and gamma parameters here and not in factory.
+    gamma_A: uint256[2] = self._unpack_2(packed_gamma_A)  # gamma is at idx 0.
 
-    self.initial_A_gamma = packed_A_gamma  # <------------------- A and gamma.
-    self.future_A_gamma = packed_A_gamma
+    assert gamma_A[0] > MIN_GAMMA-1
+    assert gamma_A[0] < MAX_GAMMA+1
+
+    assert gamma_A[1] > MIN_A-1
+    assert gamma_A[1] < MAX_A+1
+
+    self.initial_A_gamma = packed_gamma_A
+    self.future_A_gamma = packed_gamma_A
+    # ------------------------------------------------------------------------
 
     self.packed_rebalancing_params = packed_rebalancing_params  # <-- Contains
     #               rebalancing params: allowed_extra_profit, adjustment_step,
@@ -246,7 +257,7 @@ def __init__(
     self.cached_price_scale = initial_price
     self.cached_price_oracle = initial_price
     self.last_prices = initial_price
-    self.last_prices_timestamp = block.timestamp
+    self.last_timestamp = self._pack_2(block.timestamp, block.timestamp)
     self.xcp_profit_a = 10**18
     self.xcp_ma_time = 62324  # <--------- 12 hours default on contract start.
 
@@ -370,16 +381,31 @@ def exchange(
     @param receiver Address to send the output coin to. Default is msg.sender
     @return uint256 Amount of tokens at index j received by the `receiver
     """
-    return self._exchange(
+    # _transfer_in updates self.balances here:
+    dx_received: uint256 = self._transfer_in(
+        i,
+        dx,
         msg.sender,
+        False
+    )
+
+    # No ERC20 token transfers occur here:
+    out: uint256[3] = self._exchange(
         i,
         j,
-        dx,
+        dx_received,
         min_dy,
-        receiver,
-        False,
     )
 
+    # _transfer_out updates self.balances here. Update to state occurs before
+    # external calls:
+    self._transfer_out(j, out[0], receiver)
+
+    # log:
+    log TokenExchange(msg.sender, i, dx_received, j, out[0], out[1], out[2])
+
+    return out[0]
+
 
 @external
 @nonreentrant('lock')
@@ -404,16 +430,31 @@ def exchange_received(
     @param receiver Address to send the output coin to
     @return uint256 Amount of tokens at index j received by the `receiver`
     """
-    return self._exchange(
+    # _transfer_in updates self.balances here:
+    dx_received: uint256 = self._transfer_in(
+        i,
+        dx,
         msg.sender,
+        True  # <---- expect_optimistic_transfer is set to True here.
+    )
+
+    # No ERC20 token transfers occur here:
+    out: uint256[3] = self._exchange(
         i,
         j,
-        dx,
+        dx_received,
         min_dy,
-        receiver,
-        True,
     )
 
+    # _transfer_out updates self.balances here. Update to state occurs before
+    # external calls:
+    self._transfer_out(j, out[0], receiver)
+
+    # log:
+    log TokenExchange(msg.sender, i, dx_received, j, out[0], out[1], out[2])
+
+    return out[0]
+
 
 @external
 @nonreentrant("lock")
@@ -462,14 +503,14 @@ def add_liquidity(
 
     xp = [
         xp[0] * PRECISIONS[0],
-        xp[1] * price_scale / PRECISION
+        unsafe_div(xp[1] * price_scale * PRECISIONS[1], PRECISION)
     ]
     xp_old = [
         xp_old[0] * PRECISIONS[0],
-        xp_old[1] * price_scale / PRECISION
+        unsafe_div(xp_old[1] * price_scale * PRECISIONS[1], PRECISION)
     ]
 
-    for i in range(N_COINS):         # TODO: optimize
+    for i in range(N_COINS):
         if amounts_received[i] > 0:
             amountsp[i] = xp[i] - xp_old[i]
 
@@ -503,6 +544,7 @@ def add_liquidity(
         d_token -= d_token_fee
         token_supply += d_token
         self.mint(receiver, d_token)
+        self.admin_lp_virtual_balance += unsafe_div(ADMIN_FEE * d_token_fee, 10**10)
 
         price_scale = self.tweak_price(A_gamma, xp, D, 0)
 
@@ -532,8 +574,6 @@ def add_liquidity(
         price_scale
     )
 
-    self._claim_admin_fees()  # <--------- Auto-claim admin fees occasionally.
-
     return d_token
 
 
@@ -600,6 +640,38 @@ def remove_liquidity(
 
     log RemoveLiquidity(msg.sender, withdraw_amounts, total_supply - _amount)
 
+    # --------------------------- Upkeep xcp oracle --------------------------
+
+    # Update xcp since liquidity was removed:
+    xp: uint256[N_COINS] = self.xp(self.balances, self.cached_price_scale)
+    last_xcp: uint256 = isqrt(xp[0] * xp[1])  # <----------- Cache it for now.
+
+    last_timestamp: uint256[2] = self._unpack_2(self.last_timestamp)
+    if last_timestamp[1] < block.timestamp:
+
+        cached_xcp_oracle: uint256 = self.cached_xcp_oracle
+        alpha: uint256 = MATH.wad_exp(
+            -convert(
+                unsafe_div(
+                    unsafe_sub(block.timestamp, last_timestamp[1]) * 10**18,
+                    self.xcp_ma_time  # <---------- xcp ma time has is longer.
+                ),
+                int256,
+            )
+        )
+
+        self.cached_xcp_oracle = unsafe_div(
+            last_xcp * (10**18 - alpha) + cached_xcp_oracle * alpha,
+            10**18
+        )
+        last_timestamp[1] = block.timestamp
+
+        # Pack and store timestamps:
+        self.last_timestamp = self._pack_2(last_timestamp[0], last_timestamp[1])
+
+    # Store last xcp
+    self.last_xcp = last_xcp
+
     return withdraw_amounts
 
 
@@ -622,6 +694,8 @@ def remove_liquidity_one_coin(
     @return Amount of tokens at index i received by the `receiver`
     """
 
+    self._claim_admin_fees()  # <--------- Auto-claim admin fees occasionally.
+
     A_gamma: uint256[2] = self._A_gamma()
 
     dy: uint256 = 0
@@ -659,8 +733,6 @@ def remove_liquidity_one_coin(
         msg.sender, token_amount, i, dy, approx_fee, packed_price_scale
     )
 
-    self._claim_admin_fees()  # <--------- Auto-claim admin fees occasionally.
-
     return dy
 
 
@@ -669,7 +741,7 @@ def remove_liquidity_one_coin(
 
 @internal
 @pure
-def _pack(x: uint256[3]) -> uint256:
+def _pack_3(x: uint256[3]) -> uint256:
     """
     @notice Packs 3 integers with values <= 10**18 into a uint256
     @param x The uint256[3] to pack
@@ -680,7 +752,7 @@ def _pack(x: uint256[3]) -> uint256:
 
 @internal
 @pure
-def _unpack(_packed: uint256) -> uint256[3]:
+def _unpack_3(_packed: uint256) -> uint256[3]:
     """
     @notice Unpacks a uint256 into 3 integers (values must be <= 10**18)
     @param val The uint256 to unpack
@@ -693,44 +765,40 @@ def _unpack(_packed: uint256) -> uint256[3]:
     ]
 
 
+@pure
+@internal
+def _pack_2(p1: uint256, p2: uint256) -> uint256:
+    return p1 | (p2 << 128)
+
+
+@pure
+@internal
+def _unpack_2(packed: uint256) -> uint256[2]:
+    return [packed & (2**128 - 1), packed >> 128]
+
+
 # ---------------------- AMM Internal Functions -------------------------------
 
 
 @internal
 def _exchange(
-    sender: address,
     i: uint256,
     j: uint256,
-    dx: uint256,
+    dx_received: uint256,
     min_dy: uint256,
-    receiver: address,
-    expect_optimistic_transfer: bool,
-) -> uint256:
+) -> uint256[3]:
 
     assert i != j  # dev: coin index out of range
-    assert dx > 0  # dev: do not exchange 0 coins
+    assert dx_received > 0  # dev: do not exchange 0 coins
 
     A_gamma: uint256[2] = self._A_gamma()
     xp: uint256[N_COINS] = self.balances
     dy: uint256 = 0
 
-    y: uint256 = xp[j]  # <----------------- if j > N_COINS, this will revert.
-    x0: uint256 = xp[i]  # <--------------- if i > N_COINS, this will  revert.
-
-    ########################## TRANSFER IN <-------
-
-    # _transfer_in updates self.balances here:
-    dx_received: uint256 = self._transfer_in(
-        i,
-        dx,
-        sender,
-        expect_optimistic_transfer  # <---- If True, pool expects dx tokens to
-    )  #                                                    be transferred in.
-
-    xp[i] = x0 + dx_received
+    y: uint256 = xp[j]
+    x0: uint256 = xp[i] - dx_received  # old xp[i]
 
     price_scale: uint256 = self.cached_price_scale
-
     xp = [
         xp[0] * PRECISIONS[0],
         unsafe_div(xp[1] * price_scale * PRECISIONS[1], PRECISION)
@@ -756,18 +824,16 @@ def _exchange(
     D: uint256 = self.D
     y_out: uint256[2] = MATH.get_y(A_gamma[0], A_gamma[1], xp, D, j)
     dy = xp[j] - y_out[0]
-    xp[j] = xp[j] - dy
-    dy = dy - 1
+    xp[j] -= dy
+    dy -= 1
 
     if j > 0:
         dy = dy * PRECISION / price_scale
-    dy = dy / PRECISIONS[j]
+    dy /= PRECISIONS[j]
 
     fee: uint256 = unsafe_div(self._fee(xp) * dy, 10**10)
-
     dy -= fee  # <--------------------- Subtract fee from the outgoing amount.
     assert dy >= min_dy, "Slippage"
-
     y -= dy
 
     y *= PRECISIONS[j]
@@ -779,17 +845,7 @@ def _exchange(
 
     price_scale = self.tweak_price(A_gamma, xp, 0, y_out[1])
 
-    # --------------------------- Do Transfers out ---------------------------
-
-    ########################## -------> TRANSFER OUT
-
-    # _transfer_out updates self.balances here. Update to state occurs before
-    # external calls:
-    self._transfer_out(j, dy, receiver)
-
-    log TokenExchange(sender, i, dx_received, j, dy, fee, price_scale)
-
-    return dy
+    return [dy, fee, price_scale]
 
 
 @internal
@@ -816,7 +872,7 @@ def tweak_price(
     price_oracle: uint256 = self.cached_price_oracle
     last_prices: uint256 = self.last_prices
     price_scale: uint256 = self.cached_price_scale
-    rebalancing_params: uint256[3] = self._unpack(self.packed_rebalancing_params)
+    rebalancing_params: uint256[3] = self._unpack_3(self.packed_rebalancing_params)
     # Contains: allowed_extra_profit, adjustment_step, ma_time. -----^
 
     total_supply: uint256 = self.totalSupply
@@ -825,8 +881,9 @@ def tweak_price(
 
     # ----------------------- Update Oracles if needed -----------------------
 
-    last_timestamp: uint256 = self.last_prices_timestamp
-    if last_timestamp < block.timestamp:  # 0th index is for price_oracle.
+    last_timestamp: uint256[2] = self._unpack_2(self.last_timestamp)
+    alpha: uint256 = 0
+    if last_timestamp[0] < block.timestamp:  # 0th index is for price_oracle.
 
         #   The moving average price oracle is calculated using the last_price
         #      of the trade at the previous block, and the price oracle logged
@@ -834,10 +891,10 @@ def tweak_price(
 
         # ------------------ Calculate moving average params -----------------
 
-        alpha: uint256 = MATH.wad_exp(
+        alpha = MATH.wad_exp(
             -convert(
                 unsafe_div(
-                    (block.timestamp - last_timestamp) * 10**18,
+                    unsafe_sub(block.timestamp, last_timestamp[0]) * 10**18,
                     rebalancing_params[2]  # <----------------------- ma_time.
                 ),
                 int256,
@@ -855,14 +912,17 @@ def tweak_price(
         )
 
         self.cached_price_oracle = price_oracle
+        last_timestamp[0] = block.timestamp
+
+    # ----------------------------------------------------- Update xcp oracle.
 
-        # ------------------------------------------------- Update xcp oracle.
+    if last_timestamp[1] < block.timestamp:
 
         cached_xcp_oracle: uint256 = self.cached_xcp_oracle
         alpha = MATH.wad_exp(
             -convert(
                 unsafe_div(
-                    (block.timestamp - last_timestamp) * 10**18,
+                    unsafe_sub(block.timestamp, last_timestamp[1]) * 10**18,
                     self.xcp_ma_time  # <---------- xcp ma time has is longer.
                 ),
                 int256,
@@ -875,7 +935,9 @@ def tweak_price(
         )
 
         # Pack and store timestamps:
-        self.last_prices_timestamp = block.timestamp
+        last_timestamp[1] = block.timestamp
+
+    self.last_timestamp = self._pack_2(last_timestamp[0], last_timestamp[1])
 
     #  `price_oracle` is used further on to calculate its vector distance from
     # price_scale. This distance is used to calculate the amount of adjustment
@@ -907,7 +969,7 @@ def tweak_price(
 
     if old_virtual_price > 0:
 
-        xcp: uint256 = isqrt(xp[0] * xp[1] / 10**18)                               # TODO: Check precision!
+        xcp: uint256 = isqrt(xp[0] * xp[1])
         virtual_price = 10**18 * xcp / total_supply
 
         xcp_profit = unsafe_div(
@@ -977,13 +1039,13 @@ def tweak_price(
             # ------------------------------------- Convert xp to real prices.
             xp = [
                 unsafe_div(D, N_COINS),
-                D * PRECISION / (N_COINS * p_new)                                   #  TODO: can use unsafediv?
+                D * PRECISION / (N_COINS * p_new)
             ]
 
             # ---------- Calculate new virtual_price using new xp and D. Reuse
             #              `old_virtual_price` (but it has new virtual_price).
             old_virtual_price = unsafe_div(
-                10**18 * isqrt(xp[0] * xp[1] / 10**18), total_supply                # TODO: can use unsafemath?
+                10**18 * isqrt(xp[0] * xp[1]), total_supply
             )  # <----- unsafe_div because we did safediv before (if vp>1e18)
 
             # ---------------------------- Proceed if we've got enough profit.
@@ -1024,17 +1086,14 @@ def _claim_admin_fees():
 
     last_claim_time: uint256 = self.last_admin_fee_claim_timestamp
     if (
-        block.timestamp - last_claim_time < MIN_ADMIN_FEE_CLAIM_INTERVAL or
+        unsafe_sub(block.timestamp, last_claim_time) < MIN_ADMIN_FEE_CLAIM_INTERVAL or
         self.future_A_gamma_time > block.timestamp
     ):
         return
 
-    A_gamma: uint256[2] = self._A_gamma()
-
     xcp_profit: uint256 = self.xcp_profit  # <---------- Current pool profits.
     xcp_profit_a: uint256 = self.xcp_profit_a  # <- Profits at previous claim.
     current_lp_token_supply: uint256 = self.totalSupply
-    D: uint256 = self.D
 
     # Do not claim admin fees if:
     # 1. insufficient profits accrued since last claim, and
@@ -1046,6 +1105,8 @@ def _claim_admin_fees():
 
     # ---------- Conditions met to claim admin fees: compute state. ----------
 
+    A_gamma: uint256[2] = self._A_gamma()
+    D: uint256 = self.D
     vprice: uint256 = self.virtual_price
     price_scale: uint256 = self.cached_price_scale
     fee_receiver: address = factory.fee_receiver()
@@ -1066,13 +1127,16 @@ def _claim_admin_fees():
 
     # ------------------------------ Claim admin fees by minting admin's share
     #                                                of the pool in LP tokens.
-    admin_share: uint256 = 0
+
+    # This is the admin fee tokens claimed in self.add_liquidity. We add it to
+    # the LP token share that the admin needs to claim:
+    admin_share: uint256 = self.admin_lp_virtual_balance
     frac: uint256 = 0
     if fee_receiver != empty(address) and fees > 0:
 
         # -------------------------------- Calculate admin share to be minted.
         frac = vprice * 10**18 / (vprice - fees) - 10**18
-        admin_share = current_lp_token_supply * frac / 10**18
+        admin_share += current_lp_token_supply * frac / 10**18
 
         # ------ Subtract fees from profits that will be used for rebalancing.
         xcp_profit -= fees * 2
@@ -1092,6 +1156,9 @@ def _claim_admin_fees():
 
     # ---------------------------- Update State ------------------------------
 
+    # Set admin virtual LP balances to zero because we claimed:
+    self.admin_lp_virtual_balance = 0
+
     self.xcp_profit = xcp_profit
     self.last_admin_fee_claim_timestamp = block.timestamp
 
@@ -1166,7 +1233,7 @@ def _A_gamma() -> uint256[2]:
 @view
 def _fee(xp: uint256[N_COINS]) -> uint256:
 
-    fee_params: uint256[3] = self._unpack(self.packed_fee_params)
+    fee_params: uint256[3] = self._unpack_3(self.packed_fee_params)
     f: uint256 = xp[0] + xp[1]
     f = fee_params[2] * 10**18 / (
         fee_params[2] + 10**18 -
@@ -1188,7 +1255,7 @@ def get_xcp(D: uint256, price_scale: uint256) -> uint256:
         D * PRECISION / (price_scale * N_COINS)
     ]
 
-    return isqrt(x[0] * x[1] / 10**18)  # <------------------- Geometric Mean.       # TODO: Check precision!
+    return isqrt(x[0] * x[1])  # <------------------- Geometric Mean.
 
 
 @view
@@ -1263,7 +1330,7 @@ def _calc_withdraw_one_coin(
 
     xp_imprecise: uint256[N_COINS] = xp
     xp_correction: uint256 = xp[i] * N_COINS * token_amount / token_supply
-    fee: uint256 = self._unpack(self.packed_fee_params)[1]  # <- self.out_fee.
+    fee: uint256 = self._unpack_3(self.packed_fee_params)[1]  # <- self.out_fee.
 
     if xp_correction < xp_imprecise[i]:
         xp_imprecise[i] -= xp_correction
@@ -1449,6 +1516,43 @@ def burnFrom(_to: address, _value: uint256) -> bool:
 # ------------------------- AMM View Functions -------------------------------
 
 
+@internal
+@view
+def internal_price_oracle() -> uint256:
+    """
+    @notice Returns the oracle price of the coin at index `k` w.r.t the coin
+            at index 0.
+    @dev The oracle is an exponential moving average, with a periodicity
+         determined by `self.ma_time`. The aggregated prices are cached state
+         prices (dy/dx) calculated AFTER the latest trade.
+    @param k The index of the coin.
+    @return uint256 Price oracle value of kth coin.
+    """
+    price_oracle: uint256 = self.cached_price_oracle
+    price_scale: uint256 = self.cached_price_scale
+    last_prices_timestamp: uint256 = self._unpack_2(self.last_timestamp)[0]
+
+    if last_prices_timestamp < block.timestamp:  # <------------ Update moving
+        #                                                   average if needed.
+
+        last_prices: uint256 = self.last_prices
+        ma_time: uint256 = self._unpack_3(self.packed_rebalancing_params)[2]
+        alpha: uint256 = MATH.wad_exp(
+            -convert(
+                unsafe_sub(block.timestamp, last_prices_timestamp) * 10**18 / ma_time,
+                int256,
+            )
+        )
+
+        # ---- We cap state price that goes into the EMA with 2 x price_scale.
+        return (
+            min(last_prices, 2 * price_scale) * (10**18 - alpha) +
+            price_oracle * alpha
+        ) / 10**18
+
+    return price_oracle
+
+
 @external
 @view
 def fee_receiver() -> address:
@@ -1459,6 +1563,16 @@ def fee_receiver() -> address:
     return factory.fee_receiver()
 
 
+@external
+@view
+def admin() -> address:
+    """
+    @notice Returns the address of the pool's admin.
+    @return address Admin.
+    """
+    return factory.admin()
+
+
 @external
 @view
 def calc_token_amount(amounts: uint256[N_COINS], deposit: bool) -> uint256:
@@ -1516,7 +1630,7 @@ def lp_price() -> uint256:
             0th index
     @return uint256 LP price.
     """
-    return 2 * self.virtual_price * isqrt(self.cached_price_oracle) / 10**18                # TODO: Check precision.
+    return 2 * self.virtual_price * isqrt(self.internal_price_oracle() * 10**18) / 10**18
 
 
 @external
@@ -1542,32 +1656,9 @@ def price_oracle() -> uint256:
     @dev The oracle is an exponential moving average, with a periodicity
          determined by `self.ma_time`. The aggregated prices are cached state
          prices (dy/dx) calculated AFTER the latest trade.
-    @param k The index of the coin.
     @return uint256 Price oracle value of kth coin.
     """
-    price_oracle: uint256 = self.cached_price_oracle
-    price_scale: uint256 = self.cached_price_scale
-    last_prices_timestamp: uint256 = self.last_prices_timestamp
-
-    if last_prices_timestamp < block.timestamp:  # <------------ Update moving
-        #                                                   average if needed.
-
-        last_prices: uint256 = self.last_prices
-        ma_time: uint256 = self._unpack(self.packed_rebalancing_params)[2]
-        alpha: uint256 = MATH.wad_exp(
-            -convert(
-                (block.timestamp - last_prices_timestamp) * 10**18 / ma_time,
-                int256,
-            )
-        )
-
-        # ---- We cap state price that goes into the EMA with 2 x price_scale.
-        return (
-            min(last_prices, 2 * price_scale) * (10**18 - alpha) +
-            price_oracle * alpha
-        ) / 10**18
-
-    return price_oracle
+    return self.internal_price_oracle()
 
 
 @external
@@ -1585,14 +1676,17 @@ def xcp_oracle() -> uint256:
     @return uint256 Oracle value of xcp.
     """
 
-    last_prices_timestamp: uint256 = self.last_prices_timestamp
+    last_prices_timestamp: uint256 = self._unpack_2(self.last_timestamp)[1]
     cached_xcp_oracle: uint256 = self.cached_xcp_oracle
 
     if last_prices_timestamp < block.timestamp:
 
         alpha: uint256 = MATH.wad_exp(
             -convert(
-                (block.timestamp - last_prices_timestamp) * 10**18 / self.xcp_ma_time,
+                unsafe_div(
+                    unsafe_sub(block.timestamp, last_prices_timestamp) * 10**18,
+                    self.xcp_ma_time
+                ),
                 int256,
             )
         )
@@ -1611,7 +1705,6 @@ def price_scale() -> uint256:
             at index 0.
     @dev Price scale determines the price band around which liquidity is
          concentrated.
-    @param k The index of the coin.
     @return uint256 Price scale of coin.
     """
     return self.cached_price_scale
@@ -1689,7 +1782,7 @@ def mid_fee() -> uint256:
     @notice Returns the current mid fee
     @return uint256 mid_fee value.
     """
-    return self._unpack(self.packed_fee_params)[0]
+    return self._unpack_3(self.packed_fee_params)[0]
 
 
 @view
@@ -1699,7 +1792,7 @@ def out_fee() -> uint256:
     @notice Returns the current out fee
     @return uint256 out_fee value.
     """
-    return self._unpack(self.packed_fee_params)[1]
+    return self._unpack_3(self.packed_fee_params)[1]
 
 
 @view
@@ -1709,7 +1802,7 @@ def fee_gamma() -> uint256:
     @notice Returns the current fee gamma
     @return uint256 fee_gamma value.
     """
-    return self._unpack(self.packed_fee_params)[2]
+    return self._unpack_3(self.packed_fee_params)[2]
 
 
 @view
@@ -1719,7 +1812,7 @@ def allowed_extra_profit() -> uint256:
     @notice Returns the current allowed extra profit
     @return uint256 allowed_extra_profit value.
     """
-    return self._unpack(self.packed_rebalancing_params)[0]
+    return self._unpack_3(self.packed_rebalancing_params)[0]
 
 
 @view
@@ -1729,7 +1822,7 @@ def adjustment_step() -> uint256:
     @notice Returns the current adjustment step
     @return uint256 adjustment_step value.
     """
-    return self._unpack(self.packed_rebalancing_params)[1]
+    return self._unpack_3(self.packed_rebalancing_params)[1]
 
 
 @view
@@ -1741,7 +1834,7 @@ def ma_time() -> uint256:
          One can expect off-by-one errors here.
     @return uint256 ma_time value.
     """
-    return self._unpack(self.packed_rebalancing_params)[2] * 694 / 1000
+    return self._unpack_3(self.packed_rebalancing_params)[2] * 694 / 1000
 
 
 @view
@@ -1879,7 +1972,7 @@ def apply_new_parameters(
     new_out_fee: uint256 = _new_out_fee
     new_fee_gamma: uint256 = _new_fee_gamma
 
-    current_fee_params: uint256[3] = self._unpack(self.packed_fee_params)
+    current_fee_params: uint256[3] = self._unpack_3(self.packed_fee_params)
 
     if new_out_fee < MAX_FEE + 1:
         assert new_out_fee > MIN_FEE - 1  # dev: fee is out of range
@@ -1895,7 +1988,7 @@ def apply_new_parameters(
     else:
         new_fee_gamma = current_fee_params[2]
 
-    self.packed_fee_params = self._pack([new_mid_fee, new_out_fee, new_fee_gamma])
+    self.packed_fee_params = self._pack_3([new_mid_fee, new_out_fee, new_fee_gamma])
 
     # ----------------- Set liquidity rebalancing parameters -----------------
 
@@ -1903,7 +1996,7 @@ def apply_new_parameters(
     new_adjustment_step: uint256 = _new_adjustment_step
     new_ma_time: uint256 = _new_ma_time
 
-    current_rebalancing_params: uint256[3] = self._unpack(self.packed_rebalancing_params)
+    current_rebalancing_params: uint256[3] = self._unpack_3(self.packed_rebalancing_params)
 
     if new_allowed_extra_profit > 10**18:
         new_allowed_extra_profit = current_rebalancing_params[0]
@@ -1916,7 +2009,7 @@ def apply_new_parameters(
     else:
         new_ma_time = current_rebalancing_params[2]
 
-    self.packed_rebalancing_params = self._pack(
+    self.packed_rebalancing_params = self._pack_3(
         [new_allowed_extra_profit, new_adjustment_step, new_ma_time]
     )
 
diff --git a/contracts/main/LiquidityGauge.vy b/contracts/main/LiquidityGauge.vy
index 6174a6e7..c30fd36f 100644
--- a/contracts/main/LiquidityGauge.vy
+++ b/contracts/main/LiquidityGauge.vy
@@ -1,5 +1,5 @@
-# @version 0.3.10
-#pragma optimize gas
+# pragma version 0.3.10
+# pragma optimize gas
 """
 @title LiquidityGaugeV6
 @author Curve.Fi
@@ -89,9 +89,9 @@ MAX_REWARDS: constant(uint256) = 8
 TOKENLESS_PRODUCTION: constant(uint256) = 40
 WEEK: constant(uint256) = 604800
 
-VERSION: constant(String[8]) = "v6.0.0"  # <- updated from v5.0.0 (adds `create_from_blueprint` pattern)
+VERSION: constant(String[8]) = "v6.1.0"  # <- updated from v6.0.0 (makes rewards semi-permissionless)
 
-EIP712_TYPEHASH: constant(bytes32) = keccak256("EIP712Domain(string name,string version,uint256 chainId,address verifyingContract)")
+EIP712_TYPEHASH: constant(bytes32) = keccak256("EIP712Domain(string name,string version,uint256 chainId,address verifyingContract,bytes32 salt)")
 EIP2612_TYPEHASH: constant(bytes32) = keccak256("Permit(address owner,address spender,uint256 value,uint256 nonce,uint256 deadline)")
 
 VERSION_HASH: constant(bytes32) = keccak256(VERSION)
@@ -112,16 +112,16 @@ balanceOf: public(HashMap[address, uint256])
 totalSupply: public(uint256)
 allowance: public(HashMap[address, HashMap[address, uint256]])
 
-name: public(String[64])
-symbol: public(String[40])
+name: public(immutable(String[64]))
+symbol: public(immutable(String[40]))
 
 # ERC2612
 nonces: public(HashMap[address, uint256])
 
 # Gauge
-factory: public(address)
+factory: public(immutable(address))
+lp_token: public(immutable(address))
 manager: public(address)
-lp_token: public(address)
 
 is_killed: public(bool)
 
@@ -170,15 +170,15 @@ def __init__(_lp_token: address):
     @notice Contract constructor
     @param _lp_token Liquidity Pool contract address
     """
-    self.lp_token = _lp_token
-    self.factory = msg.sender
+    lp_token = _lp_token
+    factory = msg.sender
     self.manager = msg.sender
 
-    symbol: String[32] = ERC20Extended(_lp_token).symbol()
-    name: String[64] = concat("Curve.fi ", symbol, " Gauge Deposit")
+    _symbol: String[32] = ERC20Extended(_lp_token).symbol()
+    _name: String[64] = concat("Curve.fi ", symbol, " Gauge Deposit")
 
-    self.name = name
-    self.symbol = concat(symbol, "-gauge")
+    name = _name
+    symbol = concat(_symbol, "-gauge")
 
     self.period_timestamp[0] = block.timestamp
     self.inflation_params = (
@@ -186,7 +186,7 @@ def __init__(_lp_token: address):
         + CRV20(CRV).rate()
     )
 
-    NAME_HASH = keccak256(name)
+    NAME_HASH = keccak256(_name)
     salt = block.prevhash
     CACHED_CHAIN_ID = chain.id
     CACHED_DOMAIN_SEPARATOR = keccak256(
@@ -426,7 +426,7 @@ def deposit(_value: uint256, _addr: address = msg.sender, _claim_rewards: bool =
 
         self._update_liquidity_limit(_addr, new_balance, total_supply)
 
-        ERC20(self.lp_token).transferFrom(msg.sender, self, _value)
+        ERC20(lp_token).transferFrom(msg.sender, self, _value)
 
         log Deposit(_addr, _value)
         log Transfer(empty(address), _addr, _value)
@@ -455,7 +455,7 @@ def withdraw(_value: uint256, _claim_rewards: bool = False):
 
         self._update_liquidity_limit(msg.sender, new_balance, total_supply)
 
-        ERC20(self.lp_token).transfer(msg.sender, _value)
+        ERC20(lp_token).transfer(msg.sender, _value)
 
     log Withdraw(msg.sender, _value)
     log Transfer(msg.sender, empty(address), _value)
@@ -478,7 +478,7 @@ def claim_rewards(_addr: address = msg.sender, _receiver: address = empty(addres
 
 @external
 @nonreentrant('lock')
-def transferFrom(_from: address, _to :address, _value: uint256) -> bool:
+def transferFrom(_from: address, _to: address, _value: uint256) -> bool:
     """
      @notice Transfer tokens from one address to another.
      @dev Transferring claims pending reward tokens for the sender and receiver
@@ -488,7 +488,9 @@ def transferFrom(_from: address, _to :address, _value: uint256) -> bool:
     """
     _allowance: uint256 = self.allowance[_from][msg.sender]
     if _allowance != max_value(uint256):
-        self.allowance[_from][msg.sender] = _allowance - _value
+        _new_allowance: uint256 = _allowance - _value
+        self.allowance[_from][msg.sender] = _new_allowance
+        log Approval(_from, msg.sender, _new_allowance)
 
     self._transfer(_from, _to, _value)
 
@@ -572,7 +574,7 @@ def permit(
     assert ecrecover(digest, _v, _r, _s) == _owner  # dev: invalid signature
 
     self.allowance[_owner][_spender] = _value
-    self.nonces[_owner] = nonce + 1
+    self.nonces[_owner] = unsafe_add(nonce, 1)
 
     log Approval(_owner, _spender, _value)
     return True
@@ -669,7 +671,7 @@ def set_gauge_manager(_gauge_manager: address):
         method, but only for the gauge which they are the manager of.
     @param _gauge_manager The account to set as the new manager of the gauge.
     """
-    assert msg.sender in [self.manager, Factory(self.factory).admin()]  # dev: only manager or factory admin
+    assert msg.sender in [self.manager, Factory(factory).admin()]  # dev: only manager or factory admin
 
     self.manager = _gauge_manager
     log SetGaugeManager(_gauge_manager)
@@ -719,7 +721,7 @@ def add_reward(_reward_token: address, _distributor: address):
     @param _reward_token The token to add as an additional reward
     @param _distributor Address permitted to fund this contract with the reward token
     """
-    assert msg.sender in [self.manager, Factory(self.factory).admin()]  # dev: only manager or factory admin
+    assert msg.sender in [self.manager, Factory(factory).admin()]  # dev: only manager or factory admin
     assert _distributor != empty(address)  # dev: distributor cannot be zero address
 
     reward_count: uint256 = self.reward_count
@@ -740,7 +742,7 @@ def set_reward_distributor(_reward_token: address, _distributor: address):
     """
     current_distributor: address = self.reward_data[_reward_token].distributor
 
-    assert msg.sender in [current_distributor, Factory(self.factory).admin(), self.manager]
+    assert msg.sender in [current_distributor, Factory(factory).admin(), self.manager]
     assert current_distributor != empty(address)
     assert _distributor != empty(address)
 
@@ -754,7 +756,7 @@ def set_killed(_is_killed: bool):
     @dev When killed, the gauge always yields a rate of 0 and so cannot mint CRV
     @param _is_killed Killed status to set
     """
-    assert msg.sender == Factory(self.factory).admin()  # dev: only owner
+    assert msg.sender == Factory(factory).admin()  # dev: only owner
 
     self.is_killed = _is_killed
 
diff --git a/contracts/mocks/ERC20Mock.vy b/contracts/mocks/ERC20Mock.vy
index a7c4b5b1..c3e50b7a 100644
--- a/contracts/mocks/ERC20Mock.vy
+++ b/contracts/mocks/ERC20Mock.vy
@@ -60,3 +60,12 @@ def approve(_spender: address, _value: uint256) -> bool:
     self.allowances[msg.sender][_spender] = _value
     log Approval(msg.sender, _spender, _value)
     return True
+
+
+@external
+def _mint_for_testing(_target: address, _value: uint256) -> bool:
+    self.totalSupply += _value
+    self.balanceOf[_target] += _value
+    log Transfer(ZERO_ADDRESS, _target, _value)
+
+    return True
diff --git a/contracts/old/CurveCryptoSwap2Math.vy b/contracts/old/CurveCryptoSwap2Math.vy
new file mode 100644
index 00000000..724c8d40
--- /dev/null
+++ b/contracts/old/CurveCryptoSwap2Math.vy
@@ -0,0 +1,193 @@
+# pragma version 0.3.10
+# pragma optimize gas
+
+N_COINS: constant(uint256) = 2
+A_MULTIPLIER: constant(uint256) = 10000
+
+MIN_GAMMA: constant(uint256) = 10**10
+MAX_GAMMA: constant(uint256) = 2 * 10**16
+
+MIN_A: constant(uint256) = N_COINS**N_COINS * A_MULTIPLIER / 10
+MAX_A: constant(uint256) = N_COINS**N_COINS * A_MULTIPLIER * 1000
+
+
+@internal
+@pure
+def geometric_mean(unsorted_x: uint256[N_COINS], sort: bool) -> uint256:
+    """
+    (x[0] * x[1] * ...) ** (1/N)
+    """
+    x: uint256[N_COINS] = unsorted_x
+    if sort and x[0] < x[1]:
+        x = [unsorted_x[1], unsorted_x[0]]
+    D: uint256 = x[0]
+    diff: uint256 = 0
+    for i in range(255):
+        D_prev: uint256 = D
+        # tmp: uint256 = 10**18
+        # for _x in x:
+        #     tmp = tmp * _x / D
+        # D = D * ((N_COINS - 1) * 10**18 + tmp) / (N_COINS * 10**18)
+        # line below makes it for 2 coins
+        D = (D + x[0] * x[1] / D) / N_COINS
+        if D > D_prev:
+            diff = D - D_prev
+        else:
+            diff = D_prev - D
+        if diff <= 1 or diff * 10**18 < D:
+            return D
+    raise "Did not converge"
+
+
+@external
+@pure
+def newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i: uint256) -> uint256:
+    """
+    Calculating x[i] given other balances x[0..N_COINS-1] and invariant D
+    ANN = A * N**N
+    """
+    # Safety checks
+    assert ANN > MIN_A - 1 and ANN < MAX_A + 1  # dev: unsafe values A
+    assert gamma > MIN_GAMMA - 1 and gamma < MAX_GAMMA + 1  # dev: unsafe values gamma
+    assert D > 10**17 - 1 and D < 10**15 * 10**18 + 1 # dev: unsafe values D
+
+    x_j: uint256 = x[1 - i]
+    y: uint256 = D**2 / (x_j * N_COINS**2)
+    K0_i: uint256 = (10**18 * N_COINS) * x_j / D
+    # S_i = x_j
+
+    # frac = x_j * 1e18 / D => frac = K0_i / N_COINS
+    assert (K0_i > 10**16*N_COINS - 1) and (K0_i < 10**20*N_COINS + 1)  # dev: unsafe values x[i]
+
+    # x_sorted: uint256[N_COINS] = x
+    # x_sorted[i] = 0
+    # x_sorted = self.sort(x_sorted)  # From high to low
+    # x[not i] instead of x_sorted since x_soted has only 1 element
+
+    convergence_limit: uint256 = max(max(x_j / 10**14, D / 10**14), 100)
+
+    for j in range(255):
+        y_prev: uint256 = y
+
+        K0: uint256 = K0_i * y * N_COINS / D
+        S: uint256 = x_j + y
+
+        _g1k0: uint256 = gamma + 10**18
+        if _g1k0 > K0:
+            _g1k0 = _g1k0 - K0 + 1
+        else:
+            _g1k0 = K0 - _g1k0 + 1
+
+        # D / (A * N**N) * _g1k0**2 / gamma**2
+        mul1: uint256 = 10**18 * D / gamma * _g1k0 / gamma * _g1k0 * A_MULTIPLIER / ANN
+
+        # 2*K0 / _g1k0
+        mul2: uint256 = 10**18 + (2 * 10**18) * K0 / _g1k0
+
+        yfprime: uint256 = 10**18 * y + S * mul2 + mul1
+        _dyfprime: uint256 = D * mul2
+        if yfprime < _dyfprime:
+            y = y_prev / 2
+            continue
+        else:
+            yfprime -= _dyfprime
+        fprime: uint256 = yfprime / y
+
+        # y -= f / f_prime;  y = (y * fprime - f) / fprime
+        # y = (yfprime + 10**18 * D - 10**18 * S) // fprime + mul1 // fprime * (10**18 - K0) // K0
+        y_minus: uint256 = mul1 / fprime
+        y_plus: uint256 = (yfprime + 10**18 * D) / fprime + y_minus * 10**18 / K0
+        y_minus += 10**18 * S / fprime
+
+        if y_plus < y_minus:
+            y = y_prev / 2
+        else:
+            y = y_plus - y_minus
+
+        diff: uint256 = 0
+        if y > y_prev:
+            diff = y - y_prev
+        else:
+            diff = y_prev - y
+        if diff < max(convergence_limit, y / 10**14):
+            frac: uint256 = y * 10**18 / D
+            assert (frac > 10**16 - 1) and (frac < 10**20 + 1)  # dev: unsafe value for y
+            return y
+
+    raise "Did not converge"
+
+
+@external
+@view
+def newton_D(ANN: uint256, gamma: uint256, x_unsorted: uint256[N_COINS]) -> uint256:
+    """
+    Finding the invariant using Newton method.
+    ANN is higher by the factor A_MULTIPLIER
+    ANN is already A * N**N
+
+    Currently uses 60k gas
+    """
+    # Safety checks
+    assert ANN > MIN_A - 1 and ANN < MAX_A + 1  # dev: unsafe values A
+    assert gamma > MIN_GAMMA - 1 and gamma < MAX_GAMMA + 1  # dev: unsafe values gamma
+
+    # Initial value of invariant D is that for constant-product invariant
+    x: uint256[N_COINS] = x_unsorted
+    if x[0] < x[1]:
+        x = [x_unsorted[1], x_unsorted[0]]
+
+    assert x[0] > 10**9 - 1 and x[0] < 10**15 * 10**18 + 1  # dev: unsafe values x[0]
+    assert x[1] * 10**18 / x[0] > 10**14 - 1  # dev: unsafe values x[i] (input)
+
+    D: uint256 = N_COINS * self.geometric_mean(x, False)
+    S: uint256 = x[0] + x[1]
+
+    for i in range(255):
+        D_prev: uint256 = D
+
+        # K0: uint256 = 10**18
+        # for _x in x:
+        #     K0 = K0 * _x * N_COINS / D
+        # collapsed for 2 coins
+        K0: uint256 = (10**18 * N_COINS**2) * x[0] / D * x[1] / D
+
+        _g1k0: uint256 = gamma + 10**18
+        if _g1k0 > K0:
+            _g1k0 = _g1k0 - K0 + 1
+        else:
+            _g1k0 = K0 - _g1k0 + 1
+
+        # D / (A * N**N) * _g1k0**2 / gamma**2
+        mul1: uint256 = 10**18 * D / gamma * _g1k0 / gamma * _g1k0 * A_MULTIPLIER / ANN
+
+        # 2*N*K0 / _g1k0
+        mul2: uint256 = (2 * 10**18) * N_COINS * K0 / _g1k0
+
+        neg_fprime: uint256 = (S + S * mul2 / 10**18) + mul1 * N_COINS / K0 - mul2 * D / 10**18
+
+        # D -= f / fprime
+        D_plus: uint256 = D * (neg_fprime + S) / neg_fprime
+        D_minus: uint256 = D*D / neg_fprime
+        if 10**18 > K0:
+            D_minus += D * (mul1 / neg_fprime) / 10**18 * (10**18 - K0) / K0
+        else:
+            D_minus -= D * (mul1 / neg_fprime) / 10**18 * (K0 - 10**18) / K0
+
+        if D_plus > D_minus:
+            D = D_plus - D_minus
+        else:
+            D = (D_minus - D_plus) / 2
+
+        diff: uint256 = 0
+        if D > D_prev:
+            diff = D - D_prev
+        else:
+            diff = D_prev - D
+        if diff * 10**14 < max(10**16, D):  # Could reduce precision for gas efficiency here
+            # Test that we are safe with the next newton_y
+            for _x in x:
+                frac: uint256 = _x * 10**18 / D
+                assert (frac >= 10**16 - 1) and (frac < 10**20 + 1)  # dev: unsafe values x[i]
+            return D
+
+    raise "Did not converge"
diff --git a/hardhat.config.js b/hardhat.config.js
deleted file mode 100644
index d2d23fda..00000000
--- a/hardhat.config.js
+++ /dev/null
@@ -1,16 +0,0 @@
-
-// See https://hardhat.org/config/ for config options.
-module.exports = {
-  networks: {
-    hardhat: {
-      hardfork: "shanghai",
-      // Base fee of 0 allows use of 0 gas price when testing
-      initialBaseFeePerGas: 0,
-      accounts: {
-        mnemonic: "test test test test test test test test test test test junk",
-        path: "m/44'/60'/0'",
-        count: 5
-      }
-    },
-  },
-};
diff --git a/requirements.txt b/requirements.txt
index b64f89be..c950036a 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -19,5 +19,5 @@ pdbpp
 hypothesis>=6.68.1
 
 # vyper and dev framework:
-git+https://github.com/vyperlang/titanoboa@40f5bfcc2afe212bb5a6f5026148f3625596ded7
+git+https://github.com/vyperlang/titanoboa@878da481a5adc3963b07b22355c80a85a9cfecb7
 vyper>=0.3.10
diff --git a/requirements_ape.txt b/requirements_ape.txt
deleted file mode 100644
index 90a7bd0b..00000000
--- a/requirements_ape.txt
+++ /dev/null
@@ -1,20 +0,0 @@
-# testing:
-ipython
-pytest
-pdbpp
-
-# vyper and dev framework:
-git+https://github.com/vyperlang/titanoboa@c1d741c26b34798ec1620859c7f3d8f42416e4be
-vyper>=0.3.9
-eth-ape==0.6.19
-ape-etherscan==0.6.10
-ape-hardhat==0.6.12
-ape-alchemy==0.6.4
-ape-vyper==0.6.10
-python-dotenv==1.0.0
-
-# prices api:
-pycoingecko
-
-# linting:
-pre-commit
diff --git a/tests/fixtures/accounts.py b/tests/fixtures/accounts.py
index f947969c..70b91065 100644
--- a/tests/fixtures/accounts.py
+++ b/tests/fixtures/accounts.py
@@ -16,8 +16,13 @@ def owner():
 
 
 @pytest.fixture(scope="module")
-def factory_admin(tricrypto_factory):
-    return tricrypto_factory.admin()
+def hacker():
+    return boa.env.generate_address()
+
+
+@pytest.fixture(scope="module")
+def factory_admin(factory):
+    return factory.admin()
 
 
 @pytest.fixture(scope="module")
@@ -32,6 +37,13 @@ def user():
     return acc
 
 
+@pytest.fixture(scope="module")
+def user_b():
+    acc = boa.env.generate_address()
+    boa.env.set_balance(acc, 10**25)
+    return acc
+
+
 @pytest.fixture(scope="module")
 def users():
     accs = [i() for i in [boa.env.generate_address] * 10]
diff --git a/tests/fixtures/factory.py b/tests/fixtures/factory.py
index 338f542e..18fe3dc5 100644
--- a/tests/fixtures/factory.py
+++ b/tests/fixtures/factory.py
@@ -37,7 +37,7 @@ def views_contract(deployer):
 
 
 @pytest.fixture(scope="module")
-def tricrypto_factory(
+def factory(
     deployer,
     fee_receiver,
     owner,
diff --git a/tests/fixtures/pool.py b/tests/fixtures/pool.py
index 82bbaaa0..de4c2f66 100644
--- a/tests/fixtures/pool.py
+++ b/tests/fixtures/pool.py
@@ -18,7 +18,11 @@ def _get_deposit_amounts(amount_per_token_usd, initial_prices, coins):
 
 
 def _crypto_swap_with_deposit(
-    coins, user, tricrypto_swap, initial_prices, dollar_amt_each_coin=10**6
+    coins,
+    user,
+    twocrypto_swap,
+    initial_prices,
+    dollar_amt_each_coin=int(1.5 * 10**6),
 ):
 
     # add 1M of each token to the pool
@@ -34,13 +38,13 @@ def _crypto_swap_with_deposit(
 
         # approve crypto_swap to trade coin for user:
         with boa.env.prank(user):
-            coin.approve(tricrypto_swap, 2**256 - 1)
+            coin.approve(twocrypto_swap, 2**256 - 1)
 
     # Very first deposit
     with boa.env.prank(user):
-        tricrypto_swap.add_liquidity(quantities, 0)
+        twocrypto_swap.add_liquidity(quantities, 0)
 
-    return tricrypto_swap
+    return twocrypto_swap
 
 
 @pytest.fixture(scope="module")
@@ -62,7 +66,7 @@ def params():
 
 @pytest.fixture(scope="module")
 def swap(
-    tricrypto_factory,
+    factory,
     amm_interface,
     coins,
     params,
@@ -70,20 +74,20 @@ def swap(
 ):
 
     with boa.env.prank(deployer):
-        swap = tricrypto_factory.deploy_pool(
-            "Curve.fi USD<>WETH",
-            "USD<>WETH",
-            [coin.address for coin in coins],
-            0,  # <-------- 0th implementation index
-            params["A"],
-            params["gamma"],
-            params["mid_fee"],
-            params["out_fee"],
-            params["fee_gamma"],
-            params["allowed_extra_profit"],
-            params["adjustment_step"],
-            params["ma_time"],  # <--- no admin_fee needed
-            params["initial_prices"][1],
+        swap = factory.deploy_pool(
+            "Curve.fi USD<>WETH",  # _name: String[64]
+            "USD<>WETH",  # _symbol: String[32]
+            [coin.address for coin in coins],  # _coins: address[N_COINS]
+            0,  # implementation_id: uint256
+            params["A"],  # A: uint256
+            params["gamma"],  # gamma: uint256
+            params["mid_fee"],  # mid_fee: uint256
+            params["out_fee"],  # out_fee: uint256
+            params["fee_gamma"],  # fee_gamma: uint256
+            params["allowed_extra_profit"],  # allowed_extra_profit: uint256
+            params["adjustment_step"],  # adjustment_step: uint256
+            params["ma_time"],  # ma_exp_time: uint256
+            params["initial_prices"][1],  # initial_price: uint256
         )
 
     return amm_interface.at(swap)
@@ -91,11 +95,10 @@ def swap(
 
 @pytest.fixture(scope="module")
 def swap_multiprecision(
-    tricrypto_factory,
+    factory,
     amm_interface,
     stgusdc,
     deployer,
-    weth,
 ):
 
     # STG/USDC pool params (on deployment)
@@ -108,42 +111,36 @@ def swap_multiprecision(
         "fee_gamma": 230000000000000,
         "adjustment_step": 146000000000000,
         "ma_time": 866,
-        "initial_prices": 500000000000000000,
+        "initial_prices": 1777655918836068423,
     }
 
-    with boa.env.prank(deployer):
-        swap = tricrypto_factory.deploy_pool(
-            "Curve.fi STG/USDC",
-            "STGUSDC",
-            [coin.address for coin in stgusdc],
-            weth,
-            0,
-            _params["A"],
-            _params["gamma"],
-            _params["mid_fee"],
-            _params["out_fee"],
-            _params["fee_gamma"],
-            _params["allowed_extra_profit"],
-            _params["adjustment_step"],
-            _params["ma_time"],
-            _params["initial_prices"],
-        )
+    swap = factory.deploy_pool(
+        "Curve.fi STG/USDC",
+        "STGUSDC",
+        [coin.address for coin in stgusdc],  # _coins: address[N_COINS]
+        0,  # implementation_id: uint256
+        _params["A"],  # A: uint256
+        _params["gamma"],  # gamma: uint256
+        _params["mid_fee"],  # mid_fee: uint256
+        _params["out_fee"],  # out_fee: uint256
+        _params["fee_gamma"],  # fee_gamma: uint256
+        _params["allowed_extra_profit"],  # allowed_extra_profit: uint256
+        _params["adjustment_step"],  # adjustment_step: uint256
+        _params["ma_time"],  # ma_exp_time: uint256
+        _params["initial_prices"],  # initial_price: uint256
+        sender=deployer,
+    )
 
     return amm_interface.at(swap)
 
 
 @pytest.fixture(scope="module")
 def swap_with_deposit(swap, coins, user):
-    yield _crypto_swap_with_deposit(coins, user, swap, INITIAL_PRICES)
-
-
-@pytest.fixture(scope="module")
-def hyper_swap_with_deposit(hyper_swap, coins, user):
-    yield _crypto_swap_with_deposit(coins, user, hyper_swap, INITIAL_PRICES)
+    return _crypto_swap_with_deposit(coins, user, swap, INITIAL_PRICES)
 
 
 @pytest.fixture(scope="module")
 def yuge_swap(swap, coins, user):
-    yield _crypto_swap_with_deposit(
+    return _crypto_swap_with_deposit(
         coins, user, swap, INITIAL_PRICES, dollar_amt_each_coin=10**10
     )
diff --git a/tests/fixtures/tokens.py b/tests/fixtures/tokens.py
index bc553440..d16e384d 100644
--- a/tests/fixtures/tokens.py
+++ b/tests/fixtures/tokens.py
@@ -50,7 +50,7 @@ def coins(usd, weth):
 
 
 @pytest.fixture(scope="module")
-def stgusdc(usdt, weth):
+def stgusdc(stg, usdc):
     yield [stg, usdc]
 
 
diff --git a/tests/unitary/math/conftest.py b/tests/unitary/math/conftest.py
new file mode 100644
index 00000000..64bcd4b4
--- /dev/null
+++ b/tests/unitary/math/conftest.py
@@ -0,0 +1,14 @@
+import boa
+import pytest
+
+
+@pytest.fixture(scope="module")
+def math_optimized(deployer):
+    with boa.env.prank(deployer):
+        return boa.load("contracts/main/CurveCryptoMathOptimized2.vy")
+
+
+@pytest.fixture(scope="module")
+def math_unoptimized(deployer):
+    with boa.env.prank(deployer):
+        return boa.load("contracts/old/CurveCryptoSwap2Math.vy")
diff --git a/tests/unitary/math/fuzz_multicoin_curve.py b/tests/unitary/math/fuzz_multicoin_curve.py
new file mode 100644
index 00000000..84cd3e2f
--- /dev/null
+++ b/tests/unitary/math/fuzz_multicoin_curve.py
@@ -0,0 +1,163 @@
+# flake8: noqa
+import unittest
+from itertools import permutations
+
+import hypothesis.strategies as st
+from hypothesis import given, settings
+
+from tests.utils.simulation_int_many import (
+    Curve,
+    geometric_mean,
+    reduction_coefficient,
+    solve_D,
+    solve_x,
+)
+
+MAX_EXAMPLES_MEAN = 20000
+MAX_EXAMPLES_RED = 20000
+MAX_EXAMPLES_D = 10000
+MAX_EXAMPLES_Y = 5000
+MAX_EXAMPLES_YD = 100000
+MAX_EXAMPLES_NOLOSS = 100000
+MIN_FEE = 5e-5
+
+MIN_XD = 10**16
+MAX_XD = 10**20
+
+N_COINS = 2
+A_MUL = 10000
+MIN_A = int(N_COINS**N_COINS * A_MUL / 10)
+MAX_A = int(N_COINS**N_COINS * A_MUL * 1000)
+
+MIN_GAMMA = 10**10
+MAX_GAMMA = 2 * 10**15
+
+
+# Test with 2 coins
+class TestCurve(unittest.TestCase):
+    @given(
+        x=st.integers(10**9, 10**15 * 10**18),
+        y=st.integers(10**9, 10**15 * 10**18),
+    )
+    @settings(max_examples=MAX_EXAMPLES_MEAN)
+    def test_geometric_mean(self, x, y):
+        val = geometric_mean([x, y])
+        assert val > 0
+        diff = abs((x * y) ** (1 / 2) - val)
+        assert diff / val <= max(1e-10, 1 / min([x, y]))
+
+    @given(
+        x=st.integers(10**9, 10**15 * 10**18),
+        y=st.integers(10**9, 10**15 * 10**18),
+        gamma=st.integers(10**10, 10**18),
+    )
+    @settings(max_examples=MAX_EXAMPLES_RED)
+    def test_reduction_coefficient(self, x, y, gamma):
+        coeff = reduction_coefficient([x, y], gamma)
+        assert coeff <= 10**18
+
+        K = 2**2 * x * y / (x + y) ** 2
+        if gamma > 0:
+            K = (gamma / 1e18) / ((gamma / 1e18) + 1 - K)
+        assert abs(coeff / 1e18 - K) <= 1e-7
+
+    @given(
+        A=st.integers(MIN_A, MAX_A),
+        x=st.integers(10**18, 10**15 * 10**18),  # 1 USD to 1e15 USD
+        yx=st.integers(
+            10**14, 10**18
+        ),  # <- ratio 1e18 * y/x, typically 1e18 * 1
+        perm=st.integers(0, 1),  # <- permutation mapping to values
+        gamma=st.integers(MIN_GAMMA, MAX_GAMMA),
+    )
+    @settings(max_examples=MAX_EXAMPLES_D)
+    def test_D_convergence(self, A, x, yx, perm, gamma):
+        # Price not needed for convergence testing
+        pmap = list(permutations(range(2)))
+
+        y = x * yx // 10**18
+        curve = Curve(A, gamma, 10**18, 2)
+        curve.x = [0] * 2
+        i, j = pmap[perm]
+        curve.x[i] = x
+        curve.x[j] = y
+        assert curve.D() > 0
+
+    @given(
+        A=st.integers(MIN_A, MAX_A),
+        x=st.integers(10**17, 10**15 * 10**18),  # $0.1 .. $1e15
+        yx=st.integers(10**15, 10**21),
+        gamma=st.integers(MIN_GAMMA, MAX_GAMMA),
+        i=st.integers(0, 1),
+        inx=st.integers(10**15, 10**21),
+    )
+    @settings(max_examples=MAX_EXAMPLES_Y)
+    def test_y_convergence(self, A, x, yx, gamma, i, inx):
+        j = 1 - i
+        in_amount = x * inx // 10**18
+        y = x * yx // 10**18
+        curve = Curve(A, gamma, 10**18, 2)
+        curve.x = [x, y]
+        out_amount = curve.y(in_amount, i, j)
+        assert out_amount > 0
+
+    @given(
+        A=st.integers(MIN_A, MAX_A),
+        x=st.integers(10**17, 10**15 * 10**18),  # 0.1 USD to 1e15 USD
+        yx=st.integers(5 * 10**14, 20 * 10**20),
+        gamma=st.integers(MIN_GAMMA, MAX_GAMMA),
+        i=st.integers(0, 1),
+        inx=st.integers(3 * 10**15, 3 * 10**20),
+    )
+    @settings(max_examples=MAX_EXAMPLES_NOLOSS)
+    def test_y_noloss(self, A, x, yx, gamma, i, inx):
+        j = 1 - i
+        y = x * yx // 10**18
+        curve = Curve(A, gamma, 10**18, 2)
+        curve.x = [x, y]
+        in_amount = x * inx // 10**18
+        try:
+            out_amount = curve.y(in_amount, i, j)
+            D1 = curve.D()
+        except ValueError:
+            return  # Convergence checked separately - we deliberately try unsafe numbers
+        is_safe = all(
+            f >= MIN_XD and f <= MAX_XD
+            for f in [xx * 10**18 // D1 for xx in curve.x]
+        )
+        curve.x[i] = in_amount
+        curve.x[j] = out_amount
+        try:
+            D2 = curve.D()
+        except ValueError:
+            return  # Convergence checked separately - we deliberately try unsafe numbers
+        is_safe &= all(
+            f >= MIN_XD and f <= MAX_XD
+            for f in [xx * 10**18 // D2 for xx in curve.x]
+        )
+        if is_safe:
+            assert (
+                2 * (D1 - D2) / (D1 + D2) < MIN_FEE
+            )  # Only loss is prevented - gain is ok
+
+    @given(
+        A=st.integers(MIN_A, MAX_A),
+        D=st.integers(10**18, 10**15 * 10**18),  # 1 USD to 1e15 USD
+        xD=st.integers(MIN_XD, MAX_XD),
+        yD=st.integers(MIN_XD, MAX_XD),
+        gamma=st.integers(MIN_GAMMA, MAX_GAMMA),
+        j=st.integers(0, 1),
+    )
+    @settings(max_examples=MAX_EXAMPLES_YD)
+    def test_y_from_D(self, A, D, xD, yD, gamma, j):
+        xp = [D * xD // 10**18, D * yD // 10**18]
+        y = solve_x(A, gamma, xp, D, j)
+        xp[j] = y
+        D2 = solve_D(A, gamma, xp)
+        assert (
+            2 * (D - D2) / (D2 + D) < MIN_FEE
+        )  # Only loss is prevented - gain is ok
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/unitary/math/misc.py b/tests/unitary/math/misc.py
new file mode 100644
index 00000000..2a36a01c
--- /dev/null
+++ b/tests/unitary/math/misc.py
@@ -0,0 +1,39 @@
+from decimal import Decimal
+
+
+def get_y_n2_dec(ANN, gamma, x, D, i):
+
+    if i == 0:
+        m = 1
+    elif i == 1:
+        m = 0
+
+    A = Decimal(ANN) / 10**4 / 4
+    gamma = Decimal(gamma) / 10**18
+    x = [Decimal(_x) / 10**18 for _x in x]
+    D = Decimal(D) / 10**18
+
+    a = Decimal(16) * x[m] ** 3 / D**3
+    b = 4 * A * gamma**2 * x[m] - (4 * (3 + 2 * gamma) * x[m] ** 2) / D
+    c = (
+        D * (3 + 4 * gamma + (1 - 4 * A) * gamma**2) * x[m]
+        + 4 * A * gamma**2 * x[m] ** 2
+    )
+    d = -(Decimal(1) / 4) * D**3 * (1 + gamma) ** 2
+
+    delta0 = b**2 - 3 * a * c
+    delta1 = 2 * b**3 - 9 * a * b * c + 27 * a**2 * d
+    sqrt_arg = delta1**2 - 4 * delta0**3
+
+    if sqrt_arg < 0:
+        return [0, {}]
+
+    sqrt = sqrt_arg ** (Decimal(1) / 2)
+    cbrt_arg = (delta1 + sqrt) / 2
+    if cbrt_arg > 0:
+        C1 = cbrt_arg ** (Decimal(1) / 3)
+    else:
+        C1 = -((-cbrt_arg) ** (Decimal(1) / 3))
+    root = -(b + C1 + delta0 / C1) / (3 * a)
+
+    return [root, (a, b, c, d)]
diff --git a/tests/unitary/math/test_cbrt.py b/tests/unitary/math/test_cbrt.py
new file mode 100644
index 00000000..a287945d
--- /dev/null
+++ b/tests/unitary/math/test_cbrt.py
@@ -0,0 +1,47 @@
+import pytest
+from boa.test import strategy
+from hypothesis import example, given, settings
+from vyper.utils import SizeLimits
+
+SETTINGS = {"max_examples": 10000, "deadline": None}
+MAX_VAL = SizeLimits.MAX_UINT256
+MAX_CBRT_PRECISE_VAL = MAX_VAL // 10**36
+
+
+def test_cbrt_expected_output(cbrt_1e18_base, math_optimized):
+    vals = [9 * 10**18, 8 * 10**18, 10**18, 1]
+    correct_cbrts = [2080083823051904114, 2 * 10**18, 10**18, 10**12]
+    for ix, val in enumerate(vals):
+        assert math_optimized.internal._cbrt(val) == correct_cbrts[ix]
+        assert cbrt_1e18_base(val) == correct_cbrts[ix]
+
+
+@given(
+    val=strategy("uint256", min_value=0, max_value=MAX_CBRT_PRECISE_VAL - 1)
+)
+@settings(**SETTINGS)
+@example(0)
+@example(1)
+def test_cbrt_exact(math_optimized, cbrt_1e18_base, val):
+
+    cbrt_python = cbrt_1e18_base(val)
+    cbrt_vyper = math_optimized.internal._cbrt(val)
+
+    try:
+        assert cbrt_python == cbrt_vyper
+    except AssertionError:
+        assert abs(cbrt_python - cbrt_vyper) == 1
+        pytest.warn(f"cbrt_python != cbrt_vyper for val = {val}")
+
+
+@given(
+    val=strategy("uint256", min_value=MAX_CBRT_PRECISE_VAL, max_value=MAX_VAL)
+)
+@settings(**SETTINGS)
+@example(MAX_VAL)
+@example(MAX_CBRT_PRECISE_VAL)
+def test_cbrt_precision_loss_gte_limit(cbrt_1e18_base, math_optimized, val):
+
+    cbrt_vyper = math_optimized.internal._cbrt(val)
+    cbrt_python = cbrt_1e18_base(val)
+    assert cbrt_vyper == pytest.approx(cbrt_python)
diff --git a/tests/unitary/math/test_exp.py b/tests/unitary/math/test_exp.py
new file mode 100644
index 00000000..d5af9589
--- /dev/null
+++ b/tests/unitary/math/test_exp.py
@@ -0,0 +1,31 @@
+import math
+
+import boa
+import pytest
+from boa.test import strategy
+from hypothesis import given, settings
+from vyper.utils import SizeLimits
+
+
+@given(
+    strategy(
+        "int256",
+        min_value=SizeLimits.MIN_INT256,
+        max_value=SizeLimits.MAX_INT256,
+    )
+)
+@settings(max_examples=10000, deadline=None)
+def test_exp(math_optimized, x):
+
+    if x >= 135305999368893231589:
+        with boa.reverts("Math: wad_exp overflow"):
+            math_optimized.wad_exp(x)
+
+    elif x <= -42139678854452767551:
+        assert math_optimized.wad_exp(x) == 0
+
+    else:
+
+        exp_ideal = int(math.exp(x / 10**18) * 10**18)
+        exp_implementation = math_optimized.wad_exp(x)
+        assert exp_ideal == pytest.approx(exp_implementation)
diff --git a/tests/unitary/math/test_get_p.py b/tests/unitary/math/test_get_p.py
new file mode 100644
index 00000000..7d0eae5f
--- /dev/null
+++ b/tests/unitary/math/test_get_p.py
@@ -0,0 +1,113 @@
+import boa
+import pytest
+from boa.test import strategy
+from hypothesis import given, settings
+
+from tests.fixtures.pool import INITIAL_PRICES
+from tests.utils import approx
+from tests.utils.tokens import mint_for_testing
+
+SETTINGS = {"max_examples": 20, "deadline": None}
+
+
+# flake8: noqa: E501
+@pytest.fixture(scope="module")
+def dydx_safemath():
+
+    get_price_impl = """
+N_COINS: constant(uint256) = 2
+A_MULTIPLIER: constant(uint256) = 10000
+
+@external
+@view
+def get_p(
+    _xp: uint256[N_COINS], _D: uint256, _A_gamma: uint256[N_COINS]
+) -> uint256:
+
+    assert _D > 10**17 - 1 and _D < 10**15 * 10**18 + 1  # dev: unsafe D values
+
+    K0: uint256 = 4 * _xp[0] * _xp[1] / _D * 10**36 / _D
+    GK0: uint256  = (
+        2 * K0 * K0 / 10**36 * K0 / 10**36
+        + (_A_gamma[1] + 10**18)**2
+        - K0**2 / 10**36 * (2 * _A_gamma[1] + 3 * 10**18) / 10**18
+    )
+    NNAG2: uint256 = _A_gamma[0] * _A_gamma[1]**2 / A_MULTIPLIER
+    denominator: uint256 = GK0 + NNAG2 * _xp[0] / _D * K0 / 10**36
+    numerator: uint256 = _xp[0] * ( GK0 + NNAG2 * _xp[1] / _D * K0 / 10**36 ) / _xp[1]
+    return numerator * 10**18 / denominator
+"""
+    return boa.loads(get_price_impl)
+
+
+def _get_dydx_vyper(swap, price_calc):
+
+    xp = swap.internal.xp(
+        swap._storage.balances.get(),
+        swap.price_scale(),
+    )
+
+    return price_calc.get_p(xp, swap.D(), swap.internal._A_gamma())
+
+
+def _get_prices_vyper(swap, price_calc):
+
+    price_token_1_wrt_0 = _get_dydx_vyper(swap, price_calc)
+    return price_token_1_wrt_0 * swap.price_scale() // 10**18
+
+
+def _get_prices_numeric_nofee(swap, views, i):
+
+    if i == 0:  # token at index 1 is being pupmed.
+
+        dx = int(0.01 * 10**36 // INITIAL_PRICES[1])
+        dolla_out = views.internal._get_dy_nofee(1, 0, dx, swap)[0]
+        price = dolla_out * 10**18 // dx
+
+    else:  # token at index 1 is being dupmed.
+
+        dx = 10**16  # 0.01 USD
+        dy = views.internal._get_dy_nofee(0, 1, dx, swap)[0]
+        price = dx * 10**18 // dy
+
+    return price
+
+
+# ----- Tests -----
+
+
+@given(
+    dollar_amount=strategy("decimal", min_value=1e-5, max_value=5e8),
+)
+@settings(**SETTINGS)
+@pytest.mark.parametrize("i", [0, 1])
+def test_dxdy_similar(
+    yuge_swap,
+    dydx_safemath,
+    views_contract,
+    user,
+    dollar_amount,
+    coins,
+    i,
+):
+
+    previous_p = yuge_swap.price_scale()
+    j = 1 - i
+
+    amount_in = int(dollar_amount * 10**36 // INITIAL_PRICES[i])
+    mint_for_testing(coins[i], user, amount_in)
+    yuge_swap.exchange(i, j, amount_in, 0, sender=user)
+
+    dxdy_vyper = _get_prices_vyper(yuge_swap, dydx_safemath)
+    dxdy_swap = yuge_swap.last_prices()  # <-- we check unsafe impl here.
+    dxdy_numeric = _get_prices_numeric_nofee(yuge_swap, views_contract, i)
+
+    if i == 0:  # token at index 1 is being pupmed, so last_price should go up
+        assert dxdy_swap > previous_p
+        assert dxdy_numeric > previous_p
+    else:  # token at index 1 is being dupmed, so last_price should go down
+        assert dxdy_swap < previous_p
+        assert dxdy_numeric < previous_p
+
+    assert dxdy_vyper == dxdy_swap
+    assert approx(dxdy_vyper, dxdy_numeric, 1e-5)
diff --git a/tests/unitary/math/test_get_y.py b/tests/unitary/math/test_get_y.py
new file mode 100644
index 00000000..a082062b
--- /dev/null
+++ b/tests/unitary/math/test_get_y.py
@@ -0,0 +1,121 @@
+# flake8: noqa
+import time
+from decimal import Decimal
+
+import boa
+import pytest
+from hypothesis import given, note, settings
+from hypothesis import strategies as st
+
+N_COINS = 2
+# MAX_SAMPLES = 1000000  # Increase for fuzzing
+MAX_SAMPLES = 300
+N_CASES = 32
+
+A_MUL = 10000
+MIN_A = int(N_COINS**N_COINS * A_MUL / 10)
+MAX_A = int(N_COINS**N_COINS * A_MUL * 1000)
+
+MIN_GAMMA = 10**10
+MAX_GAMMA = 2 * 10**15
+
+pytest.current_case_id = 0
+pytest.negative_sqrt_arg = 0
+pytest.gas_original = 0
+pytest.gas_new = 0
+pytest.t_start = time.time()
+
+
+def inv_target_decimal_n2(A, gamma, x, D):
+    N = len(x)
+
+    x_prod = Decimal(1)
+    for x_i in x:
+        x_prod *= x_i
+    K0 = x_prod / (Decimal(D) / N) ** N
+    K0 *= 10**18
+
+    if gamma > 0:
+        # K = gamma**2 * K0 / (gamma + 10**18*(Decimal(1) - K0))**2
+        K = gamma**2 * K0 / (gamma + 10**18 - K0) ** 2 / 10**18
+    K *= A
+
+    f = (
+        K * D ** (N - 1) * sum(x)
+        + x_prod
+        - (K * D**N + (Decimal(D) / N) ** N)
+    )
+
+    return f
+
+
+def test_get_y_revert(math_contract):
+    a = 1723894848
+    gamma = 24009999997600
+    x = [112497148627520223862735198942112, 112327102289152450435452075003508]
+    D = 224824250915890636214130540882688
+    i = 0
+
+    with boa.reverts():
+        math_contract.newton_y(a, gamma, x, D, i)
+
+    with boa.reverts():
+        math_contract.get_y(a, gamma, x, D, i)
+
+
+@pytest.mark.parametrize(
+    "_tmp", range(N_CASES)
+)  # Create N_CASES independent test instances.
+@given(
+    A=st.integers(min_value=MIN_A, max_value=MAX_A),
+    D=st.integers(
+        min_value=10**18, max_value=10**14 * 10**18
+    ),  # 1 USD to 100T USD
+    xD=st.integers(
+        min_value=5 * 10**16, max_value=10**19
+    ),  # <- ratio 1e18 * x/D, typically 1e18 * 1
+    yD=st.integers(
+        min_value=5 * 10**16, max_value=10**19
+    ),  # <- ratio 1e18 * y/D, typically 1e18 * 1
+    gamma=st.integers(min_value=MIN_GAMMA, max_value=MAX_GAMMA),
+    j=st.integers(min_value=0, max_value=1),
+)
+@settings(max_examples=MAX_SAMPLES, deadline=None)
+def test_get_y(math_unoptimized, math_optimized, A, D, xD, yD, gamma, j, _tmp):
+    pytest.current_case_id += 1
+    X = [D * xD // 10**18, D * yD // 10**18]
+
+    A_dec = Decimal(A) / 10000 / 4
+
+    def calculate_F_by_y0(y0):
+        new_X = X[:]
+        new_X[j] = y0
+        return inv_target_decimal_n2(A_dec, gamma, new_X, D)
+
+    result_original = math_unoptimized.newton_y(A, gamma, X, D, j)
+    pytest.gas_original += math_unoptimized._computation.get_gas_used()
+
+    result_get_y, K0 = math_optimized.get_y(A, gamma, X, D, j)
+    pytest.gas_new += math_optimized._computation.get_gas_used()
+
+    note(
+        "{"
+        f"'ANN': {A}, 'GAMMA': {gamma}, 'x': {X}, 'D': {D}, 'index': {j}"
+        "}\n"
+    )
+
+    if K0 == 0:
+        pytest.negative_sqrt_arg += 1
+        return
+
+    if pytest.current_case_id % 1000 == 0:
+        print(
+            f"--- {pytest.current_case_id}\nPositive dy frac: {100*pytest.negative_sqrt_arg/pytest.current_case_id:.1f}%\t{time.time() - pytest.t_start:.1f} seconds.\n"
+            f"Gas advantage per call: {pytest.gas_original//pytest.current_case_id} {pytest.gas_new//pytest.current_case_id}\n"
+        )
+
+    assert abs(result_original - result_get_y) <= max(
+        10**4, result_original / 1e8
+    ) or abs(calculate_F_by_y0(result_get_y)) <= abs(
+        calculate_F_by_y0(result_original)
+    )
diff --git a/tests/unitary/math/test_log2.py b/tests/unitary/math/test_log2.py
new file mode 100644
index 00000000..4c94d9f6
--- /dev/null
+++ b/tests/unitary/math/test_log2.py
@@ -0,0 +1,28 @@
+import math
+
+from boa.test import strategy
+from hypothesis import given, settings
+from vyper.utils import SizeLimits
+
+
+@given(
+    strategy(
+        "uint256",
+        min_value=0,
+        max_value=SizeLimits.MAX_UINT256,
+    )
+)
+@settings(max_examples=10000, deadline=None)
+def test_log2(math_optimized, x):
+
+    if x == 0:
+        assert math_optimized.internal._snekmate_log_2(x, False) == 0
+        return
+
+    log2_ideal = int(math.log2(x))
+    log2_implementation = math_optimized.internal._snekmate_log_2(x, False)
+
+    try:
+        assert log2_ideal == log2_implementation
+    except:  # noqa: E722; there will be off-by-one cases
+        assert abs(log2_ideal - log2_implementation) == 1
diff --git a/tests/unitary/math/test_newton_D.py b/tests/unitary/math/test_newton_D.py
new file mode 100644
index 00000000..d1cfe5b7
--- /dev/null
+++ b/tests/unitary/math/test_newton_D.py
@@ -0,0 +1,221 @@
+# flake8: noqa
+import sys
+import time
+from decimal import Decimal
+
+import pytest
+from boa.vyper.contract import BoaError
+from hypothesis import given, settings
+from hypothesis import strategies as st
+
+import tests.utils.simulation_int_many as sim
+
+# Uncomment to be able to print when parallelized
+# sys.stdout = sys.stderr
+
+
+def inv_target_decimal_n2(A, gamma, x, D):
+    N = len(x)
+
+    x_prod = Decimal(1)
+    for x_i in x:
+        x_prod *= x_i
+    K0 = x_prod / (Decimal(D) / N) ** N
+    K0 *= 10**18
+
+    if gamma > 0:
+        # K = gamma**2 * K0 / (gamma + 10**18*(Decimal(1) - K0))**2
+        K = gamma**2 * K0 / (gamma + 10**18 - K0) ** 2 / 10**18
+    K *= A
+
+    f = (
+        K * D ** (N - 1) * sum(x)
+        + x_prod
+        - (K * D**N + (Decimal(D) / N) ** N)
+    )
+
+    return f
+
+
+N_COINS = 2
+# MAX_SAMPLES = 1000000  # Increase for fuzzing
+MAX_SAMPLES = 300  # Increase for fuzzing
+N_CASES = 1
+
+A_MUL = 10000
+MIN_A = int(N_COINS**N_COINS * A_MUL / 10)
+MAX_A = int(N_COINS**N_COINS * A_MUL * 1000)
+
+# gamma from 1e-8 up to 0.05
+MIN_GAMMA = 10**10
+MAX_GAMMA = 5 * 10**16
+
+MIN_XD = 10**17
+MAX_XD = 10**19
+
+pytest.progress = 0
+pytest.actually_tested = 0
+pytest.t_start = time.time()
+
+
+@pytest.mark.parametrize(
+    "_tmp", range(N_CASES)
+)  # Create N_CASES independent test instances.
+@given(
+    A=st.integers(min_value=MIN_A, max_value=MAX_A),
+    D=st.integers(
+        min_value=10**18, max_value=10**14 * 10**18
+    ),  # 1 USD to 100T USD
+    xD=st.integers(
+        min_value=MIN_XD, max_value=MAX_XD
+    ),  # <- ratio 1e18 * x/D, typically 1e18 * 1
+    yD=st.integers(
+        min_value=MIN_XD, max_value=MAX_XD
+    ),  # <- ratio 1e18 * y/D, typically 1e18 * 1
+    gamma=st.integers(min_value=MIN_GAMMA, max_value=MAX_GAMMA),
+    j=st.integers(min_value=0, max_value=1),
+    btcScalePrice=st.integers(min_value=10**2, max_value=10**7),
+    ethScalePrice=st.integers(min_value=10, max_value=10**5),
+    mid_fee=st.sampled_from(
+        [
+            int(0.7e-3 * 10**10),
+            int(1e-3 * 10**10),
+            int(1.2e-3 * 10**10),
+            int(4e-3 * 10**10),
+        ]
+    ),
+    out_fee=st.sampled_from([int(4.0e-3 * 10**10), int(10.0e-3 * 10**10)]),
+    fee_gamma=st.sampled_from([int(1e-2 * 1e18), int(2e-6 * 1e18)]),
+)
+@settings(max_examples=MAX_SAMPLES, deadline=None)
+def test_newton_D(
+    math_optimized,
+    math_unoptimized,
+    A,
+    D,
+    xD,
+    yD,
+    gamma,
+    j,
+    btcScalePrice,
+    ethScalePrice,
+    mid_fee,
+    out_fee,
+    fee_gamma,
+    _tmp,
+):
+    _test_newton_D(
+        math_optimized,
+        math_unoptimized,
+        A,
+        D,
+        xD,
+        yD,
+        gamma,
+        j,
+        btcScalePrice,
+        ethScalePrice,
+        mid_fee,
+        out_fee,
+        fee_gamma,
+        _tmp,
+    )
+
+
+def _test_newton_D(
+    math_optimized,
+    math_unoptimized,
+    A,
+    D,
+    xD,
+    yD,
+    gamma,
+    j,
+    btcScalePrice,
+    ethScalePrice,
+    mid_fee,
+    out_fee,
+    fee_gamma,
+    _tmp,
+):
+
+    is_safe = all(
+        f >= MIN_XD and f <= MAX_XD
+        for f in [xx * 10**18 // D for xx in [xD, yD]]
+    )
+
+    pytest.progress += 1
+    if pytest.progress % 1000 == 0 and pytest.actually_tested != 0:
+        print(
+            f"{pytest.progress} | {pytest.actually_tested} cases processed in {time.time()-pytest.t_start:.1f} seconds."
+        )
+    X = [D * xD // 10**18, D * yD // 10**18]
+
+    result_get_y = 0
+    get_y_failed = False
+    try:
+        (result_get_y, K0) = math_optimized.get_y(A, gamma, X, D, j)
+    except:
+        get_y_failed = True
+
+    if get_y_failed:
+        newton_y_failed = False
+        try:
+            math_optimized.newton_y(A, gamma, X, D, j)
+        except:
+            newton_y_failed = True
+
+    if get_y_failed and newton_y_failed:
+        return  # both canonical and new method fail, so we ignore.
+
+    if get_y_failed and not newton_y_failed and is_safe:
+        raise  # this is a problem
+
+    # dy should be positive
+    if result_get_y < X[j]:
+
+        price_scale = (btcScalePrice, ethScalePrice)
+        y = X[j]
+        dy = X[j] - result_get_y
+        dy -= 1
+
+        if j > 0:
+            dy = dy * 10**18 // price_scale[j - 1]
+
+        fee = sim.get_fee(X, fee_gamma, mid_fee, out_fee)
+        dy -= fee * dy // 10**10
+        y -= dy
+
+        if dy / X[j] <= 0.95:
+
+            pytest.actually_tested += 1
+            X[j] = y
+
+            case = (
+                "{"
+                f"'ANN': {A}, 'D': {D}, 'xD': {xD}, 'yD': {yD}, 'GAMMA': {gamma}, 'j': {j}, 'btcScalePrice': {btcScalePrice}, 'ethScalePrice': {ethScalePrice}, 'mid_fee': {mid_fee}, 'out_fee': {out_fee}, 'fee_gamma': {fee_gamma}"
+                "},\n"
+            )
+
+            result_sim = math_unoptimized.newton_D(A, gamma, X)
+            try:
+                result_contract = math_optimized.newton_D(A, gamma, X, K0)
+            except Exception as e:
+                # with open("log/newton_D_fail.txt", "a") as f:
+                #     f.write(case)
+                # with open("log/newton_D_fail_trace.txt", "a") as f:
+                #     f.write(str(e))
+                return
+
+            A_dec = Decimal(A) / 10000 / 4
+
+            def calculate_D_polynome(d):
+                d = Decimal(d)
+                return abs(inv_target_decimal_n2(A_dec, gamma, X, d))
+
+            assert abs(result_sim - result_contract) <= max(
+                10000, result_sim / 1e12
+            )
+
+            # with open("log/newton_D_pass.txt", "a") as f:
+            #     f.write(case)
diff --git a/tests/unitary/math/test_newton_D_ref.py b/tests/unitary/math/test_newton_D_ref.py
new file mode 100644
index 00000000..6852c39d
--- /dev/null
+++ b/tests/unitary/math/test_newton_D_ref.py
@@ -0,0 +1,192 @@
+# flake8: noqa
+import sys
+from decimal import Decimal
+
+import pytest
+from boa.vyper.contract import BoaError
+from hypothesis import given, settings
+from hypothesis import strategies as st
+
+import tests.utils.simulation_int_many as sim
+
+# sys.stdout = sys.stderr
+
+
+def inv_target_decimal_n2(A, gamma, x, D):
+    N = len(x)
+
+    x_prod = Decimal(1)
+    for x_i in x:
+        x_prod *= x_i
+    K0 = x_prod / (Decimal(D) / N) ** N
+    K0 *= 10**18
+
+    if gamma > 0:
+        # K = gamma**2 * K0 / (gamma + 10**18*(Decimal(1) - K0))**2
+        K = gamma**2 * K0 / (gamma + 10**18 - K0) ** 2 / 10**18
+    K *= A
+
+    f = (
+        K * D ** (N - 1) * sum(x)
+        + x_prod
+        - (K * D**N + (Decimal(D) / N) ** N)
+    )
+
+    return f
+
+
+N_COINS = 2
+# MAX_SAMPLES = 300000  # Increase for fuzzing
+MAX_SAMPLES = 300  # Increase for fuzzing
+N_CASES = 1
+
+A_MUL = 10000
+MIN_A = int(N_COINS**N_COINS * A_MUL / 10)
+MAX_A = int(N_COINS**N_COINS * A_MUL * 1000)
+
+# gamma from 1e-8 up to 0.05
+MIN_GAMMA = 10**10
+MAX_GAMMA = 2 * 10**15
+
+MIN_XD = 10**17
+MAX_XD = 10**19
+
+pytest.cases = 0
+
+
+@pytest.mark.parametrize(
+    "_tmp", range(N_CASES)
+)  # Create N_CASES independent test instances.
+@given(
+    A=st.integers(min_value=MIN_A, max_value=MAX_A),
+    D=st.integers(
+        min_value=10**18, max_value=10**14 * 10**18
+    ),  # 1 USD to 100T USD
+    xD=st.integers(
+        min_value=MIN_XD, max_value=MAX_XD
+    ),  # <- ratio 1e18 * x/D, typically 1e18 * 1
+    yD=st.integers(
+        min_value=MIN_XD, max_value=MAX_XD
+    ),  # <- ratio 1e18 * y/D, typically 1e18 * 1
+    gamma=st.integers(min_value=MIN_GAMMA, max_value=MAX_GAMMA),
+    j=st.integers(min_value=0, max_value=1),
+    asset_x_scale_price=st.integers(min_value=10**2, max_value=10**7),
+    asset_y_scale_price=st.integers(min_value=10, max_value=10**5),
+    mid_fee=st.sampled_from(
+        [
+            int(0.7e-3 * 10**10),
+            int(1e-3 * 10**10),
+            int(1.2e-3 * 10**10),
+            int(4e-3 * 10**10),
+        ]
+    ),
+    out_fee=st.sampled_from([int(4.0e-3 * 10**10), int(10.0e-3 * 10**10)]),
+    fee_gamma=st.sampled_from([int(1e-2 * 1e18), int(2e-6 * 1e18)]),
+)
+@settings(max_examples=MAX_SAMPLES, deadline=None)
+def test_newton_D(
+    math_optimized,
+    math_unoptimized,
+    A,
+    D,
+    xD,
+    yD,
+    gamma,
+    j,
+    asset_x_scale_price,
+    asset_y_scale_price,
+    mid_fee,
+    out_fee,
+    fee_gamma,
+    _tmp,
+):
+    _test_newton_D(
+        math_optimized,
+        math_unoptimized,
+        A,
+        D,
+        xD,
+        yD,
+        gamma,
+        j,
+        asset_x_scale_price,
+        asset_y_scale_price,
+        mid_fee,
+        out_fee,
+        fee_gamma,
+        _tmp,
+    )
+
+
+def _test_newton_D(
+    math_optimized,
+    math_unoptimized,
+    A,
+    D,
+    xD,
+    yD,
+    gamma,
+    j,
+    asset_x_scale_price,
+    asset_y_scale_price,
+    mid_fee,
+    out_fee,
+    fee_gamma,
+    _tmp,
+):
+
+    pytest.cases += 1
+    X = [D * xD // 10**18, D * yD // 10**18]
+    is_safe = all(
+        f >= MIN_XD and f <= MAX_XD
+        for f in [xx * 10**18 // D for xx in [xD, yD]]
+    )
+    try:
+        newton_y_output = math_unoptimized.newton_y(A, gamma, X, D, j)
+    except BoaError as e:
+        if is_safe:
+            raise
+        else:
+            return
+
+    (result_get_y, K0) = math_optimized.get_y(A, gamma, X, D, j)
+
+    # dy should be positive
+    if result_get_y < X[j]:
+
+        price_scale = (asset_x_scale_price, asset_y_scale_price)
+        y = X[j]
+        dy = X[j] - result_get_y
+        dy -= 1
+
+        if j > 0:
+            dy = dy * 10**18 // price_scale[j - 1]
+
+        fee = sim.get_fee(X, fee_gamma, mid_fee, out_fee)
+        dy -= fee * dy // 10**10
+        y -= dy
+
+        if dy / X[j] <= 0.95:
+
+            X[j] = y
+
+            # if pytest.cases % 1000 == 0:
+            #     print(f'> {pytest.cases}')
+            try:
+                result_sim = math_unoptimized.newton_D(A, gamma, X)
+            except:
+                # breakpoint()
+                raise  # this is a problem
+
+            try:
+                result_contract = math_optimized.newton_D(A, gamma, X, K0)
+            # except BoaError:
+            except:
+                raise
+
+            try:
+                assert abs(result_sim - result_contract) <= max(
+                    10000, result_sim / 1e12
+                )
+            except AssertionError:
+                raise
diff --git a/tests/unitary/math/test_packing.py b/tests/unitary/math/test_packing.py
new file mode 100644
index 00000000..de444c48
--- /dev/null
+++ b/tests/unitary/math/test_packing.py
@@ -0,0 +1,23 @@
+from boa.test import strategy
+from hypothesis import given, settings
+
+
+@given(val=strategy("uint256[3]", max_value=10**18))
+@settings(max_examples=10000, deadline=None)
+def test_pack_unpack_three_integers(swap, factory, val):
+    for contract in [swap, factory]:
+        packed = contract.internal._pack_3(val)
+        unpacked = swap.internal._unpack_3(packed)  # swap unpacks
+        for i in range(3):
+            assert unpacked[i] == val[i]
+
+
+@given(val=strategy("uint256[2]", max_value=2**128 - 1))
+@settings(max_examples=10000, deadline=None)
+def test_pack_unpack_2_integers(swap, val):
+
+    packed = swap.internal._pack_2(val[0], val[1])
+    unpacked = swap.internal._unpack_2(packed)  # swap unpacks
+
+    assert unpacked[0] == val[0]
+    assert unpacked[1] == val[1]
diff --git a/tests/unitary/pool/admin/test_commit_params.py b/tests/unitary/pool/admin/test_commit_params.py
new file mode 100644
index 00000000..215ec7e6
--- /dev/null
+++ b/tests/unitary/pool/admin/test_commit_params.py
@@ -0,0 +1,133 @@
+import copy
+
+import boa
+
+
+def _apply_new_params(swap, params):
+    swap.apply_new_parameters(
+        params["mid_fee"],
+        params["out_fee"],
+        params["fee_gamma"],
+        params["allowed_extra_profit"],
+        params["adjustment_step"],
+        params["ma_time"],
+        params["xcp_ma_time"],
+    )
+
+
+def test_commit_accept_mid_fee(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["mid_fee"] = p["mid_fee"] + 1
+    with boa.env.prank(factory_admin):
+        _apply_new_params(swap, p)
+
+    mid_fee = swap.internal._unpack_3(swap._storage.packed_fee_params.get())[0]
+    assert mid_fee == p["mid_fee"]
+
+
+def test_commit_accept_out_fee(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["out_fee"] = p["out_fee"] + 1
+    with boa.env.prank(factory_admin):
+        _apply_new_params(swap, p)
+
+    out_fee = swap.internal._unpack_3(swap._storage.packed_fee_params.get())[1]
+    assert out_fee == p["out_fee"]
+
+
+def test_commit_accept_fee_gamma(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["fee_gamma"] = 10**17
+    with boa.env.prank(factory_admin):
+        _apply_new_params(swap, p)
+
+    fee_gamma = swap.internal._unpack_3(swap._storage.packed_fee_params.get())[
+        2
+    ]
+    assert fee_gamma == p["fee_gamma"]
+
+
+def test_commit_accept_fee_params(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["mid_fee"] += 1
+    p["out_fee"] += 1
+    p["fee_gamma"] = 10**17
+
+    with boa.env.prank(swap.admin()):
+        _apply_new_params(swap, p)
+
+    fee_params = swap.internal._unpack_3(swap._storage.packed_fee_params.get())
+    assert fee_params[0] == p["mid_fee"]
+    assert fee_params[1] == p["out_fee"]
+    assert fee_params[2] == p["fee_gamma"]
+
+
+def test_commit_accept_allowed_extra_profit(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["allowed_extra_profit"] = 10**17
+    with boa.env.prank(factory_admin):
+        _apply_new_params(swap, p)
+
+    allowed_extra_profit = swap.internal._unpack_3(
+        swap._storage.packed_rebalancing_params.get()
+    )[0]
+    assert allowed_extra_profit == p["allowed_extra_profit"]
+
+
+def test_commit_accept_adjustment_step(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["adjustment_step"] = 10**17
+    with boa.env.prank(factory_admin):
+        _apply_new_params(swap, p)
+
+    adjustment_step = swap.internal._unpack_3(
+        swap._storage.packed_rebalancing_params.get()
+    )[1]
+    assert adjustment_step == p["adjustment_step"]
+
+
+def test_commit_accept_ma_time(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["ma_time"] = 872
+    with boa.env.prank(factory_admin):
+        _apply_new_params(swap, p)
+
+    ma_time = swap.internal._unpack_3(
+        swap._storage.packed_rebalancing_params.get()
+    )[2]
+    assert ma_time == p["ma_time"]
+
+
+def test_commit_accept_xcp_ma_time(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["xcp_ma_time"] = 872541
+    with boa.env.prank(factory_admin):
+        _apply_new_params(swap, p)
+
+    assert swap.xcp_ma_time() == p["xcp_ma_time"]
+
+
+def test_commit_accept_rebalancing_params(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["allowed_extra_profit"] = 10**17
+    p["adjustment_step"] = 10**17
+    p["ma_time"] = 1000
+
+    with boa.env.prank(factory_admin):
+        _apply_new_params(swap, p)
+
+    rebalancing_params = swap.internal._unpack_3(
+        swap._storage.packed_rebalancing_params.get()
+    )
+    assert rebalancing_params[0] == p["allowed_extra_profit"]
+    assert rebalancing_params[1] == p["adjustment_step"]
+    assert rebalancing_params[2] == p["ma_time"]
diff --git a/tests/unitary/pool/admin/test_ramp_A_gamma.py b/tests/unitary/pool/admin/test_ramp_A_gamma.py
new file mode 100644
index 00000000..f06126ca
--- /dev/null
+++ b/tests/unitary/pool/admin/test_ramp_A_gamma.py
@@ -0,0 +1,49 @@
+import copy
+
+import boa
+
+
+def test_ramp_A_gamma_up(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    future_A = p["A"] + 10000
+    future_gamma = p["gamma"] + 10000
+    future_time = boa.env.vm.state.timestamp + 86400
+
+    initial_A_gamma = [swap.A(), swap.gamma()]
+    swap.ramp_A_gamma(
+        future_A, future_gamma, future_time, sender=factory_admin
+    )
+
+    boa.env.time_travel(10000)
+    current_A_gamma = [swap.A(), swap.gamma()]
+    for i in range(2):
+        assert current_A_gamma[i] > initial_A_gamma[i]
+
+    boa.env.time_travel(76400)
+    current_A_gamma = [swap.A(), swap.gamma()]
+    assert current_A_gamma[0] == future_A
+    assert current_A_gamma[1] == future_gamma
+
+
+def test_ramp_A_gamma_down(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    future_A = p["A"] - 10000
+    future_gamma = p["gamma"] - 10000
+    future_time = boa.env.vm.state.timestamp + 86400
+
+    initial_A_gamma = [swap.A(), swap.gamma()]
+    swap.ramp_A_gamma(
+        future_A, future_gamma, future_time, sender=factory_admin
+    )
+
+    boa.env.time_travel(10000)
+    current_A_gamma = [swap.A(), swap.gamma()]
+    for i in range(2):
+        assert current_A_gamma[i] < initial_A_gamma[i]
+
+    boa.env.time_travel(76400)
+    current_A_gamma = [swap.A(), swap.gamma()]
+    assert current_A_gamma[0] == future_A
+    assert current_A_gamma[1] == future_gamma
diff --git a/tests/unitary/pool/admin/test_revert_commit_params.py b/tests/unitary/pool/admin/test_revert_commit_params.py
new file mode 100644
index 00000000..48228166
--- /dev/null
+++ b/tests/unitary/pool/admin/test_revert_commit_params.py
@@ -0,0 +1,81 @@
+import copy
+
+import boa
+
+
+def _apply_new_params(swap, params):
+    swap.apply_new_parameters(
+        params["mid_fee"],
+        params["out_fee"],
+        params["fee_gamma"],
+        params["allowed_extra_profit"],
+        params["adjustment_step"],
+        params["ma_time"],
+        params["xcp_ma_time"],
+    )
+
+
+def test_commit_incorrect_fee_params(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["mid_fee"] = p["out_fee"] + 1
+    with boa.env.prank(factory_admin):
+        with boa.reverts("mid-fee is too high"):
+            _apply_new_params(swap, p)
+
+        p["out_fee"] = 0
+        with boa.reverts("fee is out of range"):
+            _apply_new_params(swap, p)
+
+        # too large out_fee revert to old out_fee:
+        p["mid_fee"] = params["mid_fee"]
+        p["out_fee"] = 10**10 + 1  # <-- MAX_FEE
+        _apply_new_params(swap, p)
+        logs = swap.get_logs()[0]
+        assert logs.args[1] == params["out_fee"]
+
+
+def test_commit_incorrect_fee_gamma(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["fee_gamma"] = 0
+
+    with boa.env.prank(factory_admin):
+
+        with boa.reverts("fee_gamma out of range [1 .. 10**18]"):
+            _apply_new_params(swap, p)
+
+        p["fee_gamma"] = 10**18 + 1
+        _apply_new_params(swap, p)
+
+    # it will not change fee_gamma as it is above 10**18
+    assert swap.get_logs()[0].args[2] == params["fee_gamma"]
+
+
+def test_commit_rebalancing_params(swap, factory_admin, params):
+
+    p = copy.deepcopy(params)
+    p["allowed_extra_profit"] = 10**18 + 1
+    p["adjustment_step"] == 10**18 + 1
+    p["ma_time"] = 872542 + 1
+
+    with boa.env.prank(factory_admin):
+
+        with boa.env.anchor():
+            _apply_new_params(swap, p)
+            logs = swap.get_logs()[0]
+
+            # values revert to contract's storage values:
+            assert logs.args[3] == params["allowed_extra_profit"]
+            assert logs.args[4] == params["adjustment_step"]
+            assert logs.args[5] == params["ma_time"]
+
+        with boa.reverts("MA time should be longer than 60/ln(2)"):
+            p["ma_time"] = 86
+            _apply_new_params(swap, p)
+
+
+def test_revert_unauthorised_commit(swap, user, params):
+
+    with boa.env.prank(user), boa.reverts(dev="only owner"):
+        _apply_new_params(swap, params)
diff --git a/tests/unitary/pool/admin/test_revert_ramp.py b/tests/unitary/pool/admin/test_revert_ramp.py
new file mode 100644
index 00000000..e485c628
--- /dev/null
+++ b/tests/unitary/pool/admin/test_revert_ramp.py
@@ -0,0 +1,41 @@
+import boa
+
+
+def test_revert_unauthorised_ramp(swap, user):
+
+    with boa.env.prank(user), boa.reverts(dev="only owner"):
+        swap.ramp_A_gamma(1, 1, 1)
+
+
+def test_revert_ramp_while_ramping(swap, factory_admin):
+
+    assert swap.initial_A_gamma_time() == 0
+
+    A_gamma = [swap.A(), swap.gamma()]
+    future_time = boa.env.vm.state.timestamp + 86400 + 1
+    with boa.env.prank(factory_admin):
+        swap.ramp_A_gamma(A_gamma[0] + 1, A_gamma[1] + 1, future_time)
+
+        with boa.reverts(dev="ramp undergoing"):
+            swap.ramp_A_gamma(A_gamma[0], A_gamma[1], future_time)
+
+
+def test_revert_fast_ramps(swap, factory_admin):
+
+    A_gamma = [swap.A(), swap.gamma()]
+    future_time = boa.env.vm.state.timestamp + 10
+    with boa.env.prank(factory_admin), boa.reverts(dev="insufficient time"):
+        swap.ramp_A_gamma(A_gamma[0] + 1, A_gamma[1] + 1, future_time)
+
+
+def test_revert_unauthorised_stop_ramp(swap, factory_admin, user):
+
+    assert swap.initial_A_gamma_time() == 0
+
+    A_gamma = [swap.A(), swap.gamma()]
+    future_time = boa.env.vm.state.timestamp + 86400 + 1
+    with boa.env.prank(factory_admin):
+        swap.ramp_A_gamma(A_gamma[0] + 1, A_gamma[1] + 1, future_time)
+
+    with boa.env.prank(user), boa.reverts(dev="only owner"):
+        swap.stop_ramp_A_gamma()
diff --git a/tests/unitary/pool/stateful/stateful_base.py b/tests/unitary/pool/stateful/stateful_base.py
new file mode 100644
index 00000000..31900774
--- /dev/null
+++ b/tests/unitary/pool/stateful/stateful_base.py
@@ -0,0 +1,270 @@
+import contextlib
+from math import log
+
+import boa
+from boa.test import strategy
+from hypothesis.stateful import RuleBasedStateMachine, invariant, rule
+
+from tests.fixtures.pool import INITIAL_PRICES
+from tests.utils.tokens import mint_for_testing
+
+MAX_SAMPLES = 20
+MAX_D = 10**12 * 10**18  # $1T is hopefully a reasonable cap for tests
+
+
+class StatefulBase(RuleBasedStateMachine):
+    exchange_amount_in = strategy("uint256", max_value=10**9 * 10**18)
+    exchange_i = strategy("uint8", max_value=1)
+    sleep_time = strategy("uint256", max_value=86400 * 7)
+    user = strategy("address")
+
+    def __init__(self):
+
+        super().__init__()
+
+        self.decimals = [int(c.decimals()) for c in self.coins]
+        self.user_balances = {u: [0] * 2 for u in self.users}
+        self.initial_prices = INITIAL_PRICES
+        self.initial_deposit = [
+            10**4 * 10 ** (18 + d) // p
+            for p, d in zip(self.initial_prices, self.decimals)
+        ]  # $10k * 2
+
+        self.xcp_profit = 10**18
+        self.xcp_profit_a = 10**18
+
+        self.total_supply = 0
+        self.previous_pool_profit = 0
+
+        self.swap_admin = self.swap.admin()
+        self.fee_receiver = self.swap.fee_receiver()
+
+        for user in self.users:
+            for coin in self.coins:
+                coin.approve(self.swap, 2**256 - 1, sender=user)
+
+        self.setup()
+
+    def setup(self, user_id=0):
+
+        user = self.users[user_id]
+        for coin, q in zip(self.coins, self.initial_deposit):
+            mint_for_testing(coin, user, q)
+
+        # Very first deposit
+        self.swap.add_liquidity(self.initial_deposit, 0, sender=user)
+
+        self.balances = self.initial_deposit[:]
+        self.total_supply = self.swap.balanceOf(user)
+        self.xcp_profit = 10**18
+
+    def convert_amounts(self, amounts):
+        prices = [10**18] + [self.swap.price_scale()]
+        return [
+            p * a // 10 ** (36 - d)
+            for p, a, d in zip(prices, amounts, self.decimals)
+        ]
+
+    def check_limits(self, amounts, D=True, y=True):
+        """
+        Should be good if within limits, but if outside - can be either
+        """
+        _D = self.swap.D()
+        prices = [10**18] + [self.swap.price_scale()]
+        xp_0 = [self.swap.balances(i) for i in range(2)]
+        xp = xp_0
+        xp_0 = [
+            x * p // 10**d for x, p, d in zip(xp_0, prices, self.decimals)
+        ]
+        xp = [
+            (x + a) * p // 10**d
+            for a, x, p, d in zip(amounts, xp, prices, self.decimals)
+        ]
+
+        if D:
+            for _xp in [xp_0, xp]:
+                if (
+                    (min(_xp) * 10**18 // max(_xp) < 10**14)
+                    or (max(_xp) < 10**9 * 10**18)
+                    or (max(_xp) > 10**15 * 10**18)
+                ):
+                    return False
+
+        if y:
+            for _xp in [xp_0, xp]:
+                if (
+                    (_D < 10**17)
+                    or (_D > 10**15 * 10**18)
+                    or (min(_xp) * 10**18 // _D < 10**16)
+                    or (max(_xp) * 10**18 // _D > 10**20)
+                ):
+                    return False
+
+        return True
+
+    @rule(
+        exchange_amount_in=exchange_amount_in,
+        exchange_i=exchange_i,
+        user=user,
+    )
+    def exchange(self, exchange_amount_in, exchange_i, user):
+        out = self._exchange(exchange_amount_in, exchange_i, user)
+        if out:
+            self.swap_out = out
+            return
+        self.swap_out = None
+
+    def _exchange(
+        self, exchange_amount_in, exchange_i, user, check_out_amount=True
+    ):
+        exchange_j = 1 - exchange_i
+        try:
+            calc_amount = self.swap.get_dy(
+                exchange_i, exchange_j, exchange_amount_in
+            )
+        except Exception:
+            _amounts = [0] * 2
+            _amounts[exchange_i] = exchange_amount_in
+            if self.check_limits(_amounts) and exchange_amount_in > 10000:
+                raise
+            return None
+        mint_for_testing(self.coins[exchange_i], user, exchange_amount_in)
+
+        d_balance_i = self.coins[exchange_i].balanceOf(user)
+        d_balance_j = self.coins[exchange_j].balanceOf(user)
+        try:
+            self.coins[exchange_i].approve(
+                self.swap, 2**256 - 1, sender=user
+            )
+            out = self.swap.exchange(
+                exchange_i, exchange_j, exchange_amount_in, 0, sender=user
+            )
+        except Exception:
+            # Small amounts may fail with rounding errors
+            if (
+                calc_amount > 100
+                and exchange_amount_in > 100
+                and calc_amount / self.swap.balances(exchange_j) > 1e-13
+                and exchange_amount_in / self.swap.balances(exchange_i) > 1e-13
+            ):
+                raise
+            return None
+
+        # This is to check that we didn't end up in a borked state after
+        # an exchange succeeded
+        self.swap.get_dy(
+            exchange_j,
+            exchange_i,
+            10**16
+            * 10 ** self.decimals[exchange_j]
+            // INITIAL_PRICES[exchange_j],
+        )
+
+        d_balance_i -= self.coins[exchange_i].balanceOf(user)
+        d_balance_j -= self.coins[exchange_j].balanceOf(user)
+
+        assert d_balance_i == exchange_amount_in
+        if check_out_amount:
+            if check_out_amount is True:
+                assert (
+                    -d_balance_j == calc_amount
+                ), f"{-d_balance_j} vs {calc_amount}"
+            else:
+                assert abs(d_balance_j + calc_amount) < max(
+                    check_out_amount * calc_amount, 3
+                ), f"{-d_balance_j} vs {calc_amount}"
+
+        self.balances[exchange_i] += d_balance_i
+        self.balances[exchange_j] += d_balance_j
+
+        return out
+
+    @rule(sleep_time=sleep_time)
+    def sleep(self, sleep_time):
+        boa.env.time_travel(sleep_time)
+
+    @invariant()
+    def balances(self):
+        balances = [self.swap.balances(i) for i in range(2)]
+        balances_of = [c.balanceOf(self.swap) for c in self.coins]
+        for i in range(2):
+            assert self.balances[i] == balances[i]
+            assert self.balances[i] == balances_of[i]
+
+    @invariant()
+    def total_supply(self):
+        assert self.total_supply == self.swap.totalSupply()
+
+    @invariant()
+    def virtual_price(self):
+        virtual_price = self.swap.virtual_price()
+        xcp_profit = self.swap.xcp_profit()
+        get_virtual_price = self.swap.get_virtual_price()
+
+        assert xcp_profit >= 10**18 - 10
+        assert virtual_price >= 10**18 - 10
+        assert get_virtual_price >= 10**18 - 10
+
+        assert (
+            xcp_profit - self.xcp_profit > -3
+        ), f"{xcp_profit} vs {self.xcp_profit}"
+        assert (virtual_price - 10**18) * 2 - (
+            xcp_profit - 10**18
+        ) >= -5, f"vprice={virtual_price}, xcp_profit={xcp_profit}"
+        assert abs(log(virtual_price / get_virtual_price)) < 1e-10
+
+        self.xcp_profit = xcp_profit
+
+    @invariant()
+    def up_only_profit(self):
+
+        current_profit = xcp_profit = self.swap.xcp_profit()
+        xcp_profit_a = self.swap.xcp_profit_a()
+        current_profit = (xcp_profit + xcp_profit_a + 1) // 2
+
+        assert current_profit >= self.previous_pool_profit
+        self.previous_pool_profit = current_profit
+
+    @contextlib.contextmanager
+    def upkeep_on_claim(self):
+
+        admin_balances_pre = [
+            c.balanceOf(self.fee_receiver) for c in self.coins
+        ]
+        pool_is_ramping = (
+            self.swap.future_A_gamma_time() > boa.env.vm.state.timestamp
+        )
+
+        try:
+
+            yield
+
+        finally:
+
+            new_xcp_profit_a = self.swap.xcp_profit_a()
+            old_xcp_profit_a = self.xcp_profit_a
+
+            claimed = False
+            if new_xcp_profit_a > old_xcp_profit_a:
+                claimed = True
+                self.xcp_profit_a = new_xcp_profit_a
+
+            admin_balances_post = [
+                c.balanceOf(self.fee_receiver) for c in self.coins
+            ]
+
+            if claimed:
+
+                for i in range(2):
+                    claimed_amount = (
+                        admin_balances_post[i] - admin_balances_pre[i]
+                    )
+                    assert (
+                        claimed_amount > 0
+                    )  # check if non zero amounts of claim
+                    assert not pool_is_ramping  # cannot claim while ramping
+
+                    # update self.balances
+                    self.balances[i] -= claimed_amount
+
+        self.xcp_profit = self.swap.xcp_profit()
diff --git a/tests/unitary/pool/stateful/test_multiprecision.py b/tests/unitary/pool/stateful/test_multiprecision.py
new file mode 100644
index 00000000..dcecf832
--- /dev/null
+++ b/tests/unitary/pool/stateful/test_multiprecision.py
@@ -0,0 +1,54 @@
+import pytest
+from boa.test import strategy
+from hypothesis.stateful import rule, run_state_machine_as_test
+
+from tests.unitary.pool.stateful.test_stateful import NumbaGoUp
+
+MAX_SAMPLES = 100
+STEP_COUNT = 100
+
+
+@pytest.fixture(scope="module")
+def coins(stgusdc):
+    return stgusdc
+
+
+@pytest.fixture(scope="module")
+def swap(swap_multiprecision):
+    return swap_multiprecision
+
+
+class Multiprecision(NumbaGoUp):
+    exchange_amount_in = strategy(
+        "uint256", min_value=10**18, max_value=50000 * 10**18
+    )
+    user = strategy("address")
+    exchange_i = strategy("uint8", max_value=1)
+
+    @rule(
+        exchange_amount_in=exchange_amount_in,
+        exchange_i=exchange_i,
+        user=user,
+    )
+    def exchange(self, exchange_amount_in, exchange_i, user):
+        exchange_amount_in = exchange_amount_in // 10 ** (
+            18 - self.decimals[exchange_i]
+        )
+        super().exchange(exchange_amount_in, exchange_i, user)
+
+
+def test_multiprecision(users, coins, swap):
+    from hypothesis import settings
+    from hypothesis._settings import HealthCheck
+
+    Multiprecision.TestCase.settings = settings(
+        max_examples=MAX_SAMPLES,
+        stateful_step_count=STEP_COUNT,
+        suppress_health_check=HealthCheck.all(),
+        deadline=None,
+    )
+
+    for k, v in locals().items():
+        setattr(Multiprecision, k, v)
+
+    run_state_machine_as_test(Multiprecision)
diff --git a/tests/unitary/pool/stateful/test_ramp.py b/tests/unitary/pool/stateful/test_ramp.py
new file mode 100644
index 00000000..ab62ac9a
--- /dev/null
+++ b/tests/unitary/pool/stateful/test_ramp.py
@@ -0,0 +1,107 @@
+import boa
+from boa.test import strategy
+from hypothesis.stateful import invariant, rule, run_state_machine_as_test
+
+from tests.unitary.pool.stateful.test_stateful import NumbaGoUp
+
+MAX_SAMPLES = 20
+STEP_COUNT = 100
+MAX_D = 10**12 * 10**18  # $1T is hopefully a reasonable cap for tests
+ALLOWED_DIFFERENCE = 0.001
+
+
+class RampTest(NumbaGoUp):
+    check_out_amount = strategy("bool")
+    exchange_amount_in = strategy(
+        "uint256", min_value=10**18, max_value=50000 * 10**18
+    )
+    token_amount = strategy(
+        "uint256", min_value=10**18, max_value=10**12 * 10**18
+    )
+    deposit_amounts = strategy(
+        "uint256[3]", min_value=10**18, max_value=10**9 * 10**18
+    )
+    user = strategy("address")
+    exchange_i = strategy("uint8", max_value=1)
+
+    def setup(self, user_id=0):
+        super().setup(user_id)
+        new_A = self.swap.A() * 2
+        new_gamma = self.swap.gamma() * 2
+        self.swap.ramp_A_gamma(
+            new_A,
+            new_gamma,
+            boa.env.vm.state.timestamp + 14 * 86400,
+            sender=self.swap_admin,
+        )
+
+    @rule(user=user, deposit_amounts=deposit_amounts)
+    def deposit(self, deposit_amounts, user):
+        deposit_amounts[1:] = [
+            deposit_amounts[0],
+            deposit_amounts[1] * 10**18 // self.swap.price_oracle(),
+        ]
+        super().deposit(deposit_amounts, user)
+
+    @rule(
+        user=user,
+        exchange_i=exchange_i,
+        exchange_amount_in=exchange_amount_in,
+        check_out_amount=check_out_amount,
+    )
+    def exchange(self, exchange_amount_in, exchange_i, user, check_out_amount):
+
+        if exchange_i > 0:
+            exchange_amount_in = (
+                exchange_amount_in * 10**18 // self.swap.price_oracle()
+            )
+            if exchange_amount_in < 1000:
+                return
+
+        super()._exchange(
+            exchange_amount_in,
+            exchange_i,
+            user,
+            ALLOWED_DIFFERENCE if check_out_amount else False,
+        )
+
+    @rule(
+        user=user,
+        token_amount=token_amount,
+        exchange_i=exchange_i,
+        check_out_amount=check_out_amount,
+    )
+    def remove_liquidity_one_coin(
+        self, token_amount, exchange_i, user, check_out_amount
+    ):
+
+        if check_out_amount:
+            super().remove_liquidity_one_coin(
+                token_amount, exchange_i, user, ALLOWED_DIFFERENCE
+            )
+        else:
+            super().remove_liquidity_one_coin(
+                token_amount, exchange_i, user, False
+            )
+
+    @invariant()
+    def virtual_price(self):
+        # Invariant is not conserved here
+        pass
+
+
+def test_ramp(users, coins, swap):
+    from hypothesis import settings
+    from hypothesis._settings import HealthCheck
+
+    RampTest.TestCase.settings = settings(
+        max_examples=MAX_SAMPLES,
+        stateful_step_count=STEP_COUNT,
+        suppress_health_check=HealthCheck.all(),
+        deadline=None,
+    )
+
+    for k, v in locals().items():
+        setattr(RampTest, k, v)
+
+    run_state_machine_as_test(RampTest)
diff --git a/tests/unitary/pool/stateful/test_ramp_nocheck.py b/tests/unitary/pool/stateful/test_ramp_nocheck.py
new file mode 100644
index 00000000..2519fbc4
--- /dev/null
+++ b/tests/unitary/pool/stateful/test_ramp_nocheck.py
@@ -0,0 +1,81 @@
+import boa
+from boa.test import strategy
+from hypothesis.stateful import invariant, rule, run_state_machine_as_test
+
+from tests.unitary.pool.stateful.test_stateful import NumbaGoUp
+
+MAX_SAMPLES = 20
+STEP_COUNT = 100
+MAX_D = 10**12 * 10**18  # $1T is hopefully a reasonable cap for tests
+ALLOWED_DIFFERENCE = 0.02
+
+
+class RampTest(NumbaGoUp):
+    future_gamma = strategy(
+        "uint256",
+        min_value=int(2.8e-4 * 1e18 / 9),
+        max_value=int(2.8e-4 * 1e18 * 9),
+    )
+    future_A = strategy(
+        "uint256",
+        min_value=90 * 2**2 * 10000 // 9,
+        max_value=90 * 2**2 * 10000 * 9,
+    )
+    check_out_amount = strategy("bool")
+    exchange_amount_in = strategy(
+        "uint256", min_value=10**18, max_value=50000 * 10**18
+    )
+    token_amount = strategy(
+        "uint256", min_value=10**18, max_value=10**12 * 10**18
+    )
+    user = strategy("address")
+    exchange_i = strategy("uint8", max_value=1)
+
+    def initialize(self, future_A, future_gamma):
+        self.swap.ramp_A_gamma(
+            future_A,
+            future_gamma,
+            boa.env.vm.state.timestamp + 14 * 86400,
+            sender=self.swap_admin,
+        )
+
+    @rule(
+        exchange_amount_in=exchange_amount_in,
+        exchange_i=exchange_i,
+        user=user,
+    )
+    def exchange(self, exchange_amount_in, exchange_i, user):
+        try:
+            super()._exchange(exchange_amount_in, exchange_i, user, False)
+        except Exception:
+            if exchange_amount_in > 10**9:
+                # Small swaps can fail at ramps
+                raise
+
+    @rule(token_amount=token_amount, exchange_i=exchange_i, user=user)
+    def remove_liquidity_one_coin(self, token_amount, exchange_i, user):
+        super().remove_liquidity_one_coin(
+            token_amount, exchange_i, user, False
+        )
+
+    @invariant()
+    def virtual_price(self):
+        # Invariant is not conserved here
+        pass
+
+
+def test_ramp(users, coins, swap):
+    from hypothesis import settings
+    from hypothesis._settings import HealthCheck
+
+    RampTest.TestCase.settings = settings(
+        max_examples=MAX_SAMPLES,
+        stateful_step_count=STEP_COUNT,
+        suppress_health_check=HealthCheck.all(),
+        deadline=None,
+    )
+
+    for k, v in locals().items():
+        setattr(RampTest, k, v)
+
+    run_state_machine_as_test(RampTest)
diff --git a/tests/unitary/pool/stateful/test_simulate.py b/tests/unitary/pool/stateful/test_simulate.py
new file mode 100644
index 00000000..3b86f37b
--- /dev/null
+++ b/tests/unitary/pool/stateful/test_simulate.py
@@ -0,0 +1,120 @@
+import boa
+from boa.test import strategy
+from hypothesis.stateful import invariant, rule, run_state_machine_as_test
+
+from tests.unitary.pool.stateful.stateful_base import StatefulBase
+from tests.utils import approx
+from tests.utils import simulation_int_many as sim
+from tests.utils.tokens import mint_for_testing
+
+MAX_SAMPLES = 20
+STEP_COUNT = 100
+
+
+class StatefulSimulation(StatefulBase):
+    exchange_amount_in = strategy(
+        "uint256", min_value=10**17, max_value=10**5 * 10**18
+    )
+    exchange_i = strategy("uint8", max_value=1)
+    user = strategy("address")
+
+    def setup(self):
+
+        super().setup()
+
+        for u in self.users[1:]:
+            for coin, q in zip(self.coins, self.initial_deposit):
+                mint_for_testing(coin, u, q)
+            for i in range(2):
+                self.balances[i] += self.initial_deposit[i]
+            self.swap.add_liquidity(self.initial_deposit, 0, sender=u)
+            self.total_supply += self.swap.balanceOf(u)
+
+        self.virtual_price = self.swap.get_virtual_price()
+
+        self.trader = sim.Trader(
+            self.swap.A(),
+            self.swap.gamma(),
+            self.swap.D(),
+            2,
+            [10**18, self.swap.price_scale()],
+            self.swap.mid_fee() / 1e10,
+            self.swap.out_fee() / 1e10,
+            self.swap.allowed_extra_profit(),
+            self.swap.fee_gamma(),
+            self.swap.adjustment_step() / 1e18,
+            int(
+                self.swap.ma_time() / 0.693
+            ),  # crypto swap returns ma time in sec
+        )
+        for i in range(2):
+            self.trader.curve.x[i] = self.swap.balances(i)
+
+        # Adjust virtual prices
+        self.trader.xcp_profit = self.swap.xcp_profit()
+        self.trader.xcp_profit_real = self.swap.virtual_price()
+        self.trader.t = boa.env.vm.state.timestamp
+        self.swap_no = 0
+
+    @rule(
+        exchange_amount_in=exchange_amount_in,
+        exchange_i=exchange_i,
+        user=user,
+    )
+    def exchange(self, exchange_amount_in, exchange_i, user):
+
+        dx = (
+            exchange_amount_in
+            * 10**18
+            // self.trader.price_oracle[exchange_i]
+        )
+        self.swap_no += 1
+        super().exchange(dx, exchange_i, user)
+
+        if not self.swap_out:
+            return  # if swap breaks, dont check.
+
+        dy_trader = self.trader.buy(dx, exchange_i, 1 - exchange_i)
+        self.trader.tweak_price(boa.env.vm.state.timestamp)
+
+        # exchange checks:
+        assert approx(self.swap_out, dy_trader, 1e-3)
+        assert approx(
+            self.swap.price_oracle(), self.trader.price_oracle[1], 1e-3
+        )
+
+        boa.env.time_travel(12)
+
+    @invariant()
+    def simulator(self):
+        if self.trader.xcp_profit / 1e18 - 1 > 1e-8:
+            assert (
+                abs(self.trader.xcp_profit - self.swap.xcp_profit())
+                / (self.trader.xcp_profit - 10**18)
+                < 0.05
+            )
+
+        price_scale = self.swap.price_scale()
+        price_trader = self.trader.curve.p[1]
+        try:
+            assert approx(price_scale, price_trader, 1e-3)
+        except Exception:
+            if self.check_limits([0, 0, 0]):
+                assert False
+
+
+def test_sim(users, coins, swap):
+    from hypothesis import settings
+    from hypothesis._settings import HealthCheck
+
+    StatefulSimulation.TestCase.settings = settings(
+        max_examples=MAX_SAMPLES,
+        stateful_step_count=STEP_COUNT,
+        suppress_health_check=HealthCheck.all(),
+        deadline=None,
+    )
+
+    for k, v in locals().items():
+        setattr(StatefulSimulation, k, v)
+
+    run_state_machine_as_test(StatefulSimulation)
diff --git a/tests/unitary/pool/stateful/test_stateful.py b/tests/unitary/pool/stateful/test_stateful.py
new file mode 100644
index 00000000..8f9653ce
--- /dev/null
+++ b/tests/unitary/pool/stateful/test_stateful.py
@@ -0,0 +1,182 @@
+import boa
+from boa.test import strategy
+from hypothesis.stateful import rule, run_state_machine_as_test
+
+from tests.fixtures.pool import INITIAL_PRICES
+from tests.unitary.pool.stateful.stateful_base import StatefulBase
+from tests.utils.tokens import mint_for_testing
+
+MAX_SAMPLES = 20
+STEP_COUNT = 100
+MAX_D = 10**12 * 10**18  # $1T is hopefully a reasonable cap for tests
+
+
+class NumbaGoUp(StatefulBase):
+    """
+    Test that profit goes up
+    """
+
+    user = strategy("address")
+    exchange_i = strategy("uint8", max_value=1)
+    deposit_amounts = strategy(
+        "uint256[2]", min_value=0, max_value=10**9 * 10**18
+    )
+    token_amount = strategy("uint256", max_value=10**12 * 10**18)
+    check_out_amount = strategy("bool")
+
+    @rule(deposit_amounts=deposit_amounts, user=user)
+    def deposit(self, deposit_amounts, user):
+
+        if self.swap.D() > MAX_D:
+            return
+
+        amounts = self.convert_amounts(deposit_amounts)
+        if sum(amounts) == 0:
+            return
+
+        new_balances = [x + y for x, y in zip(self.balances, amounts)]
+
+        for coin, q in zip(self.coins, amounts):
+            mint_for_testing(coin, user, q)
+
+        try:
+
+            tokens = self.swap.balanceOf(user)
+            self.swap.add_liquidity(amounts, 0, sender=user)
+            tokens = self.swap.balanceOf(user) - tokens
+            self.total_supply += tokens
+            self.balances = new_balances
+
+        except Exception:
+
+            if self.check_limits(amounts):
+                raise
+            return
+
+        # This is to check that we didn't end up in a borked state after
+        # an exchange succeeded
+        try:
+            self.swap.get_dy(0, 1, 10 ** (self.decimals[0] - 2))
+        except Exception:
+            self.swap.get_dy(
+                1,
+                0,
+                10**16 * 10 ** self.decimals[1] // self.swap.price_scale(),
+            )
+
+    @rule(token_amount=token_amount, user=user)
+    def remove_liquidity(self, token_amount, user):
+        if self.swap.balanceOf(user) < token_amount or token_amount == 0:
+            with boa.reverts():
+                self.swap.remove_liquidity(token_amount, [0] * 2, sender=user)
+        else:
+            amounts = [c.balanceOf(user) for c in self.coins]
+            tokens = self.swap.balanceOf(user)
+            with self.upkeep_on_claim():
+                self.swap.remove_liquidity(token_amount, [0] * 2, sender=user)
+            tokens -= self.swap.balanceOf(user)
+            self.total_supply -= tokens
+            amounts = [
+                (c.balanceOf(user) - a) for c, a in zip(self.coins, amounts)
+            ]
+            self.balances = [b - a for a, b in zip(amounts, self.balances)]
+
+            # Virtual price resets if everything is withdrawn
+            if self.total_supply == 0:
+                self.virtual_price = 10**18
+
+    @rule(
+        token_amount=token_amount,
+        exchange_i=exchange_i,
+        user=user,
+        check_out_amount=check_out_amount,
+    )
+    def remove_liquidity_one_coin(
+        self, token_amount, exchange_i, user, check_out_amount
+    ):
+
+        try:
+            calc_out_amount = self.swap.calc_withdraw_one_coin(
+                token_amount, exchange_i
+            )
+        except Exception:
+            if (
+                self.check_limits([0] * 2)
+                and not (token_amount > self.total_supply)
+                and token_amount > 10000
+            ):
+                self.swap.calc_withdraw_one_coin(
+                    token_amount, exchange_i, sender=user
+                )
+            return
+
+        d_token = self.swap.balanceOf(user)
+        if d_token < token_amount:
+            with boa.reverts():
+                self.swap.remove_liquidity_one_coin(
+                    token_amount, exchange_i, 0, sender=user
+                )
+            return
+
+        d_balance = self.coins[exchange_i].balanceOf(user)
+        try:
+            with self.upkeep_on_claim():
+                self.swap.remove_liquidity_one_coin(
+                    token_amount, exchange_i, 0, sender=user
+                )
+        except Exception:
+            # Small amounts may fail with rounding errors
+            if (
+                calc_out_amount > 100
+                and token_amount / self.total_supply > 1e-10
+                and calc_out_amount / self.swap.balances(exchange_i) > 1e-10
+            ):
+                raise
+            return
+
+        # This is to check that we didn't end up in a borked state after
+        # an exchange succeeded
+        _deposit = [0] * 2
+        _deposit[exchange_i] = (
+            10**16
+            * 10 ** self.decimals[exchange_i]
+            // ([10**18] + INITIAL_PRICES)[exchange_i]
+        )
+        assert self.swap.calc_token_amount(_deposit, True)
+
+        d_balance = self.coins[exchange_i].balanceOf(user) - d_balance
+        d_token = d_token - self.swap.balanceOf(user)
+
+        if check_out_amount:
+            if check_out_amount is True:
+                assert (
+                    calc_out_amount == d_balance
+                ), f"{calc_out_amount} vs {d_balance} for {token_amount}"
+            else:
+                assert abs(calc_out_amount - d_balance) <= max(
+                    check_out_amount * calc_out_amount, 5
+                ), f"{calc_out_amount} vs {d_balance} for {token_amount}"
+
+        self.balances[exchange_i] -= d_balance
+        self.total_supply -= d_token
+
+        # Virtual price resets if everything is withdrawn
+        if self.total_supply == 0:
+            self.virtual_price = 10**18
+
+
+def test_numba_go_up(users, coins, swap):
+    from hypothesis import settings
+    from hypothesis._settings import HealthCheck
+
+    NumbaGoUp.TestCase.settings = settings(
+        max_examples=MAX_SAMPLES,
+        stateful_step_count=STEP_COUNT,
+        suppress_health_check=HealthCheck.all(),
+        deadline=None,
+    )
+
+    for k, v in locals().items():
+        setattr(NumbaGoUp, k, v)
+
+    run_state_machine_as_test(NumbaGoUp)
diff --git a/tests/unitary/pool/test_a_gamma.py b/tests/unitary/pool/test_a_gamma.py
index 07966400..33167a5a 100644
--- a/tests/unitary/pool/test_a_gamma.py
+++ b/tests/unitary/pool/test_a_gamma.py
@@ -1,10 +1,7 @@
 import boa
 
 
-def test_A_gamma(swap, math_contract, params):
-
-    breakpoint()
-
+def test_A_gamma(swap, params):
     A = swap.A()
     gamma = swap.gamma()
 
diff --git a/tests/unitary/pool/test_admin_fee.py b/tests/unitary/pool/test_admin_fee.py
new file mode 100644
index 00000000..a314907d
--- /dev/null
+++ b/tests/unitary/pool/test_admin_fee.py
@@ -0,0 +1,45 @@
+import boa
+from hypothesis import given, settings
+from hypothesis import strategies as st
+
+from tests.fixtures.pool import INITIAL_PRICES
+from tests.utils.tokens import mint_for_testing
+
+
+@given(ratio=st.floats(min_value=0.0001, max_value=0.1))
+@settings(max_examples=10, deadline=None)
+def test_admin_fee_after_deposit(
+    swap, coins, fee_receiver, user, user_b, ratio
+):
+    quantities = [10**42 // p for p in INITIAL_PRICES]
+
+    for coin, q in zip(coins, quantities):
+        for u in [user, user_b]:
+            mint_for_testing(coin, u, q)
+            with boa.env.prank(u):
+                coin.approve(swap, 2**256 - 1)
+
+    split_quantities = [quantities[0] // 100, quantities[1] // 100]
+
+    with boa.env.prank(user):
+        swap.add_liquidity(split_quantities, 0)
+
+    with boa.env.prank(user_b):
+        for _ in range(100):
+            before = coins[1].balanceOf(user_b)
+            swap.exchange(0, 1, split_quantities[0] // 100, 0)
+            after = coins[1].balanceOf(user_b)
+            to_swap = after - before
+            swap.exchange(1, 0, to_swap, 0)
+
+    balances = [swap.balances(i) for i in range(2)]
+    split_quantities = [int(balances[0] * ratio), int(balances[1] * ratio)]
+    with boa.env.prank(user):
+        swap.add_liquidity(split_quantities, 0)
+
+    assert (
+        coins[0].balanceOf(fee_receiver) + coins[1].balanceOf(fee_receiver)
+        == 0
+    )
+
+    return swap
diff --git a/tests/unitary/pool/test_deposit_withdraw.py b/tests/unitary/pool/test_deposit_withdraw.py
new file mode 100644
index 00000000..b70bffee
--- /dev/null
+++ b/tests/unitary/pool/test_deposit_withdraw.py
@@ -0,0 +1,373 @@
+import boa
+import pytest
+from boa.test import strategy
+from hypothesis import given, settings
+
+from tests.fixtures.pool import INITIAL_PRICES
+from tests.utils import approx
+from tests.utils import simulation_int_many as sim
+from tests.utils.tokens import mint_for_testing
+
+SETTINGS = {"max_examples": 100, "deadline": None}
+
+
+def assert_string_contains(string, substrings):
+    assert any(substring in string for substring in substrings)
+
+
+@pytest.fixture(scope="module")
+def test_1st_deposit_and_last_withdraw(swap, coins, user, fee_receiver):
+
+    quantities = [10**36 // p for p in INITIAL_PRICES]  # $3M worth
+
+    for coin, q in zip(coins, quantities):
+        mint_for_testing(coin, user, q)
+        with boa.env.prank(user):
+            coin.approve(swap, 2**256 - 1)
+
+    bal_before = boa.env.get_balance(swap.address)
+    with boa.env.prank(user):
+        swap.add_liquidity(quantities, 0)
+
+    # test if eth wasnt deposited:
+    assert boa.env.get_balance(swap.address) == bal_before
+
+    token_balance = swap.balanceOf(user)
+    assert (
+        token_balance == swap.totalSupply() - swap.balanceOf(fee_receiver) > 0
+    )
+    assert abs(swap.get_virtual_price() / 1e18 - 1) < 1e-3
+
+    # Empty the contract
+    with boa.env.prank(user):
+        swap.remove_liquidity(token_balance, [0] * len(coins))
+
+    # nothing is left except admin balances
+    assert (
+        swap.balanceOf(user)
+        == swap.totalSupply() - swap.balanceOf(fee_receiver)
+        == 0
+    )
+
+    return swap
+
+
+@pytest.fixture(scope="module")
+def test_first_deposit_full_withdraw_second_deposit(
+    test_1st_deposit_and_last_withdraw, user, coins, fee_receiver
+):
+    swap = test_1st_deposit_and_last_withdraw
+    quantities = [10**36 // p for p in INITIAL_PRICES]  # $2M worth
+
+    for coin, q in zip(coins, quantities):
+        mint_for_testing(coin, user, q)
+        with boa.env.prank(user):
+            coin.approve(swap, 2**256 - 1)
+
+    # Second deposit
+    eth_bal_before = boa.env.get_balance(swap.address)
+    swap_balances_before = [swap.balances(i) for i in range(3)]
+    with boa.env.prank(user):
+        swap.add_liquidity(quantities, 0)
+
+    assert swap.xcp_profit_a() >= 10**18
+    assert swap.xcp_profit() >= 10**18
+    assert swap.virtual_price() >= 10**18
+
+    # test if eth was not deposited:
+    assert boa.env.get_balance(swap.address) == eth_bal_before
+
+    for i in range(len(coins)):
+        assert swap.balances(i) == quantities[i] + swap_balances_before[i]
+
+    token_balance = swap.balanceOf(user)
+    assert (
+        token_balance + swap.balanceOf(fee_receiver) == swap.totalSupply() > 0
+    )
+    assert abs(swap.get_virtual_price() / 1e18 - 1) < 1e-3
+
+    return swap
+
+
+@given(
+    i=strategy("uint", min_value=0, max_value=1),
+    amount=strategy("uint256", max_value=2000, min_value=1),
+)
+@settings(**SETTINGS)
+def test_second_deposit_single_token(
+    swap_with_deposit, coins, user, i, amount
+):
+
+    # deposit single token:
+    quantities = [0, 0]
+    for j in range(len(coins)):
+        if j == i:
+            quantities[j] = amount * 10 ** coins[i].decimals()
+            mint_for_testing(coins[j], user, quantities[j])
+
+    # Single sided deposit
+    with boa.env.prank(user):
+        swap_with_deposit.add_liquidity(quantities, 0)
+
+
+@given(
+    values=strategy(
+        "uint256[2]", min_value=10**16, max_value=10**9 * 10**18
+    )
+)
+@settings(**SETTINGS)
+def test_second_deposit(
+    swap_with_deposit,
+    coins,
+    user,
+    values,
+):
+
+    amounts = [v * 10**18 // p for v, p in zip(values, INITIAL_PRICES)]
+
+    # get simmed D value here:
+    xp = [10**6 * 10**18] * len(coins)  # initial D
+
+    # D after second deposit:
+    for i in range(len(coins)):
+        xp[i] += int(values[i] * 10**18)
+
+    _A, _gamma = [swap_with_deposit.A(), swap_with_deposit.gamma()]
+    _D = sim.solve_D(_A, _gamma, xp)
+
+    safe = all(
+        f >= 1.1e16 and f <= 0.9e20 for f in [_x * 10**18 // _D for _x in xp]
+    )
+
+    for coin, q in zip(coins, amounts):
+        mint_for_testing(coin, user, 10**30)
+        with boa.env.prank(user):
+            coin.approve(swap_with_deposit, 2**256 - 1)
+
+    try:
+
+        calculated = swap_with_deposit.calc_token_amount(amounts, True)
+        measured = swap_with_deposit.balanceOf(user)
+        d_balances = [swap_with_deposit.balances(i) for i in range(len(coins))]
+        claimed_fees = [0] * len(coins)
+
+        with boa.env.prank(user):
+            swap_with_deposit.add_liquidity(amounts, int(calculated * 0.999))
+
+        logs = swap_with_deposit.get_logs()
+        for log in logs:
+            if log.event_type.name == "ClaimAdminFee":
+                claimed_fees = log.args[0]
+
+        d_balances = [
+            swap_with_deposit.balances(i) - d_balances[i] + claimed_fees[i]
+            for i in range(len(coins))
+        ]
+        measured = swap_with_deposit.balanceOf(user) - measured
+
+        assert calculated <= measured
+        assert tuple(amounts) == tuple(d_balances)
+
+    except Exception:
+
+        if safe:
+            raise
+
+    # This is to check that we didn't end up in a borked state after
+    # a deposit succeeded
+    swap_with_deposit.get_dy(0, 1, 10**16)
+
+
+@given(
+    value=strategy(
+        "uint256", min_value=10**16, max_value=10**6 * 10**18
+    ),
+    i=strategy("uint", min_value=0, max_value=1),
+)
+@settings(**SETTINGS)
+def test_second_deposit_one(
+    swap_with_deposit,
+    views_contract,
+    coins,
+    user,
+    value,
+    i,
+):
+
+    amounts = [0] * len(coins)
+    amounts[i] = value * 10**18 // (INITIAL_PRICES)[i]
+    mint_for_testing(coins[i], user, amounts[i])
+
+    try:
+
+        calculated = views_contract.calc_token_amount(
+            amounts, True, swap_with_deposit
+        )
+        measured = swap_with_deposit.balanceOf(user)
+        d_balances = [swap_with_deposit.balances(i) for i in range(len(coins))]
+        claimed_fees = [0] * len(coins)
+
+        with boa.env.prank(user):
+            swap_with_deposit.add_liquidity(amounts, int(calculated * 0.999))
+
+        logs = swap_with_deposit.get_logs()
+        for log in logs:
+            if log.event_type.name == "ClaimAdminFee":
+                claimed_fees = log.args[0]
+
+        d_balances = [
+            swap_with_deposit.balances(i) - d_balances[i] + claimed_fees[i]
+            for i in range(len(coins))
+        ]
+        measured = swap_with_deposit.balanceOf(user) - measured
+
+        assert calculated <= measured
+        assert tuple(amounts) == tuple(d_balances)
+
+    except boa.BoaError as b_error:
+
+        assert_string_contains(
+            b_error.stack_trace.last_frame.pretty_vm_reason,
+            ["Unsafe value for y", "Unsafe values x[i]"],
+        )
+
+
+@given(
+    token_amount=strategy(
+        "uint256", min_value=10**12, max_value=4000 * 10**18
+    )
+)  # supply is 2400 * 1e18
+@settings(**SETTINGS)
+def test_immediate_withdraw(
+    swap_with_deposit,
+    views_contract,
+    coins,
+    user,
+    token_amount,
+):
+
+    f = token_amount / swap_with_deposit.totalSupply()
+    if f <= 1:
+        expected = [
+            int(f * swap_with_deposit.balances(i)) for i in range(len(coins))
+        ]
+        measured = [c.balanceOf(user) for c in coins]
+        token_amount_calc = views_contract.calc_token_amount(
+            expected, False, swap_with_deposit
+        )
+        assert abs(token_amount_calc - token_amount) / token_amount < 1e-3
+        d_balances = [swap_with_deposit.balances(i) for i in range(len(coins))]
+
+        with boa.env.prank(user):
+            swap_with_deposit.remove_liquidity(
+                token_amount,
+                [int(0.999 * e) for e in expected],
+            )
+
+        d_balances = [
+            d_balances[i] - swap_with_deposit.balances(i)
+            for i in range(len(coins))
+        ]
+        measured = [c.balanceOf(user) - m for c, m in zip(coins, measured)]
+
+        for e, m in zip(expected, measured):
+            assert abs(e - m) / e < 1e-3
+
+        assert tuple(d_balances) == tuple(measured)
+
+    else:
+        with boa.reverts(), boa.env.prank(user):
+            swap_with_deposit.remove_liquidity(token_amount, [0] * len(coins))
+
+
+@given(
+    token_amount=strategy(
+        "uint256", min_value=10**12, max_value=4 * 10**6 * 10**18
+    ),  # Can be more than we have
+    i=strategy("uint", min_value=0, max_value=1),
+)
+@settings(**SETTINGS)
+def test_immediate_withdraw_one(
+    swap_with_deposit,
+    views_contract,
+    coins,
+    user,
+    token_amount,
+    i,
+):
+
+    if token_amount >= swap_with_deposit.totalSupply():
+        with boa.reverts():
+            swap_with_deposit.calc_withdraw_one_coin(token_amount, i)
+
+    else:
+
+        # Test if we are safe
+        xp = [10**6 * 10**18] * len(coins)
+        _supply = swap_with_deposit.totalSupply()
+        _A, _gamma = [swap_with_deposit.A(), swap_with_deposit.gamma()]
+        _D = swap_with_deposit.D() * (_supply - token_amount) // _supply
+
+        xp[i] = sim.solve_x(_A, _gamma, xp, _D, i)
+
+        safe = all(
+            f >= 1.1e16 and f <= 0.9e20
+            for f in [_x * 10**18 // _D for _x in xp]
+        )
+
+        try:
+            calculated = swap_with_deposit.calc_withdraw_one_coin(
+                token_amount, i
+            )
+        except Exception:
+            if safe:
+                raise
+            return
+
+        measured = coins[i].balanceOf(user)
+        d_balances = [swap_with_deposit.balances(k) for k in range(len(coins))]
+        claimed_fees = [0] * len(coins)
+        try:
+            with boa.env.prank(user):
+                swap_with_deposit.remove_liquidity_one_coin(
+                    token_amount, i, int(0.999 * calculated)
+                )
+
+            logs = swap_with_deposit.get_logs()
+            for log in logs:
+                if log.event_type.name == "ClaimAdminFee":
+                    claimed_fees = log.args[0]
+
+        except Exception:
+
+            # Check if it could fall into unsafe region here
+            frac = (
+                (d_balances[i] - calculated)
+                * (INITIAL_PRICES)[i]
+                // swap_with_deposit.D()
+            )
+
+            if frac > 1.1e16 and frac < 0.9e20:
+                raise
+            return  # dont continue tests
+
+        d_balances = [
+            d_balances[k] - swap_with_deposit.balances(k)
+            for k in range(len(coins))
+        ]
+        measured = coins[i].balanceOf(user) - measured
+
+        assert calculated >= measured
+        assert calculated - (0.01 / 100) * calculated < measured
+        assert approx(calculated, measured, 1e-3)
+
+        for k in range(len(coins)):
+            claimed_tokens = claimed_fees[k]
+            if k == i:
+                assert d_balances[k] == measured + claimed_tokens
+            else:
+                assert d_balances[k] == claimed_tokens
+
+        # This is to check that we didn't end up in a borked state after
+        # a withdrawal succeeded
+        views_contract.get_dy(0, 1, 10**16, swap_with_deposit)
diff --git a/tests/unitary/pool/test_exchange.py b/tests/unitary/pool/test_exchange.py
new file mode 100644
index 00000000..058589a3
--- /dev/null
+++ b/tests/unitary/pool/test_exchange.py
@@ -0,0 +1,150 @@
+import boa
+from boa.test import strategy
+from hypothesis import given, settings  # noqa
+
+from tests.fixtures.pool import INITIAL_PRICES
+from tests.utils.tokens import mint_for_testing
+
+SETTINGS = {"max_examples": 100, "deadline": None}
+
+
+def test_exchange_reverts(user, views_contract, swap_with_deposit):
+
+    with boa.reverts():
+        views_contract.get_dy(0, 2, 10**6, swap_with_deposit)
+
+    with boa.reverts():
+        views_contract.get_dy(2, 1, 10**6, swap_with_deposit)
+
+    with boa.reverts():
+        swap_with_deposit.exchange(1, 3, 10**6, 0, sender=user)
+
+    with boa.reverts():
+        swap_with_deposit.exchange(2, 2, 10**6, 0, sender=user)
+
+
+@given(
+    amount=strategy(
+        "uint256", min_value=10**10, max_value=10**6 * 10**18
+    ),  # Can be more than we have
+    i=strategy("uint", min_value=0, max_value=1),
+    j=strategy("uint", min_value=0, max_value=1),
+)
+@settings(**SETTINGS)
+def test_exchange_all(
+    swap_with_deposit,
+    views_contract,
+    coins,
+    user,
+    amount,
+    i,
+    j,
+):
+
+    if i == j:
+        return
+
+    amount = amount * 10**18 // INITIAL_PRICES[i]
+    mint_for_testing(coins[i], user, amount)
+
+    calculated = views_contract.get_dy(i, j, amount, swap_with_deposit)
+
+    measured_i = coins[i].balanceOf(user)
+    measured_j = coins[j].balanceOf(user)
+    d_balance_i = swap_with_deposit.balances(i)
+    d_balance_j = swap_with_deposit.balances(j)
+
+    with boa.env.prank(user):
+        swap_with_deposit.exchange(i, j, amount, int(0.999 * calculated))
+
+    measured_i -= coins[i].balanceOf(user)
+    measured_j = coins[j].balanceOf(user) - measured_j
+    d_balance_i = swap_with_deposit.balances(i) - d_balance_i
+    d_balance_j = swap_with_deposit.balances(j) - d_balance_j
+
+    assert amount == measured_i
+    assert calculated == measured_j
+
+    assert d_balance_i == amount
+    assert -d_balance_j == measured_j
+
+
+@given(
+    amount=strategy(
+        "uint256", min_value=10**10, max_value=10**6 * 10**18
+    ),  # Can be more than we have
+    i=strategy("uint", min_value=0, max_value=1),
+    j=strategy("uint", min_value=0, max_value=1),
+)
+@settings(**SETTINGS)
+def test_exchange_received_success(
+    swap_with_deposit,
+    views_contract,
+    coins,
+    user,
+    amount,
+    i,
+    j,
+):
+
+    if i == j:
+        return
+
+    amount = amount * 10**18 // INITIAL_PRICES[i]
+    mint_for_testing(coins[i], user, amount)
+
+    calculated = views_contract.get_dy(i, j, amount, swap_with_deposit)
+
+    measured_i = coins[i].balanceOf(user)
+    measured_j = coins[j].balanceOf(user)
+    d_balance_i = swap_with_deposit.balances(i)
+    d_balance_j = swap_with_deposit.balances(j)
+
+    with boa.env.prank(user):
+        coins[i].transfer(swap_with_deposit, amount)
+        out = swap_with_deposit.exchange_received(
+            i, j, amount, int(0.999 * calculated), user
+        )
+
+    measured_i -= coins[i].balanceOf(user)
+    measured_j = coins[j].balanceOf(user) - measured_j
+    d_balance_i = swap_with_deposit.balances(i) - d_balance_i
+    d_balance_j = swap_with_deposit.balances(j) - d_balance_j
+
+    assert amount == measured_i
+    assert calculated == measured_j == out
+
+    assert d_balance_i == amount
+    assert -d_balance_j == measured_j == out
+
+
+@given(
+    amount=strategy(
+        "uint256", min_value=10**10, max_value=10**6 * 10**18
+    ),  # Can be more than we have
+    i=strategy("uint", min_value=0, max_value=1),
+    j=strategy("uint", min_value=0, max_value=1),
+)
+@settings(**SETTINGS)
+def test_exchange_received_revert_on_no_transfer(
+    swap_with_deposit,
+    views_contract,
+    coins,
+    user,
+    amount,
+    i,
+    j,
+):
+
+    if i == j:
+        return
+
+    amount = amount * 10**18 // INITIAL_PRICES[i]
+    mint_for_testing(coins[i], user, amount)
+
+    calculated = views_contract.get_dy(i, j, amount, swap_with_deposit)
+
+    with boa.env.prank(user), boa.reverts(dev="user didn't give us coins"):
+        swap_with_deposit.exchange_received(
+            i, j, amount, int(0.999 * calculated), user
+        )
diff --git a/tests/unitary/pool/test_oracles.py b/tests/unitary/pool/test_oracles.py
new file mode 100644
index 00000000..80666a0c
--- /dev/null
+++ b/tests/unitary/pool/test_oracles.py
@@ -0,0 +1,219 @@
+from math import exp, log2, sqrt
+
+import boa
+import pytest
+from boa.test import strategy
+from hypothesis import given, settings
+
+from tests.fixtures.pool import INITIAL_PRICES
+from tests.utils import approx
+from tests.utils.tokens import mint_for_testing
+
+SETTINGS = {"max_examples": 1000, "deadline": None}
+
+
+def norm(price_oracle, price_scale):
+    norm = 0
+    ratio = price_oracle * 10**18 / price_scale
+    if ratio > 10**18:
+        ratio -= 10**18
+    else:
+        ratio = 10**18 - ratio
+    norm += ratio**2
+    return sqrt(norm)
+
+
+def test_initial(swap_with_deposit):
+    assert swap_with_deposit.price_scale() == INITIAL_PRICES[1]
+    assert swap_with_deposit.price_oracle() == INITIAL_PRICES[1]
+
+
+@given(
+    token_frac=strategy("uint256", min_value=10**6, max_value=10**16),
+    i=strategy("uint8", max_value=1),
+)
+@settings(**SETTINGS)
+def test_last_price_remove_liq(swap_with_deposit, user, token_frac, i):
+
+    prices = INITIAL_PRICES
+    token_amount = token_frac * swap_with_deposit.totalSupply() // 10**18
+
+    with boa.env.prank(user):
+        swap_with_deposit.remove_liquidity_one_coin(token_amount, i, 0)
+
+    oracle_price = swap_with_deposit.last_prices()
+    assert abs(log2(oracle_price / prices[1])) < 0.1
+
+
+@given(
+    amount=strategy(
+        "uint256", min_value=10**10, max_value=2 * 10**6 * 10**18
+    ),  # Can be more than we have
+    i=strategy("uint8", min_value=0, max_value=1),
+    t=strategy("uint256", min_value=10, max_value=10 * 86400),
+)
+@settings(**SETTINGS)
+def test_ma(swap_with_deposit, coins, user, amount, i, t):
+
+    prices1 = INITIAL_PRICES
+    amount = amount * 10**18 // prices1[i]
+    mint_for_testing(coins[i], user, amount)
+
+    rebal_params = swap_with_deposit.internal._unpack_3(
+        swap_with_deposit._storage.packed_rebalancing_params.get()
+    )
+    ma_time = rebal_params[2]
+
+    # here we dont mine because we're time travelling later
+    with boa.env.prank(user):
+        swap_with_deposit.exchange(i, 1 - i, amount, 0)
+
+    prices2 = swap_with_deposit.last_prices()
+    boa.env.time_travel(t)
+
+    with boa.env.prank(user):
+        swap_with_deposit.remove_liquidity_one_coin(10**15, 0, 0)
+
+    prices3 = swap_with_deposit.price_oracle()
+
+    # cap new price by 2x previous price oracle value:
+    new_price = min(prices2, 2 * prices1[1])
+
+    alpha = exp(-1 * t / ma_time)
+    theory = prices1[1] * alpha + new_price * (1 - alpha)
+    assert abs(log2(theory / prices3)) < 0.001
+
+
+@given(
+    amount=strategy(
+        "uint256", min_value=10**10, max_value=2 * 10**6 * 10**18
+    ),  # Can be more than we have
+    i=strategy("uint8", min_value=0, max_value=1),
+    t=strategy("uint256", min_value=10, max_value=10 * 86400),
+)
+@settings(**SETTINGS)
+def test_xcp_ma(swap_with_deposit, coins, user, amount, i, t):
+
+    price_scale = swap_with_deposit.price_scale()
+    D0 = swap_with_deposit.D()
+    xp = [0, 0]
+    xp[0] = D0 // 2  # N_COINS = 2
+    xp[1] = D0 * 10**18 // (2 * price_scale)
+
+    xcp0 = boa.eval(f"isqrt({xp[0]*xp[1]})")
+
+    # after first deposit anf before any swaps:
+    # xcp oracle is equal to totalSupply
+    assert xcp0 == swap_with_deposit.totalSupply()
+
+    amount = amount * 10**18 // INITIAL_PRICES[i]
+    mint_for_testing(coins[i], user, amount)
+
+    ma_time = swap_with_deposit.xcp_ma_time()
+
+    # swap to populate
+    with boa.env.prank(user):
+        swap_with_deposit.exchange(i, 1 - i, amount, 0)
+
+    xcp1 = swap_with_deposit.last_xcp()
+    tvl = (
+        swap_with_deposit.virtual_price()
+        * swap_with_deposit.totalSupply()
+        // 10**18
+    )
+    assert approx(xcp1, tvl, 1e-10)
+
+    boa.env.time_travel(t)
+
+    with boa.env.prank(user):
+        swap_with_deposit.remove_liquidity_one_coin(10**15, 0, 0)
+
+    xcp2 = swap_with_deposit.xcp_oracle()
+
+    alpha = exp(-1 * t / ma_time)
+    theory = xcp0 * alpha + xcp1 * (1 - alpha)
+
+    assert approx(theory, xcp2, 1e-10)
+
+
+# Sanity check for price scale
+@given(
+    amount=strategy(
+        "uint256", min_value=10**10, max_value=2 * 10**6 * 10**18
+    ),  # Can be more than we have
+    i=strategy("uint8", min_value=0, max_value=1),
+    t=strategy("uint256", max_value=10 * 86400),
+)
+@settings(**SETTINGS)
+def test_price_scale_range(swap_with_deposit, coins, user, amount, i, t):
+
+    prices1 = INITIAL_PRICES
+    amount = amount * 10**18 // prices1[i]
+    mint_for_testing(coins[i], user, amount)
+
+    with boa.env.prank(user):
+        swap_with_deposit.exchange(i, 1 - i, amount, 0)
+
+    prices2 = swap_with_deposit.last_prices()
+    boa.env.time_travel(seconds=t)
+
+    with boa.env.prank(user):
+        swap_with_deposit.remove_liquidity_one_coin(10**15, 0, 0)
+
+    prices3 = swap_with_deposit.price_scale()
+
+    if prices1[1] > prices2:
+        assert prices3 <= prices1[1] and prices3 >= prices2
+    else:
+        assert prices3 >= prices1[1] and prices3 <= prices2
+
+
+@pytest.mark.parametrize("i", [0, 1])
+def test_price_scale_change(swap_with_deposit, i, coins, users):
+    j = 1 - i
+    amount = 10**6 * 10**18
+    t = 86400
+    user = users[1]
+    prices1 = INITIAL_PRICES
+    amount = amount * 10**18 // prices1[i]
+    mint_for_testing(coins[i], user, amount)
+    mint_for_testing(coins[j], user, amount)
+    coins[i].approve(swap_with_deposit, 2**256 - 1, sender=user)
+    coins[j].approve(swap_with_deposit, 2**256 - 1, sender=user)
+
+    out = swap_with_deposit.exchange(i, j, amount, 0, sender=user)
+    swap_with_deposit.exchange(j, i, int(out * 0.95), 0, sender=user)
+    price_scale_1 = swap_with_deposit.price_scale()
+
+    boa.env.time_travel(t)
+
+    swap_with_deposit.exchange(0, 1, coins[0].balanceOf(user), 0, sender=user)
+
+    price_oracle = swap_with_deposit.price_oracle()
+    rebal_params = swap_with_deposit.internal._unpack_3(
+        swap_with_deposit._storage.packed_rebalancing_params.get()
+    )
+    _norm = norm(price_oracle, price_scale_1)
+    step = max(rebal_params[1], _norm / 5)
+    price_scale_2 = swap_with_deposit.price_scale()
+
+    price_diff = abs(price_scale_2 - price_scale_1)
+    adjustment = int(step * abs(price_oracle - price_scale_1) / _norm)
+    assert price_diff > 0
+    assert approx(adjustment, price_diff, 0.01)
+    assert approx(
+        swap_with_deposit.virtual_price(),
+        swap_with_deposit.get_virtual_price(),
+        1e-10,
+    )
+
+
+def test_lp_price(swap_with_deposit):
+    tvl = (
+        swap_with_deposit.balances(0)
+        + swap_with_deposit.balances(1)
+        * swap_with_deposit.price_scale()
+        // 10**18
+    )
+    naive_price = tvl * 10**18 // swap_with_deposit.totalSupply()
+    assert approx(naive_price, swap_with_deposit.lp_price(), 0)
diff --git a/tests/unitary/pool/token/conftest.py b/tests/unitary/pool/token/conftest.py
new file mode 100644
index 00000000..f4149cb2
--- /dev/null
+++ b/tests/unitary/pool/token/conftest.py
@@ -0,0 +1,58 @@
+from copy import deepcopy
+
+import boa
+import pytest
+from eth_account import Account as EthAccount
+from eth_account._utils.structured_data.hashing import (
+    hash_domain,
+    hash_message,
+)
+from eth_account.messages import SignableMessage
+from hexbytes import HexBytes
+
+
+@pytest.fixture(scope="module")
+def sign_permit():
+    def _sign_permit(swap, owner, spender, value, deadline):
+
+        PERMIT_STRUCT = {
+            "types": {
+                "EIP712Domain": [
+                    {"name": "name", "type": "string"},
+                    {"name": "version", "type": "string"},
+                    {"name": "chainId", "type": "uint256"},
+                    {"name": "verifyingContract", "type": "address"},
+                    {"name": "salt", "type": "bytes32"},
+                ],
+                "Permit": [
+                    {"name": "owner", "type": "address"},
+                    {"name": "spender", "type": "address"},
+                    {"name": "value", "type": "uint256"},
+                    {"name": "nonce", "type": "uint256"},
+                    {"name": "deadline", "type": "uint256"},
+                ],
+            },
+            "primaryType": "Permit",
+        }
+
+        struct = deepcopy(PERMIT_STRUCT)
+        struct["domain"] = dict(
+            name=swap.name(),
+            version=swap.version(),
+            chainId=boa.env.vm.chain_context.chain_id,
+            verifyingContract=swap.address,
+            salt=HexBytes(swap.salt()),
+        )
+        struct["message"] = dict(
+            owner=owner.address,
+            spender=spender,
+            value=value,
+            nonce=swap.nonces(owner.address),
+            deadline=deadline,
+        )
+        signable_message = SignableMessage(
+            b"\x01", hash_domain(struct), hash_message(struct)
+        )
+        return EthAccount.sign_message(signable_message, owner._private_key)
+
+    return _sign_permit
diff --git a/tests/unitary/pool/token/test_approve.py b/tests/unitary/pool/token/test_approve.py
new file mode 100644
index 00000000..fd4e8546
--- /dev/null
+++ b/tests/unitary/pool/token/test_approve.py
@@ -0,0 +1,81 @@
+import boa
+import pytest
+
+from tests.utils.tokens import mint_for_testing
+
+
+@pytest.mark.parametrize("idx", range(5))
+def test_initial_approval_is_zero(swap, alice, users, idx):
+    assert swap.allowance(alice, users[idx]) == 0
+
+
+def test_approve(swap, alice, bob):
+
+    with boa.env.prank(alice):
+        swap.approve(bob, 10**19)
+
+    assert swap.allowance(alice, bob) == 10**19
+
+
+def test_modify_approve_zero_nonzero(swap, alice, bob):
+
+    with boa.env.prank(alice):
+        swap.approve(bob, 10**19)
+        swap.approve(bob, 0)
+        swap.approve(bob, 12345678)
+
+    assert swap.allowance(alice, bob) == 12345678
+
+
+def test_revoke_approve(swap, alice, bob):
+
+    with boa.env.prank(alice):
+        swap.approve(bob, 10**19)
+        swap.approve(bob, 0)
+
+    assert swap.allowance(alice, bob) == 0
+
+
+def test_approve_self(swap, alice):
+
+    with boa.env.prank(alice):
+        swap.approve(alice, 10**19)
+
+    assert swap.allowance(alice, alice) == 10**19
+
+
+def test_only_affects_target(swap, alice, bob):
+    with boa.env.prank(alice):
+        swap.approve(bob, 10**19)
+
+    assert swap.allowance(bob, alice) == 0
+
+
+def test_returns_true(swap, alice, bob):
+    with boa.env.prank(alice):
+        assert swap.approve(bob, 10**19)
+
+
+def test_approval_event_fires(swap, alice, bob):
+
+    with boa.env.prank(alice):
+        swap.approve(bob, 10**19)
+
+    logs = swap.get_logs()
+
+    assert len(logs) == 1
+    assert logs[0].event_type.name == "Approval"
+    assert logs[0].topics[0].lower() == alice.lower()
+    assert logs[0].topics[1].lower() == bob.lower()
+    assert logs[0].args[0] == 10**19
+
+
+def test_infinite_approval(swap, alice, bob):
+    with boa.env.prank(alice):
+        swap.approve(bob, 2**256 - 1)
+
+    mint_for_testing(swap, alice, 10**18)
+    with boa.env.prank(bob):
+        swap.transferFrom(alice, bob, 10**18)
+
+    assert swap.allowance(alice, bob) == 2**256 - 1
diff --git a/tests/unitary/pool/token/test_permit.py b/tests/unitary/pool/token/test_permit.py
new file mode 100644
index 00000000..83493b21
--- /dev/null
+++ b/tests/unitary/pool/token/test_permit.py
@@ -0,0 +1,92 @@
+import boa
+from hexbytes import HexBytes
+
+ZERO_ADDRESS = "0x0000000000000000000000000000000000000000"
+
+
+# tests inspired by:
+# https://github.com/yearn/yearn-vaults/blob/master/tests/functional/vault/test_permit.py  # noqa: E501
+# https://github.com/curvefi/curve-stablecoin/blob/5b6708138d82419917328e8042f3857eac034796/tests/stablecoin/test_approve.py  # noqa: E501
+
+
+def test_permit_success(eth_acc, bob, swap, sign_permit):
+
+    value = 2**256 - 1
+    deadline = boa.env.vm.state.timestamp + 600
+
+    sig = sign_permit(
+        swap=swap,
+        owner=eth_acc,
+        spender=bob,
+        value=value,
+        deadline=deadline,
+    )
+
+    assert swap.allowance(eth_acc.address, bob) == 0
+    with boa.env.prank(bob):
+        assert swap.permit(
+            eth_acc.address,
+            bob,
+            value,
+            deadline,
+            sig.v,
+            HexBytes(sig.r),
+            HexBytes(sig.s),
+        )
+
+    logs = swap.get_logs()
+
+    assert swap.allowance(eth_acc.address, bob) == value
+    assert len(logs) == 1
+    assert swap.nonces(eth_acc.address) == 1
+    assert logs[0].event_type.name == "Approval"
+    assert logs[0].topics[0].lower() == eth_acc.address.lower()
+    assert logs[0].topics[1].lower() == bob.lower()
+    assert logs[0].args[0] == value
+
+
+def test_permit_reverts_owner_is_invalid(bob, swap):
+    with boa.reverts(dev="invalid owner"), boa.env.prank(bob):
+        swap.permit(
+            ZERO_ADDRESS,
+            bob,
+            2**256 - 1,
+            boa.env.vm.state.timestamp + 600,
+            27,
+            b"\x00" * 32,
+            b"\x00" * 32,
+        )
+
+
+def test_permit_reverts_deadline_is_invalid(bob, swap):
+    with boa.reverts(dev="permit expired"), boa.env.prank(bob):
+        swap.permit(
+            bob,
+            bob,
+            2**256 - 1,
+            boa.env.vm.state.timestamp - 600,
+            27,
+            b"\x00" * 32,
+            b"\x00" * 32,
+        )
+
+
+def test_permit_reverts_signature_is_invalid(bob, swap):
+    with boa.reverts(dev="invalid signature"), boa.env.prank(bob):
+        swap.permit(
+            bob,
+            bob,
+            2**256 - 1,
+            boa.env.vm.state.timestamp + 600,
+            27,
+            b"\x00" * 32,
+            b"\x00" * 32,
+        )
+
+
+def test_domain_separator_updates_when_chain_id_updates(swap):
+
+    domain_separator = swap.DOMAIN_SEPARATOR()
+    with boa.env.anchor():
+        boa.env.vm.patch.chain_id = 42
+        assert domain_separator != swap.DOMAIN_SEPARATOR()
diff --git a/tests/unitary/pool/token/test_transfer.py b/tests/unitary/pool/token/test_transfer.py
new file mode 100644
index 00000000..4cbdb48f
--- /dev/null
+++ b/tests/unitary/pool/token/test_transfer.py
@@ -0,0 +1,91 @@
+import boa
+
+
+def test_sender_balance_decreases(loaded_alice, bob, swap):
+    sender_balance = swap.balanceOf(loaded_alice)
+    amount = sender_balance // 4
+
+    with boa.env.prank(loaded_alice):
+        swap.transfer(bob, amount)
+
+    assert swap.balanceOf(loaded_alice) == sender_balance - amount
+
+
+def test_receiver_balance_increases(loaded_alice, bob, swap):
+    receiver_balance = swap.balanceOf(bob)
+    amount = swap.balanceOf(loaded_alice) // 4
+
+    with boa.env.prank(loaded_alice):
+        swap.transfer(bob, amount)
+
+    assert swap.balanceOf(bob) == receiver_balance + amount
+
+
+def test_total_supply_not_affected(loaded_alice, bob, swap):
+    total_supply = swap.totalSupply()
+    amount = swap.balanceOf(loaded_alice)
+
+    with boa.env.prank(loaded_alice):
+        swap.transfer(bob, amount)
+
+    assert swap.totalSupply() == total_supply
+
+
+def test_returns_true(loaded_alice, bob, swap):
+    amount = swap.balanceOf(loaded_alice)
+
+    with boa.env.prank(loaded_alice):
+        assert swap.transfer(bob, amount)
+
+
+def test_transfer_full_balance(loaded_alice, bob, swap):
+    amount = swap.balanceOf(loaded_alice)
+    receiver_balance = swap.balanceOf(bob)
+
+    with boa.env.prank(loaded_alice):
+        swap.transfer(bob, amount)
+
+    assert swap.balanceOf(loaded_alice) == 0
+    assert swap.balanceOf(bob) == receiver_balance + amount
+
+
+def test_transfer_zero_tokens(loaded_alice, bob, swap):
+    sender_balance = swap.balanceOf(loaded_alice)
+    receiver_balance = swap.balanceOf(bob)
+
+    with boa.env.prank(loaded_alice):
+        swap.transfer(bob, 0)
+
+    assert swap.balanceOf(loaded_alice) == sender_balance
+    assert swap.balanceOf(bob) == receiver_balance
+
+
+def test_transfer_to_self(loaded_alice, swap):
+    sender_balance = swap.balanceOf(loaded_alice)
+    amount = sender_balance // 4
+
+    with boa.env.prank(loaded_alice):
+        swap.transfer(loaded_alice, amount)
+
+    assert swap.balanceOf(loaded_alice) == sender_balance
+
+
+def test_insufficient_balance(loaded_alice, bob, swap):
+    balance = swap.balanceOf(loaded_alice)
+
+    with boa.reverts(), boa.env.prank(loaded_alice):
+        swap.transfer(bob, balance + 1)
+
+
+def test_transfer_event_fires(loaded_alice, bob, swap):
+    amount = swap.balanceOf(loaded_alice)
+    with boa.env.prank(loaded_alice):
+        swap.transfer(bob, amount)
+
+    logs = swap.get_logs()
+
+    assert len(logs) == 1
+    assert logs[0].event_type.name == "Transfer"
+    assert logs[0].args[0] == amount
+    assert logs[0].topics[0].lower() == loaded_alice.lower()
+    assert logs[0].topics[1].lower() == bob.lower()
diff --git a/tests/unitary/pool/token/test_transferFrom.py b/tests/unitary/pool/token/test_transferFrom.py
new file mode 100644
index 00000000..29060699
--- /dev/null
+++ b/tests/unitary/pool/token/test_transferFrom.py
@@ -0,0 +1,212 @@
+import boa
+
+
+def test_sender_balance_decreases(loaded_alice, bob, charlie, swap):
+    sender_balance = swap.balanceOf(loaded_alice)
+    amount = sender_balance // 4
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, amount)
+
+    with boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, amount)
+
+    assert swap.balanceOf(loaded_alice) == sender_balance - amount
+
+
+def test_receiver_balance_increases(loaded_alice, bob, charlie, swap):
+    receiver_balance = swap.balanceOf(charlie)
+    amount = swap.balanceOf(loaded_alice) // 4
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, amount)
+
+    with boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, amount)
+
+    assert swap.balanceOf(charlie) == receiver_balance + amount
+
+
+def test_caller_balance_not_affected(loaded_alice, bob, charlie, swap):
+    caller_balance = swap.balanceOf(bob)
+    amount = swap.balanceOf(loaded_alice)
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, amount)
+
+    with boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, amount)
+
+    assert swap.balanceOf(bob) == caller_balance
+
+
+def test_caller_approval_affected(alice, bob, charlie, swap):
+    approval_amount = swap.balanceOf(alice)
+    transfer_amount = approval_amount // 4
+
+    with boa.env.prank(alice):
+        swap.approve(bob, approval_amount)
+
+    with boa.env.prank(bob):
+        swap.transferFrom(alice, charlie, transfer_amount)
+
+    assert swap.allowance(alice, bob) == approval_amount - transfer_amount
+
+
+def test_receiver_approval_not_affected(loaded_alice, bob, charlie, swap):
+    approval_amount = swap.balanceOf(loaded_alice)
+    transfer_amount = approval_amount // 4
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, approval_amount)
+        swap.approve(charlie, approval_amount)
+
+    with boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, transfer_amount)
+
+    assert swap.allowance(loaded_alice, charlie) == approval_amount
+
+
+def test_total_supply_not_affected(loaded_alice, bob, charlie, swap):
+    total_supply = swap.totalSupply()
+    amount = swap.balanceOf(loaded_alice)
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, amount)
+
+    with boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, amount)
+
+    assert swap.totalSupply() == total_supply
+
+
+def test_returns_true(loaded_alice, bob, charlie, swap):
+    amount = swap.balanceOf(loaded_alice)
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, amount)
+
+    with boa.env.prank(bob):
+        assert swap.transferFrom(loaded_alice, charlie, amount)
+
+
+def test_transfer_full_balance(loaded_alice, bob, charlie, swap):
+    amount = swap.balanceOf(loaded_alice)
+    receiver_balance = swap.balanceOf(charlie)
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, amount)
+
+    with boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, amount)
+
+    assert swap.balanceOf(loaded_alice) == 0
+    assert swap.balanceOf(charlie) == receiver_balance + amount
+
+
+def test_transfer_zero_tokens(loaded_alice, bob, charlie, swap):
+    sender_balance = swap.balanceOf(loaded_alice)
+    receiver_balance = swap.balanceOf(charlie)
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, sender_balance)
+
+    with boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, 0)
+
+    assert swap.balanceOf(loaded_alice) == sender_balance
+    assert swap.balanceOf(charlie) == receiver_balance
+
+
+def test_transfer_zero_tokens_without_approval(
+    loaded_alice, bob, charlie, swap
+):
+    sender_balance = swap.balanceOf(loaded_alice)
+    receiver_balance = swap.balanceOf(charlie)
+
+    with boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, 0)
+
+    assert swap.balanceOf(loaded_alice) == sender_balance
+    assert swap.balanceOf(charlie) == receiver_balance
+
+
+def test_insufficient_balance(loaded_alice, bob, charlie, swap):
+    balance = swap.balanceOf(loaded_alice)
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, balance + 1)
+
+    with boa.reverts(), boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, balance + 1)
+
+
+def test_insufficient_approval(loaded_alice, bob, charlie, swap):
+    balance = swap.balanceOf(loaded_alice)
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, balance - 1)
+
+    with boa.reverts(), boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, balance)
+
+
+def test_no_approval(loaded_alice, bob, charlie, swap):
+    balance = swap.balanceOf(loaded_alice)
+
+    with boa.reverts(), boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, balance)
+
+
+def test_revoked_approval(loaded_alice, bob, charlie, swap):
+    balance = swap.balanceOf(loaded_alice)
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, balance)
+        swap.approve(bob, 0)
+
+    with boa.reverts(), boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, balance)
+
+
+def test_transfer_to_self(loaded_alice, swap):
+    sender_balance = swap.balanceOf(loaded_alice)
+    amount = sender_balance // 4
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(loaded_alice, sender_balance)
+        swap.transferFrom(loaded_alice, loaded_alice, amount)
+
+    assert swap.balanceOf(loaded_alice) == sender_balance
+    assert (
+        swap.allowance(loaded_alice, loaded_alice) == sender_balance - amount
+    )
+
+
+def test_transfer_to_self_no_approval(loaded_alice, swap):
+    amount = swap.balanceOf(loaded_alice)
+
+    with boa.reverts(), boa.env.prank(loaded_alice):
+        swap.transferFrom(loaded_alice, loaded_alice, amount)
+
+
+def test_transfer_event_fires(loaded_alice, bob, charlie, swap):
+    amount = swap.balanceOf(loaded_alice)
+
+    with boa.env.prank(loaded_alice):
+        swap.approve(bob, amount)
+
+    with boa.env.prank(bob):
+        swap.transferFrom(loaded_alice, charlie, amount)
+
+    logs = swap.get_logs()
+
+    assert len(logs) == 2
+    assert logs[0].event_type.name == "Approval"
+    assert logs[0].args[0] == 0  # since everything got transferred
+    assert logs[0].topics[0].lower() == loaded_alice.lower()
+    assert logs[0].topics[1].lower() == bob.lower()
+
+    assert logs[1].event_type.name == "Transfer"
+    assert logs[1].args[0] == amount
+    assert logs[1].topics[0].lower() == loaded_alice.lower()
+    assert logs[1].topics[1].lower() == charlie.lower()
diff --git a/tests/utils/__init__.py b/tests/utils/__init__.py
index e69de29b..c10431db 100644
--- a/tests/utils/__init__.py
+++ b/tests/utils/__init__.py
@@ -0,0 +1,5 @@
+import math
+
+
+def approx(x1, x2, precision):
+    return abs(math.log(x1 / x2)) <= precision
diff --git a/tests/utils/simulation_int_many.py b/tests/utils/simulation_int_many.py
new file mode 100644
index 00000000..d7d4d9e8
--- /dev/null
+++ b/tests/utils/simulation_int_many.py
@@ -0,0 +1,353 @@
+#!/usr/bin/env python3
+# flake8: noqa
+from math import exp
+
+from tests.unitary.math.misc import get_y_n2_dec
+
+A_MULTIPLIER = 10000
+
+
+def geometric_mean(x):
+    N = len(x)
+    x = sorted(x, reverse=True)  # Presort - good for convergence
+    D = x[0]
+    for i in range(255):
+        D_prev = D
+        tmp = 10**18
+        for _x in x:
+            tmp = tmp * _x // D
+        D = D * ((N - 1) * 10**18 + tmp) // (N * 10**18)
+        diff = abs(D - D_prev)
+        if diff <= 1 or diff * 10**18 < D:
+            return D
+    raise ValueError("Did not converge")
+
+
+def reduction_coefficient(x, gamma):
+    N = len(x)
+    x_prod = 10**18
+    K = 10**18
+    S = sum(x)
+    for x_i in x:
+        x_prod = x_prod * x_i // 10**18
+        K = K * N * x_i // S
+    if gamma > 0:
+        K = gamma * 10**18 // (gamma + 10**18 - K)
+    return K
+
+
+def get_fee(x, fee_gamma, mid_fee, out_fee):
+    f = reduction_coefficient(x, fee_gamma)
+    return (mid_fee * f + out_fee * (10**18 - f)) // 10**18
+
+
+def newton_D(A, gamma, x, D0):
+    D = D0
+    i = 0
+
+    S = sum(x)
+    x = sorted(x, reverse=True)
+    N = len(x)
+
+    assert N == 2
+
+    for i in range(255):
+        D_prev = D
+
+        K0 = 10**18
+        for _x in x:
+            K0 = K0 * _x * N // D
+
+        _g1k0 = abs(gamma + 10**18 - K0)
+
+        # D / (A * N**N) * _g1k0**2 / gamma**2
+        mul1 = (
+            10**18 * D // gamma * _g1k0 // gamma * _g1k0 * A_MULTIPLIER // A
+        )
+
+        # 2*N*K0 / _g1k0
+        mul2 = (2 * 10**18) * N * K0 // _g1k0
+
+        neg_fprime = (
+            (S + S * mul2 // 10**18) + mul1 * N // K0 - mul2 * D // 10**18
+        )
+        assert neg_fprime > 0  # Python only: -f' > 0
+
+        # D -= f / fprime
+        D = (D * neg_fprime + D * S - D**2) // neg_fprime - D * (
+            mul1 // neg_fprime
+        ) // 10**18 * (10**18 - K0) // K0
+
+        if D < 0:
+            D = -D // 2
+        if abs(D - D_prev) <= max(100, D // 10**14):
+            return D
+
+    raise ValueError("Did not converge")
+
+
+def newton_y(A, gamma, x, D, i):
+    N = len(x)
+
+    assert N == 2
+
+    y = D // N
+    K0_i = 10**18
+    S_i = 0
+    x_sorted = sorted(_x for j, _x in enumerate(x) if j != i)
+    convergence_limit = max(max(x_sorted) // 10**14, D // 10**14, 100)
+    for _x in x_sorted:
+        y = y * D // (_x * N)  # Small _x first
+        S_i += _x
+    for _x in x_sorted[::-1]:
+        K0_i = K0_i * _x * N // D  # Large _x first
+
+    for j in range(255):
+        y_prev = y
+
+        K0 = K0_i * y * N // D
+        S = S_i + y
+
+        _g1k0 = abs(gamma + 10**18 - K0)
+
+        # D / (A * N**N) * _g1k0**2 / gamma**2
+        mul1 = (
+            10**18 * D // gamma * _g1k0 // gamma * _g1k0 * A_MULTIPLIER // A
+        )
+
+        # 2*K0 / _g1k0
+        mul2 = 10**18 + (2 * 10**18) * K0 // _g1k0
+
+        yfprime = 10**18 * y + S * mul2 + mul1 - D * mul2
+        fprime = yfprime // y
+        assert fprime > 0  # Python only: f' > 0
+
+        # y -= f / f_prime;  y = (y * fprime - f) / fprime
+        y = (
+            yfprime + 10**18 * D - 10**18 * S
+        ) // fprime + mul1 // fprime * (10**18 - K0) // K0
+
+        if j > 100:  # Just logging when doesn't converge
+            print(j, y, D, x)
+        if y < 0 or fprime < 0:
+            y = y_prev // 2
+        if abs(y - y_prev) <= max(convergence_limit, y // 10**14):
+            return y
+
+    raise Exception("Did not converge")
+
+
+def solve_x(A, gamma, x, D, i):
+    return int(get_y_n2_dec(A, gamma, x, D, i)[0] * 10**18)
+    # return newton_y(A, gamma, x, D, i)
+
+
+def solve_D(A, gamma, x):
+    D0 = len(x) * geometric_mean(x)  # <- fuzz to make sure it's ok XXX
+    return newton_D(A, gamma, x, D0)
+
+
+class Curve:
+    def __init__(self, A, gamma, D, n, p):
+        self.A = A
+        self.gamma = gamma
+        self.n = n
+        self.p = p
+        self.x = [D // n * 10**18 // self.p[i] for i in range(n)]
+
+    def xp(self):
+        return [x * p // 10**18 for x, p in zip(self.x, self.p)]
+
+    def D(self):
+        xp = self.xp()
+        if any(x <= 0 for x in xp):
+            raise ValueError
+        return solve_D(self.A, self.gamma, xp)
+
+    def y(self, x, i, j):
+        xp = self.xp()
+        xp[i] = x * self.p[i] // 10**18
+        yp = solve_x(self.A, self.gamma, xp, self.D(), j)
+        return yp * 10**18 // self.p[j]
+
+    def get_p(self):
+
+        A = self.A
+        gamma = self.gamma
+        xp = self.xp()
+        D = self.D()
+
+        K0 = xp[0] * xp[1] * 4 // D * 10**36 // D
+        gK0 = (
+            2 * K0 * K0 // 10**36 * K0 // 10**36
+            + (gamma + 10**18) ** 2
+            - (K0 * K0 // 10**36 * (2 * gamma + 3 * 10**18) // 10**18)
+        )
+        NNAG2 = A * gamma**2 // A_MULTIPLIER
+        numerator = (
+            xp[0] * (gK0 + NNAG2 * xp[1] // D * K0 // 10**36) // xp[1]
+        )
+        denominator = gK0 + NNAG2 * xp[0] // D * K0 // 10**36
+
+        return numerator * 10**18 // denominator
+
+
+class Trader:
+    def __init__(
+        self,
+        A,
+        gamma,
+        D,
+        n,
+        p0,
+        mid_fee=1e-3,
+        out_fee=3e-3,
+        allowed_extra_profit=2 * 10**13,
+        fee_gamma=None,
+        adjustment_step=0.003,
+        ma_time=866,
+        log=True,
+    ):
+        # allowed_extra_profit is actually not used
+        self.p0 = p0[:]
+        self.price_oracle = self.p0[:]
+        self.last_price = self.p0[:]
+        self.curve = Curve(A, gamma, D, n, p=p0[:])
+        self.dx = int(D * 1e-8)
+        self.mid_fee = int(mid_fee * 1e10)
+        self.out_fee = int(out_fee * 1e10)
+        self.D0 = self.curve.D()
+        self.xcp_0 = self.get_xcp()
+        self.xcp_profit = 10**18
+        self.xcp_profit_real = 10**18
+        self.xcp = self.xcp_0
+        self.allowed_extra_profit = allowed_extra_profit
+        self.adjustment_step = int(10**18 * adjustment_step)
+        self.log = log
+        self.fee_gamma = fee_gamma or gamma
+        self.total_vol = 0.0
+        self.ma_time = ma_time
+        self.ext_fee = 0  # 0.03e-2
+        self.slippage = 0
+        self.slippage_count = 0
+
+    def fee(self):
+        f = reduction_coefficient(self.curve.xp(), self.fee_gamma)
+        return (self.mid_fee * f + self.out_fee * (10**18 - f)) // 10**18
+
+    def get_xcp(self):
+        # First calculate the ideal balance
+        # Then calculate, what the constant-product would be
+        D = self.curve.D()
+        N = len(self.curve.x)
+        X = [D * 10**18 // (N * p) for p in self.curve.p]
+
+        return geometric_mean(X)
+
+    def update_xcp(self, only_real=False):
+        xcp = self.get_xcp()
+        self.xcp_profit_real = self.xcp_profit_real * xcp // self.xcp
+        if not only_real:
+            self.xcp_profit = self.xcp_profit * xcp // self.xcp
+        self.xcp = xcp
+
+    def buy(self, dx, i, j, max_price=1e100):
+        """
+        Buy y for x
+        """
+        try:
+            x_old = self.curve.x[:]
+            x = self.curve.x[i] + dx
+            y = self.curve.y(x, i, j)
+            dy = self.curve.x[j] - y
+            self.curve.x[i] = x
+            self.curve.x[j] = y
+            fee = self.fee()
+            self.curve.x[j] += dy * fee // 10**10
+            dy = dy * (10**10 - fee) // 10**10
+            if dx * 10**18 // dy > max_price or dy < 0:
+                self.curve.x = x_old
+                return False
+            self.update_xcp()
+            return dy
+        except ValueError:
+            return False
+
+    def sell(self, dy, i, j, min_price=0):
+        """
+        Sell y for x
+        """
+        try:
+            x_old = self.curve.x[:]
+            y = self.curve.x[j] + dy
+            x = self.curve.y(y, j, i)
+            dx = self.curve.x[i] - x
+            self.curve.x[i] = x
+            self.curve.x[j] = y
+            fee = self.fee()
+            self.curve.x[i] += dx * fee // 10**10
+            dx = dx * (10**10 - fee) // 10**10
+            if dx * 10**18 // dy < min_price or dx < 0:
+                self.curve.x = x_old
+                return False
+            self.update_xcp()
+            return dx
+        except ValueError:
+            return False
+
+    def _ma_multiplier(self, t):
+        return int(10**18 * exp(-1 * (t - self.t) / self.ma_time))
+
+    def ma_recorder(self, t, price_vector):
+        # need to convert this to exp!
+        # XXX what if every block only has p_b being last
+        N = len(price_vector)
+        if t > self.t:
+            alpha = self._ma_multiplier(t)
+            last_price = min(price_vector[1], 2 * self.curve.p[1])
+            self.price_oracle[1] = (
+                int(
+                    last_price * (10**18 - alpha)
+                    + self.price_oracle[1] * alpha
+                )
+                // 10**18
+            )
+            self.t = t
+
+    def tweak_price(self, t):
+
+        self.ma_recorder(t, self.last_price)
+        self.last_price[1] = self.curve.get_p() * self.curve.p[1] // 10**18
+
+        # update price_scale:
+        norm = int(
+            sum(
+                (p_real * 10**18 // p_target - 10**18) ** 2
+                for p_real, p_target in zip(self.price_oracle, self.curve.p)
+            )
+            ** 0.5
+        )
+        adjustment_step = max(self.adjustment_step, norm // 5)
+        if norm <= adjustment_step:
+            # Already close to the target price
+            return norm
+
+        price_scale_adjustment = (
+            adjustment_step * (self.price_oracle[1] - self.curve.p[1]) // norm
+        )
+        p_new = [10**18, self.curve.p[1] + price_scale_adjustment]
+
+        old_p = self.curve.p[:]
+        old_profit = self.xcp_profit_real
+        old_xcp = self.xcp
+
+        self.curve.p = p_new
+        self.update_xcp(only_real=True)
+
+        if 2 * (self.xcp_profit_real - 10**18) <= self.xcp_profit - 10**18:
+            # If real profit is less than half of maximum - revert params back
+            self.curve.p = old_p
+            self.xcp_profit_real = old_profit
+            self.xcp = old_xcp
+
+        return norm