From adfed8dec631d9dec1fe208ec532eba96cab14dd Mon Sep 17 00:00:00 2001 From: bout3fiddy <11488427+bout3fiddy@users.noreply.github.com> Date: Fri, 24 Nov 2023 15:12:45 +0100 Subject: [PATCH] add more newton_D checks --- contracts/main/CurveCryptoMathOptimized2.vy | 6 +- contracts/old/CurveCryptoSwap2Math.vy | 193 ++++++++++++++++++++ tests/unitary/math/conftest.py | 2 +- tests/unitary/math/test_newton_D_ref.py | 190 +++++++++++++++++++ tests/utils/simulation_int_many.py | 5 + 5 files changed, 392 insertions(+), 4 deletions(-) create mode 100644 contracts/old/CurveCryptoSwap2Math.vy create mode 100644 tests/unitary/math/test_newton_D_ref.py diff --git a/contracts/main/CurveCryptoMathOptimized2.vy b/contracts/main/CurveCryptoMathOptimized2.vy index 29c1ad20..9b155388 100644 --- a/contracts/main/CurveCryptoMathOptimized2.vy +++ b/contracts/main/CurveCryptoMathOptimized2.vy @@ -214,7 +214,7 @@ def newton_y(ANN: uint256, gamma: uint256, x: uint256[N_COINS], D: uint256, i: u 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 + assert (frac >= 10**16 - 1) and (frac < 10**20 + 1) # dev: unsafe value for y return y @@ -357,7 +357,7 @@ def get_y( 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 + assert (frac >= 10**16 - 1) and (frac < 10**20 + 1) # dev: unsafe value for y return y_out @@ -445,7 +445,7 @@ def newton_D(ANN: uint256, gamma: uint256, x_unsorted: uint256[N_COINS], K0_prev 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] + assert (frac >= 10**16 - 1) and (frac < 10**20 + 1) # dev: unsafe values x[i] return D raise "Did not converge" 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/tests/unitary/math/conftest.py b/tests/unitary/math/conftest.py index 35e87c3e..64bcd4b4 100644 --- a/tests/unitary/math/conftest.py +++ b/tests/unitary/math/conftest.py @@ -11,4 +11,4 @@ def math_optimized(deployer): @pytest.fixture(scope="module") def math_unoptimized(deployer): with boa.env.prank(deployer): - return boa.load("contracts/experimental/n=2.vy") + return boa.load("contracts/old/CurveCryptoSwap2Math.vy") 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..77ebc9d5 --- /dev/null +++ b/tests/unitary/math/test_newton_D_ref.py @@ -0,0 +1,190 @@ +# flake8: noqa +import sys +from decimal import Decimal + +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 = 3000000 # Increase for fuzzing +MAX_SAMPLES = 300 # Increase for fuzzing + +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**16 - 1 +MAX_XD = 10**20 + 1 + + +@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=int(1.001e16), max_value=int(0.999e20) + ), # <- ratio 1e18 * x/D, typically 1e18 * 1 + yD=st.integers( + min_value=int(1.001e16), max_value=int(0.999e20) + ), # <- 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, +): + _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, + ) + + +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, +): + + 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 + + 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: + 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}" + "}" + ) + print("broken at:", case) + raise + + try: + assert abs(result_sim - result_contract) <= max( + 10000, result_sim / 1e12 + ) + except AssertionError: + 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" + ) + print("broken at:", case) + raise diff --git a/tests/utils/simulation_int_many.py b/tests/utils/simulation_int_many.py index 179558c4..1b56ea25 100644 --- a/tests/utils/simulation_int_many.py +++ b/tests/utils/simulation_int_many.py @@ -37,6 +37,11 @@ def reduction_coefficient(x, gamma): 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