Curve v2 CryptoSwap: math lib

Curve v2 CryptoSwap 详解系列：

1. Curve v2 CryptoSwap: white paper
2. Curve v2 CryptoSwap: add and remove liquidity
3. Curve v2 CryptoSwap: exchange and fee
4. Curve v2 CryptoSwap: repegging
5. Curve v2 CryptoSwap: math lib

CurveCryptoMath 是 CurveCryptoSwap 的数学依赖库，其中包含了 curve 运算关键的牛顿法迭代求解，以及几何平均、 非整数幂运算等算法实现，非常值得学习。

newton_D & newton_y

关于牛顿法的详细解析及推导过程，请查看我们的另一篇文章

https://0xreviews.xyz/2022/02/28/curve-newton-method/

geometric_mean

计算几何平均的函数，利用牛顿法求解。

常规的几何平均计算方式为 (x[0] * x[1] * ...) ** (1/N), 但 EVM 并不支持开根号运算。于是这里利用牛顿法求解。

令解 D^N = x[0] * x[1] * ... = prod(x), 写成 f(D) = 0 的形式然后代入牛顿法公式中。

• 函数 f(D)=D^N - prod(x)
• 导函数 f'(D)=N*D^(N-1)

于是牛顿法迭代的公式如下

$D=D-\frac{(D^N-\prod{x_i})}{ND^{N-1}}$

通分

$\frac{ND^N-D^N+\prod{x_i}}{ND^{N-1}}=\frac{(N-1)D^N+\prod{x_i}}{ND^{N-1}}$

分子分母同乘 D 再化简

$D\frac{(N-1)D^N+\prod{x_i}}{ND^N}=D((N-1)+\frac{\prod{x_i}}{ND^N})$

最后得出迭代使用的公式

• 需要注意的是 D 和 x_i 的精度都是 10**18，所以 (N-1) 和 N 也要乘以 10**18

@internal
@view
def _geometric_mean(unsorted_x: uint256[N_COINS], sort: bool = True) -> uint256:
    """
    (x[0] * x[1] * ...) ** (1/N)
    """
    x: uint256[N_COINS] = unsorted_x
    if sort:
        x = self.sort(x)
    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-1) + prod/D^N) / N
        # (N-1)'s and N's decimal must be increased to 10**18,
        # then it can be calculate with tmp
        D = D * ((N_COINS - 1) * 10**18 + tmp) / (N_COINS * 10**18)
        if D > D_prev:
            diff = D - D_prev
        else:
            diff = D_prev - D
        # Exit conditions for Newton's method
        if diff <= 1 or diff * 10**18 < D:
            return D
    raise "Did not converge"


@external
@view
def geometric_mean(unsorted_x: uint256[N_COINS], sort: bool = True) -> uint256:
    return self._geometric_mean(unsorted_x, sort)

reduction_coefficient

用于计算 Dynamic fees 公式中的 g 的函数。

Dynamic fees:

$f=g \cdot f_{mid}+(1-g)\cdot f_{out}$

$g=\frac{\gamma_{fee}}{\gamma_{fee}+1-\frac{\prod{x_i}}{(\sum{x_i}/N)^N}}$

• K = prod(x) / (sum(x) / N)**N 这里的 K 与平衡方程的 K 不同
  • 实际代码中的运算顺序与表达式不同
  • 考虑到 EVM 中除法的精度问题，需要先乘后除，但又因为分母有 N 次方运算，可能会有数值溢出的风险，于是调整为如下形式
  • x_0 * N / sum + x_1 * N / sum + x_2 * N / sum
• K 的初值为 10**18 ，不仅为了和 fee_gamma 统一精度，再做上述运算时也可以提高除法的精度
• 最终结果 g 直接使用 K 作为变量，节省本地变量数量，防止 EVM 本地变量过多的报错

@external
@view
def reduction_coefficient(x: uint256[N_COINS], fee_gamma: uint256) -> uint256:
    """
    fee_gamma / (fee_gamma + (1 - K))
    where
    K = prod(x) / (sum(x) / N)**N
    (all normalized to 1e18)
    """
    K: uint256 = 10**18
    S: uint256 = 0    # sum(xp)
    for x_i in x:
        S += x_i

    # Could be good to pre-sort x, but it is used only for dynamic fee,
    # so that is not so important
    for x_i in x:
        K = K * N_COINS * x_i / S
    if fee_gamma > 0:  # return k when fee_gamma is zero
        K = fee_gamma * 10**18 / (fee_gamma + 10**18 - K)
    return K

sort

利用冒泡排序法将数组的数值从高到低排列，详细原理参见 Bubble sort wiki。

@internal
@pure
def sort(A0: uint256[N_COINS]) -> uint256[N_COINS]:
    """
    Insertion sort from high to low
    """
    A: uint256[N_COINS] = A0
    for i in range(1, N_COINS):
        x: uint256 = A[i]
        cur: uint256 = i
        for j in range(N_COINS):
            y: uint256 = A[cur-1]
            if y > x:
                break
            A[cur] = y
            cur -= 1
            if cur == 0:
                break
        A[cur] = x
    return A

halfpower

curve 实现了一个可以高效求解非整数幂的函数, 用于计算 oracle 公式中的 alpha，(EVM 只能做整数幂的运算)。

oracle.alpha: $\alpha=2^{-\frac{t}{T_{1/2}}}=(\frac{1}{2})^u$

该函数将幂分成整数部分和小数部分运算。整数幂部分直接使用 vyper 的指数运算实现。

小数部分运用了二项式定理来实现近似计算：

$(a+b)^u=\sum_{k=0}^{+\infty}\frac{u(u-1)...(u-k+1)}{k!}a^kb^{u-k}$

这里我们令 a=0.5, b=0，关于 b 的项 b^(u-k)， 都是是 0 的指数形式，即都为 1，可以省略，则展开式为：

$(0.5+0)^u=1*(0.5)^0+\frac{u}{1!}(0.5)^1+\frac{u(u-2+1)}{2!}(0.5)^2+\frac{u(u-2+1)(u-3+1)}{3!}(0.5)^3+...$

• 注意当 k=0 时，系数部分 $\sum_{k=0}^{+\infty}\frac{u(u-1)...(u-k+1)}{k!}$ 应该为 1；因为系数分子中 u-k+1 这一项在 k=0 时是 u+1，显然不符合二项式定理公式中阶乘的形式，我们在编程时需要把 k=0 的情况直接让系数赋值为 1；
• 循环实现阶乘时，需要考虑符号问题
• 理论上项越多，结果将越逼近真实解，但实际情况中我们可以设置一个精度上限，超过该精度的项可以忽略，没有继续运算的必要

@external
@view
def halfpow(power: uint256, precision: uint256) -> uint256:
    """
    1e18 * 0.5 ** (power/1e18)

    Inspired by: https://github.com/balancer-labs/balancer-core/blob/master/contracts/BNum.sol#L128
    """
    intpow: uint256 = power / 10**18            # Separate the integer part of the exponent
    otherpow: uint256 = power - intpow * 10**18 # Separating the fractional part of the exponent is the exponent u
    # If intpow > 59, can view it as zero
    # because (1/2)**59=1.734723476 * 10**-18 that's very small
    if intpow > 59:
        return 0
    result: uint256 = 10**18 / (2**intpow)  # The integer part of the result is cached first
    if otherpow == 0:                       # If the decimal part is 0, return the result.
        return result

    term: uint256 = 10**18              # Caches the iterated items of the expansion, starting at 1 (first term)
    x: uint256 = 5 * 10**17             # X in the formula, 0.5 of 18 presision
    S: uint256 = 10**18                 # Because the first term is 1, the initial value of S is also 1
    neg: bool = False                   # Indicates whether there is a minus sign in the final factorial expansion. The default is positive

    # Iterate over the binomial expansion and sum them
    for i in range(1, 256):
        K: uint256 = i * 10**18        # The Angle k in the formula
        # (u-k+1) in the factorial, initialize it as k - 1
        c: uint256 = K - 10**18
        # (u-(k-1)) need to judge the sign
        if otherpow > c:            # u > (k-1)
            c = otherpow - c        # c = u - (k-1) = (u-k+1)
            # Set neg to true, and the rest of the loop will remain true
            # Because u is less than 1 and after the first loop, (u-k+1) will always be negative
            neg = not neg
        else:                       # u <= (k-1)
            c -= otherpow           # c = (k-1) - u = abs(u-k+1)
        # calculate i+1 term
        # The logic of *= corresponds to the factorial operation
        # term *= (u-k+1)/k * x
        # c and x and k are all 18 precision, and if you don't fill in the division of 18 precision,
        # it will lead to inconsistent accuracy
        # And term *= is getting bigger and bigger, so the precision division is performed before /k
        # to prevent overflow during the loop.
        term = term * (c * x / 10**18) / K
        # Add up the results based on the plus and minus signs
        if neg:
            S -= term
        else:
            S += term
        # Stop iterating if the iteration item is too small (less than the minimum precision required)
        # and return the result
        if term < precision:
            return result * S / 10**18

    raise "Did not converge"

sqrt_int

牛顿法求解开方运算

$z^2=x$

其中 x 是已知量，写成牛顿法函数形式

$z^2-x=0$

迭代公式 z -= f(z)/f'(z)

$z=z-\frac{z^2-x}{2z}=\frac{2z^2-z^2+x}{2z}$

$z=(\frac{x}{z}+z)*\frac{1}{2}$

@external
@view
def sqrt_int(x: uint256) -> uint256:
    """
    Originating from: https://github.com/vyperlang/vyper/issues/1266
    """

    if x == 0:
        return 0

    z: uint256 = (x + 10**18) / 2
    # By assigning y to x, if x = 1, then we can
    # directly get result in the first loop.
    # sqrt(1) = 1
    y: uint256 = x

    for i in range(256):
        # If the results of two iterations are the same, the optimal solution is obtained
        if z == y:
            return y
        y = z # y is used to cache the results of the last iteration
        z = (x * 10**18 / z + z) / 2

    raise "Did not converge"

至此，我们已经将 Curve CryptoSwap 的主要代码逻辑全部解析完毕！