diff --git a/contracts/mocks/newton_y_large_gamma.vy b/contracts/mocks/newton_y_large_gamma.vy
new file mode 100644
index 00000000..13d55ae9
--- /dev/null
+++ b/contracts/mocks/newton_y_large_gamma.vy
@@ -0,0 +1,102 @@
+# Minimized version of the math contracts before the gamma value expansion.
+# Additionally to the final value it also returns the number of iterations it took to find the value.
+# For testing purposes only.
+# From commit: 6dec22f6956cc04fb865d93c1e521f146e066cab
+
+N_COINS: constant(uint256) = 2
+A_MULTIPLIER: constant(uint256) = 10000
+
+MIN_GAMMA: constant(uint256) = 10**10
+MAX_GAMMA_SMALL: constant(uint256) = 2 * 10**16
+MAX_GAMMA: constant(uint256) = 3 * 10**17
+
+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, lim_mul: uint256) -> (uint256, uint256):
+    """
+    Calculating x[i] given other balances x[0..N_COINS-1] and invariant D
+    ANN = A * N**N
+    This is computationally expensive.
+    """
+
+    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
+
+    assert (K0_i >= unsafe_div(10**36, lim_mul)) and (K0_i <= lim_mul)  # dev: unsafe values x[i]
+
+    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):
+            return y, j
+
+    raise "Did not converge"
+
+
+@external
+@pure
+def newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i: uint256) -> (uint256, uint256):
+
+    # 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
+    lim_mul: uint256 = 100 * 10**18  # 100.0
+    if gamma > MAX_GAMMA_SMALL:
+        lim_mul = unsafe_div(unsafe_mul(lim_mul, MAX_GAMMA_SMALL), gamma)  # smaller than 100.0
+
+    iterations: uint256 = 0
+    y: uint256 = 0
+
+    y, iterations = self._newton_y(ANN, gamma, x, D, i, lim_mul)
+    frac: uint256 = y * 10**18 / D
+    assert (frac >= unsafe_div(10**36 / N_COINS, lim_mul)) and (frac <= unsafe_div(lim_mul, N_COINS))  # dev: unsafe value for y
+
+    return y, iterations
diff --git a/contracts/mocks/newton_y_small_gamma.vy b/contracts/mocks/newton_y_small_gamma.vy
new file mode 100644
index 00000000..a1c5a2b5
--- /dev/null
+++ b/contracts/mocks/newton_y_small_gamma.vy
@@ -0,0 +1,95 @@
+# Minimized version of the math contracts before the gamma value expansion.
+# Additionally to the final value it also returns the number of iterations it took to find the value.
+# For testing purposes only.
+# From commit: 1c800bd7937f63a9c278a220af846d322f356dd5
+
+N_COINS: constant(uint256) = 2
+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
+    This is computationally expensive.
+    """
+
+    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
+
+    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)
+
+    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):
+            return y
+
+    raise "Did not converge"
+
+
+@external
+@pure
+def newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i: uint256) -> uint256:
+
+    # 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
diff --git a/tests/unitary/math/test_newton_y.py b/tests/unitary/math/test_newton_y.py
new file mode 100644
index 00000000..e587c3ea
--- /dev/null
+++ b/tests/unitary/math/test_newton_y.py
@@ -0,0 +1,53 @@
+import boa
+import pytest
+from hypothesis import given, settings
+from hypothesis import strategies as st
+
+N_COINS = 2
+# MAX_SAMPLES = 1000000  # Increase for fuzzing
+MAX_SAMPLES = 10000
+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 = 3 * 10**17
+
+
+@pytest.fixture(scope="module")
+def math_large_gamma():
+    return boa.load("contracts/mocks/newton_y_large_gamma.vy")
+
+
+@pytest.fixture(scope="module")
+def math_small_gamma():
+    return boa.load("contracts/mock/newton_y_small_gamma.vy")
+
+
+@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=10**17 // 2, max_value=10**19 // 2
+    ),  # <- ratio 1e18 * x/D, typically 1e18 * 1
+    yD=st.integers(
+        min_value=10**17 // 2, max_value=10**19 // 2
+    ),  # <- 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_iteration_diff(math_large_gamma, A, D, xD, yD, gamma, j):
+    pass
+    # TODO: make a test that:
+    # - measures how many iterations it takes for the
+    #   old value to converge between the two versions
+    # - makes sure that we're converging to the correct value
+    # - use hypothesis.note to have some clear statistics about
+    #   the differences in divergence
+    # X = [D * xD // 10**18, D * yD // 10**18]
+    # math_large_gamma.newton_y(A, gamma, X, D, j)