Curve Newton's Method (EN)

Reference and Innovation in Newton’s Method Code

1. The Necessity and Importance of Newton’s Method in Solidity & Vyper

In Web 1.0 and Web 2.0, the user-facing application layer ultimately captured much of the value. The blockchain world is not short of increasingly homogeneous infrastructure; it is short of application layers with differentiated functions. If Web 3.0 truly takes off, value capture may gradually shift from the base layer toward the application layer. In my view, current EVM performance limits, gas costs, and some characteristics of Solidity and Vyper discourage developers from building computation-heavy applications. One of the most difficult problems today is the tradeoff between numerical performance and resource consumption. Algorithms that are fast and efficient in Python, C++, or JavaScript often require serious compromises when implemented in Solidity or Vyper. To ship a specific feature, protocol teams may need to implement difficult numerical algorithms themselves. These implementations are rarely reusable libraries, so developers need to understand both the math and the engineering tradeoffs behind them.

Curve is exactly this kind of project. It is a protocol with many original ideas and engineering innovations, driven largely by a single founder. Curve’s importance in the blockchain world needs no elaboration.

Here I salute Michael Egorov, the founder of Curve, and thank him for the innovations he has brought to the blockchain industry. In this article, I focus on one small but important topic with three goals:

• Explain how Newton’s method is used in Curve.
• Discuss how Newton’s method can be transformed for different business scenarios.
• Encourage further research into numerical optimization algorithms.

The third point is the most important to me, for the following reason:

Newton’s method is clearly not the universally fastest or cheapest algorithm for numerical optimization. Just as no DeFi project can guarantee that it has no bugs, no numerical optimization algorithm can guarantee that it is always the fastest and cheapest choice in every DeFi system. I hope later researchers can continue to study the questions raised here. If possible, I also hope the Curve project can support research into questions such as:

• If Curve’s Newton method is not optimal, which algorithm is faster and more gas efficient?
• Can we build reusable Solidity libraries for Newton’s method, or for variants of Newton’s method, so that other developers can meet basic numerical-computation needs more easily?

2. The Code Implementation Idea of Newton’s Method in Curve V1

2.1 Introduction of Newton’s Method

Since the equation cannot be solved directly inside the contract, Curve uses Newton’s method to iteratively find an approximate solution. $x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}$ In simple terms, Newton’s method repeatedly applies the formula above to make the current value move closer and closer to the true solution.

2.2 Curve V1: Derivation of get_D with Newton’s Method

The function $f(D)$ used in Newton’s method is derived from Curve’s core invariant. The derivation is briefly described below:

1. First, transform the core formula into the form f(D) = 0.

$f(D)=An^n\Sigma x_i + D-ADn^n-\frac{D^{n+1}}{n^n\Pi x_i}$

2. D_new = D - f(D)/f'(D)

$D_{new}=D-\frac{An^n\Sigma x_i + D-ADn^n-\frac{D^{n+1}}{n^n\Pi x_i}}{1-An^n-(n+1)\frac{D^n}{n^n\Pi x_i}}$

3. It eventually becomes the form used in the code.

$D_{new}=\frac{An^n\Sigma x_i-\frac{D^{n+1}}{n^n\Pi x_i}}{An^n+(n+1)\frac{D^n}{n^n\Pi x_i}-1}$

4. Use a for loop to iterate toward the solution for D. When the difference between the new D and the previous value is less than or equal to 1, stop iterating.

5. In practice, the local optimum is usually found within about four iterations. If no suitable solution is found after 255 iterations, the program terminates and the transaction reverts. Users can still use remove_liquidity to withdraw assets.

Note: One term is missing from the Newton’s method formula in the comments, while the correct simplified formula is implemented in the code.

• curve-contract SwapTemplateBase.vy#L214

@pure
@internal
def _get_D(_xp: uint256[N_COINS], _amp: uint256) -> uint256:
    """
    D invariant calculation in non-overflowing integer operations
    iteratively

    A * sum(x_i) * n**n + D = A * D * n**n + D**(n+1) / (n**n * prod(x_i))

    Converging solution:
    D[j+1] = (A * n**n * sum(x_i) - D[j]**(n+1) / (n**n prod(x_i))) / (A * n**n + (n+1)*D[j]**n/(n**n prod(x_i)) - 1)

    Here the original code comment denominator misses one item
    (n+1)*D[j]**n/(n**n prod(x_i))

    Link of Original code :
    https://github.com/curvefi/curve-contract/blob/master/contracts/pool-templates/base/SwapTemplateBase.vy#L214
    """
    S: uint256 = 0
    Dprev: uint256 = 0

    for _x in _xp:
        S += _x
    if S == 0:
        return 0

    D: uint256 = S
    Ann: uint256 = _amp * N_COINS
    for _i in range(255):
        D_P: uint256 = D
        for _x in _xp:
            D_P = D_P * D / (_x * N_COINS)  # If division by 0, this will be borked: only withdrawal will work. And that is good
            # Only when liquidity is removed is it possible that _x will be 0
            # at which point the program will crash, which is also the result we expect
            # The value of xp in the pool cannot be 0 when adding liquidity
        Dprev = D
        D = (Ann * S / A_PRECISION + D_P * N_COINS) * D / ((Ann - A_PRECISION) * D / A_PRECISION + (N_COINS + 1) * D_P)
        # Equality with the precision of 1
        # When the difference between the new value of the iteration and the previous value is less than or equal to 1, that is, 0 or 1
        # The local optimal solution is considered to be found
        if D > Dprev:
            if D - Dprev <= 1:
                return D
        else:
            if Dprev - D <= 1:
                return D
    # convergence typically occurs in 4 rounds or less, this should be unreachable!
    # if it does happen the pool is borked and LPs can withdraw via `remove_liquidity`
    raise


@view
@internal
def _get_D_mem(_balances: uint256[N_COINS], _amp: uint256) -> uint256:
    return self._get_D(self._xp_mem(_balances), _amp)

2.3 Exploration, Innovation, and Lessons from Historical Comments on get_D

Sometimes mistakes are the best way to understand the correct version of Newton’s method. When refactoring Curve’s Vyper code into Solidity, we found that the comments here were inconsistent with the source code. Missing or inaccurate comments are relatively common in Curve’s source code.

The comments were originally written by Sam Werner, a member of the Curve team and a PhD student at Imperial College London’s Centre for Cryptocurrency Research and Engineering. At the time, we thought Sam might have proposed a better method than the one used by Curve’s founder, so we implemented a version of Newton’s method based on his comment. Due to space limitations, this variant will not be discussed further in this article.

According to Sam’s comment, the numerator of Newton’s method does not need to change. By scaling and simplifying the denominator, we change the slope of each tangent in the Newton iteration. In theory, this may still converge to the same numerical result.

The experiment showed that changing the denominator, and therefore the tangent slope, does not shift the final result very much. However, it greatly affects the convergence speed of Newton’s method. This is dangerous because both Curve V1 and Curve V2 require the iteration to reach a certain precision within 255 steps.

Convergence of Newton’s method: The global convergence of Newton’s method was not proven until 2018. Near the optimum, Newton’s method has quadratic convergence. The result is that Newton’s method needs at most about $log log (1/\varepsilon)$ steps for the distance between the approximate solution and the true optimum to become less than $\varepsilon$. A common precision standard in DeFi is $10^{18}$. Calculating $log log (1/10^{-18})$ gives a result of about 4, which is the basis for the idea that many Newton iterations can converge in roughly four steps.

The experiment shows that changing the slope in Newton’s method can damage the convergence rate enough to make the method unusable in practice.

There are many variants of Newton’s method, including direct improvements of Newton’s method (AN, HN, MN), damped Newton’s method, and quasi-Newton methods such as DFP, BFGS, and L-BFGS. These methods are safer than directly scaling the derivative in Newton’s method. Direct scaling of Newton’s method is rarely well verified.

We also considered partial scaling in Curve V2 and found that the resulting convergence speed was not ideal.

In short, although Sam Werner’s comment was imperfect, it unexpectedly pushed us toward a deeper understanding of Newton’s method and its application in Curve.

2.4 Curve V1: Derivation of get_y with Newton’s Method

Given input asset i, its updated quantity x, and the previous _xp values, the goal is to calculate the updated quantity of output asset j.

As with the calculation of D, Newton’s method is also used here to calculate y, that is, $x_j$. The derivation of $f(x_j)$ is as follows: $An^n\sum x_i + D = ADn^n + \frac{D^{n+1}}{n^n\prod x_i}$

1. Let sum' and prod' be the sum and product that exclude the output asset, respectively.

   • sum' = sum(x) - x_j
   • prod' = prod(x) / x_j

2. Divide the core formula by A*n**n

   $x_j\sum x_i+\frac{D}{An^n}-D=\frac{D^{n+1}}{An^nn^n\prod x_i}$

3. Multiply by x_j, then substitute sum' and prod'.

   $x_j(x_j+sum')-x_j(An^n-1)\frac{D}{An^n}=\frac{D^{n+1}x_j}{An^{2n}(prod'*x_j)}$

4. Expand the polynomial. This is the f(x_j) used by Newton’s method.

   $x_j^2+x_j(sum'-(An^n-1)\frac{D}{An^n})=\frac{D^{n+1}}{An^{2n}prod'}$

5. Write it in the form x_j**2 + b*x_j = c, then substitute it into Newton’s method: x = x - f(x) / f'(x).

   • x_j = (x_j**2 + c) / (2*x_j + b)

@view
@internal
def _get_y(i: int128, j: int128, x: uint256, _xp: uint256[N_COINS]) -> uint256:
    """
    Calculate x[j] if one makes x[i] = x

    Done by solving quadratic equation iteratively.
    x_1**2 + x_1 * (sum' - (A*n**n - 1) * D / (A * n**n)) = D ** (n + 1) / (n ** (2 * n) * prod' * A)
    x_1**2 + b*x_1 = c

    x_1 = (x_1**2 + c) / (2*x_1 + b)
    """
    # x in the input is converted to the same price/precision

    assert i != j       # dev: same coin
    assert j >= 0       # dev: j below zero
    assert j < N_COINS  # dev: j above N_COINS

    # should be unreachable, but good for safety
    assert i >= 0
    assert i < N_COINS

    A: uint256 = self._A()
    D: uint256 = self._get_D(_xp, A)
    Ann: uint256 = A * N_COINS
    c: uint256 = D
    S: uint256 = 0
    _x: uint256 = 0
    y_prev: uint256 = 0

    for _i in range(N_COINS):
        if _i == i:
            _x = x
        elif _i != j:
            _x = _xp[_i]
        else:
            continue
        S += _x
        c = c * D / (_x * N_COINS)
    c = c * D * A_PRECISION / (Ann * N_COINS)
    b: uint256 = S + D * A_PRECISION / Ann  # - D
    y: uint256 = D
    for _i in range(255):
        y_prev = y
        y = (y*y + c) / (2 * y + b - D)
        # Equality with the precision of 1
        if y > y_prev:
            if y - y_prev <= 1:
                return y
        else:
            if y_prev - y <= 1:
                return y
    raise


@view
@external
def get_dy(i: int128, j: int128, _dx: uint256) -> uint256:
    xp: uint256[N_COINS] = self._xp()
    rates: uint256[N_COINS] = RATES

    x: uint256 = xp[i] + (_dx * rates[i] / PRECISION)
    y: uint256 = self._get_y(i, j, x, xp)
    dy: uint256 = xp[j] - y - 1
    fee: uint256 = self.fee * dy / FEE_DENOMINATOR
    return (dy - fee) * PRECISION / rates[j]

There is also an _get_y_D() function. The difference is that it allows the caller to provide a custom value of D when solving for y.

@pure
@internal
def _get_y_D(A: uint256, i: int128, _xp: uint256[N_COINS], D: uint256) -> uint256:

3. The Code Implementation Idea of Newton’s Method in Curve V2

Newton’s method in Curve V2 is transformed in several ways, and the implementation is difficult to understand because the code has limited comments. Here, instead of introducing Newton’s method from scratch, we look directly at how the author converts mathematical ideas into code step by step. I hope this can inspire future developers to design more complex algorithms. For readability, the code shown here comes from the Python version in the test file.

3.1 Curve V2: Derivation of geometric_mean with Newton’s Method

Here, we input an array x and calculate its geometric mean as D. Let D in the code be d, and let N be n, as shown in the formula:

$d=\left(\prod _{i=0}^{n-1} x(i)\right){}^{1/n}$

Here, d is the target solved by Newton’s method. To construct a convenient function, raise both sides of the equation to the nth power:

$d^n=\prod _{i=0}^{n-1} x(i)$

Move the right side of the equation to the left to construct $f(d)$:

$f(d)=d^n-\prod _{i=0}^{n-1} x(i)$

Differentiate $f(d)$:

$f'(d)=n d^{n-1}$

According to Newton’s method, $d_{new}=d_{prev}-f(d_{prev})/f'(d_{prev})$.

$d_{\text{new}}=d-\frac{d^n-\prod _{i=0}^{n-1} x(i)}{n d^{n-1}}$

Cancel $(d^{n-1})$ in the denominator:

$d_{\text{new}}=d-\frac{d-\frac{\prod _{i=0}^{n-1} x(i)}{d^{n-1}}}{n}$

Next, rewrite the equation as a numerator divided by a denominator. This makes it easier to simplify the numerator first. The principle of multiplying before dividing also helps control precision in Vyper and Solidity.

$d_{\text{new}}=\frac{nd-d+\frac{\prod _{i=0}^{n-1} x(i)}{d^{n-1}}}{n}$

The term $\frac{\prod _{i=0}^{n-1} x(i)}{d^{n-1}}$ appears here. The numerator contains n terms from the x array, while the denominator contains n - 1 terms of d. To make both sides contain n terms and allow them to be iterated together in code, multiply the numerator and denominator by d:

$d_{\text{new}}=\frac{nd-d+\frac{d \left(\prod _{i=0}^{n-1} x(i)\right)}{d^n}}{n}$

Then simplify the numerator and extract the common factor d. This kind of simplification reduces computational cost.

$d_{\text{new}}=\frac{d \left((n-1)+\frac{\prod _{i=0}^{n-1} x(i)}{d^n}\right)}{n}$

This directly explains the code:

D = D * ((N - 1) * 10 ** 18 + tmp) // (N * 10**18)

Here, tmp is exactly the term $\frac{\prod _{i=0}^{n-1} x(i)}{d^{n}}$ from the derivation. It corresponds to the value obtained after running the loop n times. After precision scaling is added, the expression matches the author’s code.

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
    print(x)
    raise ValueError("Did not converge")

3.2 Curve V2: Derivation of newton_D with Newton’s Method

The relevant formulas in the whitepaper are as follows. The initial value for Newton’s method is $D_0$.

Equilibrium equation: $KD^{N-1}{\Large\sum}{x_i}+{\Large\prod} x_i=KD^N+\left(\dfrac{D}{N}\right)^N$

$K_0 = \dfrac{\prod x_iN^N}{D^N}$ , $K=AK_0\dfrac{\gamma^2}{(\gamma+1-K_0)^2}$

$F({x},D)=K({x},D)D^{N-1}{\Large\sum}{x_i}+{\Large\prod}(x_i)-K({x},D)D^N-\left(\dfrac{D}{N}\right)^N$

$D_0=N\left({\Large\prod}x_k\right)^{\frac{1}{N}}$

$x_{i,0}=\dfrac{D^{N-1}}{\prod_{k \neq i}x_kN^{N-1}}$

How does Newton’s method lead to intermediate variables such as mul1 and mul2? Replace D with d and follow the author’s derivation. First, fully expand the Curve V2 equilibrium equation:

$f(d)=-\frac{A \gamma ^2 n^n \prod _{i=1}^n x(i)}{\left(-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1\right){}^2}+\frac{A \gamma ^2 n^n \left(\prod _{i=1}^n x(i)\right) \sum _{i=1}^n x(i)}{d \left(-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1\right){}^2}+\prod _{i=1}^n x(i)-\left(\frac{d}{n}\right)^n$

Differentiate the above formula with respect to d:

$f'(d)=\frac{2 A \gamma ^2 n^{2 n+1} d^{-n-1} \left(\prod _{i=1}^n x(i)\right){}^2}{\left(-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1\right){}^3}-\frac{2 A \gamma ^2 n^{2 n+1} d^{-n-2} \left(\prod _{i=1}^n x(i)\right){}^2 \left(\sum _{i=1}^n x(i)\right)}{\left(-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1\right){}^3}-\frac{A \gamma ^2 n^n \left(\prod _{i=1}^n x(i)\right) \sum _{i=1}^n x(i)}{d^2 \left(-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1\right){}^2}-\left(\frac{d}{n}\right)^{n-1}$

According to Newton’s method:

$d_{new}=d-f(d)/f'(d)$

Expanding the formula gives:

$d_{new}=d-\frac{-\frac{A \gamma ^2 n^n \prod _{i=1}^n x(i)}{\left(-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1\right){}^2}+\frac{A \gamma ^2 n^n \left(\prod _{i=1}^n x(i)\right) \sum _{i=1}^n x(i)}{d \left(-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1\right){}^2}+\prod _{i=1}^n x(i)-\left(\frac{d}{n}\right)^n}{\frac{2 A \gamma ^2 n^{2 n+1} d^{-n-1} \left(\prod _{i=1}^n x(i)\right){}^2}{\left(-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1\right){}^3}-\frac{2 A \gamma ^2 n^{2 n+1} d^{-n-2} \left(\prod _{i=1}^n x(i)\right){}^2 \left(\sum _{i=1}^n x(i)\right)}{\left(-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1\right){}^3}-\frac{A \gamma ^2 n^n \left(\prod _{i=1}^n x(i)\right) \sum _{i=1}^n x(i)}{d^2 \left(-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1\right){}^2}-\left(\frac{d}{n}\right)^{n-1}}$

The denominator that appears most often is $-n^n d^{-n} \left(\prod _{i=1}^n x(i)\right)+\gamma +1$, which is the expanded form of $\gamma+1-K_0$ in the equilibrium equation. For readability and simplification, denote $\gamma+1-K_0$ as g1k0. After this replacement, the equation becomes:

$d-\frac{\frac{A \gamma ^2 n^n \left(\prod _{i=1}^n x(i)\right) \sum _{i=1}^n x(i)}{d\ \text{g1k0}^2}-\frac{A \gamma ^2 n^n \prod _{i=1}^n x(i)}{\text{g1k0}^2}+\prod _{i=1}^n x(i)-\left(\frac{d}{n}\right)^n}{\frac{2 A \gamma ^2 n^{2 n+1} d^{-n-1} \left(\prod _{i=1}^n x(i)\right){}^2}{\text{g1k0}^3}-\frac{2 A \gamma ^2 n^{2 n+1} d^{-n-2} \left(\prod _{i=1}^n x(i)\right){}^2 \left(\sum _{i=1}^n x(i)\right)}{\text{g1k0}^3}-\frac{A \gamma ^2 n^n \left(\prod _{i=1}^n x(i)\right) \sum _{i=1}^n x(i)}{d^2 \text{g1k0}^2}-\left(\frac{d}{n}\right)^{n-1}}$

At the same time, abbreviate $\prod _{i=1}^n x(i)$ as Prod, and abbreviate $\sum _{i=1}^n x(i)$ as S. The equation can then be simplified to:

$d-\frac{\frac{A \gamma ^2 n^n \text{Prod} S}{d \text{g1k0}^2}-\frac{A \gamma ^2 n^n \text{Prod}}{\text{g1k0}^2}-\left(\frac{d}{n}\right)^n+\text{Prod}}{\frac{2 A \gamma ^2 n^{2 n+1} \text{Prod}^2 d^{-n-1}}{\text{g1k0}^3}-\frac{2 A \gamma ^2 n^{2 n+1} \text{Prod}^2 S d^{-n-2}}{\text{g1k0}^3}-\frac{A \gamma ^2 n^n \text{Prod} S}{d^2 \text{g1k0}^2}-\left(\frac{d}{n}\right)^{n-1}}$

Denote the numerator and denominator as:

$\text{numerator}=\frac{A \gamma ^2 n^n \text{Prod} S}{d \text{g1k0}^2}-\frac{A \gamma ^2 n^n \text{Prod}}{\text{g1k0}^2}-\left(\frac{d}{n}\right)^n+\text{Prod}$

$\text{denominator}=\frac{2 A \gamma ^2 n^{2 n+1} \text{Prod}^2 d^{-n-1}}{\text{g1k0}^3}-\frac{2 A \gamma ^2 n^{2 n+1} \text{Prod}^2 S d^{-n-2}}{\text{g1k0}^3}-\frac{A \gamma ^2 n^n \text{Prod} S}{d^2 \text{g1k0}^2}-\left(\frac{d}{n}\right)^{n-1}$

Next, observe the complex formula. Most terms in the numerator and denominator contain $\frac{A \gamma ^2 n^n \text{Prod}}{\text{g1k0}^2}$ as a common factor. Meanwhile, $K=\frac{A \gamma ^2 n^n \text{Prod} d^{-n}}{\text{g1k0}^2}$ in the equilibrium equation is very similar to this factor. Since K contains K0, and K0 contains both D and an N-step loop, we can force K into the expression as a new common factor:

$denominator= \frac{A \gamma ^2 n^n d^{-n} \text{Prod}}{\text{g1k0}^2}*(-\frac{\text{g1k0}^2 n^{1-n} d^{n-1}}{A \gamma ^2 n^n \text{Prod}}-\frac{2 n^{n+1} \text{Prod} S d^{-n-2}}{\text{g1k0}}+\frac{2 n^{n+1} \text{Prod} d^{-n-1}}{\text{g1k0}}-\frac{S}{d^2})$ The factor $\text{Prod}*n^n d^{-n}=K0$ appears in the denominator, so the similar K0 structure in the denominator can be rewritten as K0:

$denominator= K*(-\frac{\text{g1k0}^2 n}{d \left(A \gamma ^2 \text{K0} n^n\right)}-\frac{2 \text{K0} n S}{d^2 \text{g1k0}}-\frac{S}{d^2}+\frac{2 \text{K0} n}{d \text{g1k0}})$

The denominator also contains d, so we consider eliminating that denominator. Notice that A and $n^n$ always appear together, so merge them into ANN = A n^n. This gives:

$denominator= Kd^2*(-\frac{d \text{g1k0}^2 n}{ANN \gamma ^2 \text{K0} }+\frac{2 d \text{K0} n}{\text{g1k0}}-\frac{2 \text{K0} n S}{\text{g1k0}}-S)$

Observe again that the formula contains many instances of $\frac{2 \text{K0} n }{\text{g1k0}}$. Mark this as mul2. At the same time, mark the first more complex part that does not involve K0 or a loop as $mul1=\frac{\text{g1k0}^2}{\text{ANN} \gamma ^2}$.

$denominator=-Kd^2*({S+S*mul2+mul1} \frac{n}{K0} - {mul2}*d)$

Denote -denominator / Kd^2 as neg_prime. Similarly, the numerator can be transformed by Kd^2, and a similar result can be obtained by substituting it into:

$d_{new}=d-f(d)/f'(d)$ $d_{new} = (d * \text{neg\_fprime} + d * S - d^2)/\text{neg\_fprime} - d * ({mul1}/\text{neg\_fprime}) * (1-K0)/ K0$

This gives the same core Newton update formula used in the code.

def newton_D(A, gamma, x, D0):
    D = D0
    i = 0

    S = sum(x)
    x = sorted(x, reverse=True)
    N = len(x)

    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")

3.3 Curve V2: Derivation of newton_y with Newton’s Method

There is an error in the initial value stated in the Curve V2 whitepaper. To make the solution converge faster, the initial values for Newton’s method on D and $x_i$ can be set as follows:

$D_0=N(\prod{x_k})^{\frac{1}{N}}$

$x_{i,0}=\frac{D^{N}}{\prod_{k!=i}{x_k}N^{N}}$

The initial value definition in the whitepaper is incorrect. The power of D in the numerator is written as N - 1, and the power of N in the denominator is also written as N - 1. Both should be N, and this is corrected in the code as x_i,0.

Incorrect initial value definition in the whitepaper: $x_{i,0}=\frac{D^{N-1}}{\prod_{k \neq i}{x_k}N^{N-1}}$

How does Newton’s method lead to intermediate variables such as mul1 and mul2? Continue following the author’s derivation. Similarly, fully expand the Curve V2 equilibrium equation. For convenience, assume that x(n) is the unknown value to be solved by newton_y. Expand the equilibrium equation and simplify the analysis as a function of $x(n)$:

$f(x(n))=-\frac{A \gamma ^2 n^n x(n) \prod _{k=1}^{n-1} x(k)}{\left(-n^n d^{-n} x(n) \left(\prod _{k=1}^{n-1} x(k)\right)+\gamma +1\right){}^2}+\frac{A \gamma ^2 n^n \left(x(n) \prod _{k=1}^{n-1} x(k)\right) \left(\sum _{k=1}^{n-1} x(k)+x(n)\right)}{d \left(-n^n d^{-n} x(n) \left(\prod _{k=1}^{n-1} x(k)\right)+\gamma +1\right){}^2}+x(n) \left(\prod _{k=1}^{n-1} x(k)\right)-\left(\frac{d}{n}\right)^n$

Differentiate $f(x(n))$ to get:

$f'(x(n))= -\frac{A \gamma ^2 n^n \prod _{k=1}^{n-1} x(k)}{\left(-n^n d^{-n} x(n) \left(\prod _{k=1}^{n-1} x(k)\right)+\gamma +1\right){}^2}+\frac{A \gamma ^2 n^n x(n) \left(\prod _{k=1}^{n-1} x(k)\right)}{d \left(-n^n d^{-n} x(n) \left(\prod _{k=1}^{n-1} x(k)\right)+\gamma +1\right){}^2}+\frac{A \gamma ^2 n^n \left(\prod _{k=1}^{n-1} x(k)\right) \left(\sum _{k=1}^{n-1} x(k)+x(n)\right)}{d \left(-n^n d^{-n} x(n) \left(\prod _{k=1}^{n-1} x(k)\right)+\gamma +1\right){}^2}-\frac{2 A \gamma ^2 n^{2 n} d^{-n} x(n) \left(\prod _{k=1}^{n-1} x(k)\right){}^2}{\left(-n^n d^{-n} x(n) \left(\prod _{k=1}^{n-1} x(k)\right)+\gamma +1\right){}^3}+\frac{2 A \gamma ^2 n^{2 n} d^{-n-1} x(n) \left(\prod _{k=1}^{n-1} x(k)\right){}^2 \left(\sum _{k=1}^{n-1} x(k)+x(n)\right)}{\left(-n^n d^{-n} x(n) \left(\prod _{k=1}^{n-1} x(k)\right)+\gamma +1\right){}^3}+\prod _{k=1}^{n-1} x(k)$

Similar to newton_D, we can improve readability by introducing substitutions.

Denote $-n^n d^{-n} x(n) \left(\prod _{k=1}^{n-1} x(k)\right)+\gamma +1$ as g1k0.

Denote $A*n^n$ as ANN.

Denote $\sum _{k=1}^{n-1} x(k)$ as S.

The formula can be simplified as:

$f'(x(n))=-\frac{2 \text{ANN} \gamma ^2 n^n d^{-n} x(n) \left(\prod _{k=1}^{n-1} x(k)\right){}^2}{\text{g1k0}^3}+\frac{2 \text{ANN} \gamma ^2 n^n d^{-n-1} x(n) \left(\prod _{k=1}^{n-1} x(k)\right){}^2 \left(\sum _{k=1}^{n-1} x(k)+x(n)\right)}{\text{g1k0}^3}+\frac{\text{ANN} \gamma ^2 x(n) \left(\prod _{k=1}^{n-1} x(k)\right)}{d \text{g1k0}^2}+\frac{\text{ANN} \gamma ^2 \left(\prod _{k=1}^{n-1} x(k)\right) \left(\sum _{k=1}^{n-1} x(k)+x(n)\right)}{d \text{g1k0}^2}-\frac{\text{ANN} \gamma ^2 \prod _{k=1}^{n-1} x(k)}{\text{g1k0}^2}+\prod _{k=1}^{n-1} x(k)$

Similar to the newton_D function, observation shows that K is useful here. Denote $\prod _{k=1}^{n-1} x(k)$ as $Prod\_i$ to obtain:

$f'(x(n))=K \left(\frac{\text{g1k0}^2 n^{-n} d^n}{\text{ANN} \gamma ^2 x(n)}+\left(\frac{\frac{n^{-n} \text{S\_i} d^n}{x(n)}+2 n^{-n} d^n}{d}-\frac{n^{-n} d^n}{x(n)}+\frac{2 \text{Prod\_i} S}{d \text{g1k0}}-\frac{2 \text{Prod\_i}}{\text{g1k0}}\right)\right)$

Notice that $Prod\_i$ often appears in the numerator, and that $Prod\_i*x(n)=Prod$. Therefore, transform the equation by $\frac{1}{x(n)}$ to simplify it:

$f'(x(n))=\frac{K}{x(n)}(\frac{\text{g1k0}^2 n^{-n} d^n}{\text{ANN} \gamma ^2}+\frac{n^{-n} S d^n+2 n^{-n} d^n x(n)}{d}-\frac{n^{-n} d^n x(n)}{d}+n^{-n} \left(-d^n\right)+\frac{2 \text{Prod} S}{d \text{g1k0}}-\frac{2 \text{Prod}}{\text{g1k0}})$

Continuing the observation, the formula above contains many $d^n n^{-n}$ terms. After transformation, we obtain:

$f'(x(n))=\frac{\left(K n^{-n} d^n\right)}{x(n)} \left(\frac{\text{g1k0}^2}{\text{ANN} \gamma ^2}+\frac{2 \text{K0} S}{d \text{g1k0}}+\frac{x(n)+S}{d}-\frac{2 \text{K0}}{\text{g1k0}}-1\right)$

The denominator contains many d terms, so move $1/d$ outside the parentheses:

$f'(x(n))=\frac{\left(K n^{-n} d^{n-1}\right)}{x(n)} \left(\frac{\text{g1k0}^2 d}{\text{ANN} \gamma ^2}+S+\frac{2 \text{K0} S}{ \text{g1k0}}+x(n)-\frac{2 \text{K0}d}{\text{g1k0}}-d\right)$

Observe again that the formula contains many S terms. Combining their coefficients gives $mul2=1+ \frac{2 \text{K0} }{\text{g1k0}}$, which we mark as mul2. At the same time, mark the first more complex part that does not involve K0 or a loop as: $mul1=d*\frac{\text{g1k0}^2}{\text{ANN} \gamma ^2}$

Denote x(n) as the y solved by newton_y.

$f(x(n))=\frac{\left(K n^{-n} d^n\right)}{y}*( (y + S * mul2 + mul1 - D * mul2))$ Denote the part in parentheses in the above formula as yfprime.

$\text{yfprime}=(y + S * mul2 + mul1 - D * mul2))$

$\text{fprime}=\text{yfprime}/y$

This is consistent with the core iteration in the code.

Similarly, the numerator can be transformed by Kd^2, and a similar result can be obtained by substituting it into:

$y_{new}=y-f(y)/f'(y)$

$y_{new} = (yfprime + D - S) / fprime + mul1 / fprime * (1 - K0) / K0$

This gives the same core Newton update formula used in the code.

def newton_y(A, gamma, x, D, i):
    N = len(x)

    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")

4. Applicability of Newton’s Method and Summary of Code Ideas

4.1 The Applicability of Newton’s Method in DeFi

Newton’s method is a broadly useful numerical algorithm. If a good initial value is available, it can solve many equations without analytical solutions in about ten iterations. Therefore, it is suitable for many complex DeFi numerical algorithms, and I believe it will be used by more protocols over time.

Promoting Newton’s method also helps more developers understand Curve V1 and V2, which is one of the original motivations for this research. As shown above, applying Newton’s method in production code has a much higher threshold than learning it from a textbook.

4.2 Summary of Newton’s Method in Vyper and Solidity

• Choose the right common factors and intermediate variables first. Some intermediate variables in Newton-method implementations are generated through N iterations. These variables should be chosen together with the forms that appear frequently during algebraic simplification. For example, in newton_D, K0 and K appear often, so they should be prioritized as intermediate variables instead of repeatedly using the $n^n$ or $d^n$ terms contained inside K0.
• Multiply before dividing. Minimize direct division in numerator/denominator form. Reduce fractions to a common denominator and turn as much of the computation as possible into multiplication. At the same time, pay close attention to EVM overflow when multiplying large values.
• Derivation can change precision requirements. Some algebraic transformations reduce precision. The code must manually adjust scaling factors so that all values remain at consistent precision. Only values at the same precision can reproduce the intended mathematical formula.
• Replace frequently occurring forms with simple symbols.
• Prevent numerical overflow and process larger values first. When executing Newton’s method in the EVM, division, precision, and continuous multiplication loops must be balanced carefully to avoid overflow. In general, large-number divisions should be performed before small-number divisions, so data often needs to be rearranged from large to small.

4.3 Difficulties in understanding the code of Curve V2

Compared with Uniswap V3, Curve V2 is harder to understand for several reasons:

1. Curve leans more heavily on mathematical optimization, and much of its DeFi logic is original in the digital-token ecosystem.
2. More than 95% of Curve V2’s code was contributed by Michael, Curve’s founder. The code has few comments, and some comments written by others contain errors.
3. Some key points in the whitepaper are misstated, which makes the code harder to understand.
4. Compared with Uniswap, Curve has less third-party research available as reference material.

4.4 Other algorithms of Curve and precision control problems

Due to space limitations, this article does not provide a complete description of every use of Newton’s method in Curve. Another difficulty in implementing Newton’s method in production code is precision control, which deserves a separate detailed article.

There are also other algorithms in the Curve V2 math library that lack sufficient comments and are difficult to understand. For example, the halfpow function is actually an application of the binomial algorithm in a non-integer domain.

Through this work, we gained a deeper understanding of the mathematical algorithms and specific implementation details in Curve V1 and V2, as well as the business logic around those contracts.

At present, 0xMC has completed a Solidity-based reconstruction of the Curve V1 core contract.

The Curve V1 and Curve V2 code review sections have also been published on 0xreviews.xyz and 0xmc.com.