Silkswap: An asymmetric automated market maker model for stablecoins

Under the CPMM model, the quantities $x$ and $y$ must satisfy the following equation:

$$x\,p\,y=K$$

(3)

with $K>0$ constant. The equilibrium point, satisfying both (2) and (3), is $(x,y)=(\sqrt{K},\frac{\sqrt{K}}{p})$. Under the CSMM model instead, the following equation holds:

$$x\,+\,py=D$$

(4)

with $D>0$ constant. The equilibrium point, satisfying (2) and (4), is $(x,y)=(\frac{D}{2},\frac{D}{2p})$. The equilibrium point is a fundamental point where the ratio of the quantities of tokens in the pool is equal to $p$, and any invariant curve should contain this point. Since this point is unique, it follows that $K=\frac{D^{2}}{4}$. In figure 1 we can see the graph of these two models, with different values of $p$.

Now let us consider the following linear combination of the CSMM and CPMM models. Let us multiply Eq. (4) by $\chi AD$, and sum it with eq. (3):

$$(\chi AD)\,(x+py)+xpy\,=\,(\chi AD)D+\frac{D^{2}}{4}$$

(5)

where $\chi$ is a function of $x$ and $y$, defined as

$$\chi:=\biggl{(}\frac{4xpy}{D^{2}}\biggr{)}^{\gamma}$$

(6)

and

$$\gamma:=\begin{cases}\gamma_{1},&\mbox{if }x\leq py\\ \gamma_{2},&\mbox{if }x>py\end{cases}$$

(7)

and the other parameters are constants satisfying $A>0$, $\gamma_{1}\geq 0$, $\gamma_{2}\geq 0$. Let us remark that $\chi$ is an adimensional quantity and $D$ has the same dimension of $x$.
We call the equation (5), together with (6) and (7), the Silkswap invariant.
The graph of the Silkswap invariant is shown in figures 3 and 3 under different set of parameters. We can see that the CPMM hyperbola is always greater than the CSMM line, and they touch each other at the equilibrium point. The Silkswap invariant graph lies in between them. In the Appendix we prove these facts in the theorems (A.1) and (A.3). The parameter $A$ tells us how close the Silkswap graph is to the CPMM or to the CSMM. From Eq. (5) we can easily see that for $A\to 0$ the invariant converges to the CPMM curve, while for $A\to\infty$ it converges to the CSMM curve.²²2The parameter $A$ has the same meaning of the parameter $A$ in the Stableswap model [Egorov, 2019].
The two dimensional function $\chi$ in (6) is a discontinuous function. We show the surface plot in Fig. (4). In the Appendix we prove that, although the Silkswap invariant contains a discontinuous function, the points of our interest are regular points where the invariant is continuously differentiable.

In order to compute the price of a token, it is convenient to express the invariant in the explicit form $y=f(x)$, where $f:\mathbb{R}_{>0}\to\mathbb{R}_{>0}$ is continuously differentiable, decreasing and convex ³³3The set $\mathbb{R}_{>0}$ is the set $\{x\in\mathbb{R}|x>0\}$.. For the CPMM the explicit form is the hyperbola $f(x)=\frac{K}{px}$, while for the CSMM the explicit form is the straight line $f(x)=\frac{1}{p}(-x+D)$, defined for $0\leq x\leq D$.

If we trade an infinitesimal quantity $dy$, using the first order Taylor approximation $dy\approx\frac{df(x)}{dx}dx$, it results equal to trade the quantity $\frac{df(x)}{dx}dx$. Let us recall that by definition $\frac{df(x)}{dx}<0$, so we need to use the absolute value to define a positive price. We define the current price of the token $X$ as function of $x$

$$P_{X}(x)\,:=\,\biggl{|}\frac{df(x)}{dx}\biggr{|}\,=\,\biggl{|}\frac{dy}{dx}\biggr{|}.$$

(8)

The price of the token $Y$ is defined as:

$$P_{Y}(x)\,:=\,\biggl{|}\frac{dx}{dy}\biggr{|}=\frac{1}{P_{X}}.$$

(9)

Under the CPMM model, the current price is $P_{Y}(x)=\frac{px^{2}}{K}$, and at equilibrium $P_{Y}(\sqrt{K})=p$. Under the CSMM model, the current price is a constant value, $P_{Y}(x)=p$, for any $0<x<D$.

If we consider the function (10) depending also on $D$, for any constant $c>0$ we have

$$F(x,y,D)=0\quad\Longrightarrow\quad F(cx,cy,cD)=0.$$

(15)

This means that the Silkswap invariant is also invariant by scaling.
This property is very useful in practice because it reduces the convergence time of some numerical methods and prevents possible overflows when the variables have very high magnitudes.
Let us introduce the variable $z:=py$. We can consider $F(cx,cy,cD)$ with $c=\frac{1}{D}$, equal to $F(\frac{x}{D},\frac{y}{D},1)$, and define the corresponding scaled function:

$$\tilde{F}(\tilde{x},\tilde{z})\,:=\,A\,\bigl{(}4\tilde{x}\tilde{z}\bigr{)}^{\gamma}\,(\tilde{x}+\tilde{z}-1)+\tilde{x}\tilde{z}-\frac{1}{4}.$$

(16)

where $\tilde{x}=\frac{x}{D}$, and $\tilde{z}=\frac{z}{D}$. We show in Fig. 6 the graph of this function. The scaled Silkswap invariant is given by

$$\tilde{F}(\tilde{x},\tilde{z})\,=\,0.$$

(17)

The function (16) is related with (10) by

$$\tilde{F}(\tilde{x},\tilde{z})=\frac{1}{D^{2}}F(x,y).$$

(18)

Using the chain rule

	$\displaystyle\frac{\partial\tilde{F}(\tilde{x},\tilde{z})}{\partial\tilde{x}}\,=\,\frac{\partial\tilde{F}(\tilde{x},\tilde{z})}{\partial x}\,\frac{dx}{d\tilde{x}}\,=\,\frac{1}{D}\frac{\partial F(x,y)}{\partial x}$		(19)
	$\displaystyle\frac{\partial\tilde{F}(\tilde{x},\tilde{z})}{\partial\tilde{z}}\,=\,\frac{\partial\tilde{F}(\tilde{x},\tilde{z})}{\partial y}\,\frac{dy}{d\tilde{z}}\,=\,\frac{1}{pD}\frac{\partial F(x,y)}{\partial y},$		(20)

we can rewrite the price equation (14) as

$$\frac{dy}{dx}=-\frac{1}{p}\;\frac{\frac{\partial\tilde{F}(\tilde{x},\tilde{z})}{\partial\tilde{x}}}{\frac{\partial\tilde{F}(\tilde{x},\tilde{z})}{\partial\tilde{z}}}\,.$$

(21)

In the numerical calculation of the swap amount, we will make use of the scaled function (16), rather than (10), because it reduces a lot the number of operations, and consequently the run time and gas fees.

Since this model is defined by an implicit equation, we need numerical methods to compute the variables of interest. If the amounts of tokens in the pool satisfy the equilibrium equation 2, then the parameter $D$ can be quickly computed from Eq. (4). However, it is quite rare that a liquidity pool is in perfect equilibrium and the equation 2 is almost never satisfied. In this cases we need to compute $D$ from Eq. (5) using a numerical method. Let us define $F(D)$ as the function (10) when we consider $D$ variable and $x$, $y$ fixed. With an abuse of notation we can call it ”invariant function dependent of $D$”, but let us recall that the name invariant can be used only when $F(D)=0$. In Fig 8, we can see that $F(D)$ is a smooth decreasing function. We compare three different numerical methods: the Newton method, Halley method and the bisection method, see Table 1.

The idea of using two times arithmetic mean (AM) and geometric mean (GM) as starting points for Newton and Halley methods comes from Theorem A.4. We can see that 2AM is a better choice. Although Newton method has more loop iterations than Halley, it performs better in terms of time. The reason is that the calculation of the second derivative is more expensive, in terms of operations, than the extra iterations. The expression of the derivatives are in Appendix B. Bisection method, as expected from theoretical results, is the slowest.

We will also make use of numerical methods to calculate the amount of tokens that are returned during a swap. When a quantity $\Delta y$ is introduced into the pool, we need to compute the quantity $\Delta x$ that is extracted from the pool, and such that the point $(x-\Delta x,y+\Delta y)$ belongs to the Silkswap invariant. The parameter $D$ is fixed, and it is irrelevant for this calculation, therefore in order to simplify the problem, we can use the scaled function (16). Let us define $\tilde{F}(\tilde{x}|\tilde{y})$ and $\tilde{F}(\tilde{y}|\tilde{x})$ the function with $\tilde{y}$ or $\tilde{x}$ fixed respectively. Again, with an abuse of notation we can call these functions ”invariant functions depending on $\tilde{x}$ or $\tilde{y}$”. In the following we consider the case when $\Delta x$ is extracted and therefore we need to find the zero of $\tilde{F}(\tilde{x}|\tilde{y})$, but the same analysis works for $\tilde{F}(\tilde{y}|\tilde{x})$. In the Figures 8 and 8, we can see that the function is discontinuous and can have different shapes depending on the choice of the parameters. So the numerical convergence to the right solution is not always guaranteed.
By definition, the parameter $\gamma$ can be any non-negative real number. In practice we restrict $\gamma$ to be a positive integer value in order to simplify the model. The value of $\gamma_{1}$ is particularly important because it controls the behavior of the left tail:

	$\displaystyle\lim_{\tilde{x}\to-\infty}\tilde{F}(\tilde{x}|\tilde{y})\to-\infty\quad\text{ for even }\gamma_{1}$
	$\displaystyle\lim_{\tilde{x}\to-\infty}\tilde{F}(\tilde{x}|\tilde{y})\to+\infty\quad\text{ for odd }\gamma_{1},$

as we can see if we compare the plot of $\tilde{F}(\tilde{x}|\tilde{y})$ in Fig. 8 with the plots in Fig. 8. When $\gamma_{1}$ is odd, $\tilde{F}(\tilde{x}|\tilde{y})$ has two zeros. Since $\tilde{x}>0$ by definition, we know that the correct solution must be positive. However, numerical methods such as Newton and Halley can converge to the wrong solution if initialized with an unlucky starting point, see Fig. 8 (LEFT). Under some set of parameters, it is possible that the function $\tilde{F}(\tilde{x}|\tilde{y})$ is increasing for $\tilde{x}>0$, and the Newton method works fine, see Fig. 8 (RIGHT). However, to avoid this uncertainty, in these cases it is better to use numerical methods that always guarantee the convergence to the right solution, such as the bisection algorithm. If we call $(\tilde{x}_{0},\tilde{y}_{0})$ the amounts in the pool before the swap, and $(\tilde{x}_{1},\tilde{y}_{1})$ the amounts after the swap such that $\tilde{F}(\tilde{x}_{1}|\tilde{y}_{1})=0$, then we know that $\tilde{x}_{1}\in[0,\tilde{x}_{0}]$. The swap amount is therefore $\Delta x=D(\tilde{x}_{0}-\tilde{x}_{1})$.
In practical applications we want to take advantage of the speed of the Newton method, and therefore we choose to use only odd values for $\gamma_{1}$. For comparisons between the three numerical methods under consideration, see Table 2. Let us comment a few points:

1. 1.

   The number of iterations of the bisection method does not depend on the size of the pool. The reason is that we are searching in the interval $[0,\tilde{x}_{0}]$, and $\tilde{x}_{0}=\frac{x}{D}$. In the example, since we are at equilibrium $D=2x$ and therefore $\tilde{x}_{0}=\frac{1}{2}$.

2. 2.

   The number of iterations of the bisection method does not depend on the size of the trade either.

3. 3.

   As for the calculation of $D$, Newton is slightly faster than Halley. The computation of second order derivatives (see Appendix B) is expensive.

4. 4.

   The number of iterations of the Newton method increases when the size of the trade is big with respect to the size of the pool. In this case the Newton method is the one that performs worse while the bisection method is the one that performs best. However, in practice it is very unlikely to see transactions bigger than the size of the pool.

5. 5.

   We used $\tilde{x}_{0}$ i.e. the value $\frac{x}{D}$ before the swap, as an initial guess. We tried also with $\frac{\text{2AM}}{D}$ and $\frac{\text{2GM}}{D}$ as initial guesses, but the performances are worse.

When considering an even $\gamma_{1}$, it turns out that the Newton method is superior for both the calculations of $D$ and the swap size. In production, we decided to use Newton as main solver, and bisection as fallback in case of failure.

We tested the model with values of $A$ ranging from 1 to $10^{5}$ and values of gamma from $1$ to $75$, under a variety of conditions including pool sizes from 1 to $10^{11}$ USDC total value, with a similar range of trade sizes. Generally, the Rust implementation can handle trade sizes several orders of magnitude above the size of the pool without overflowing, unless the values for $\gamma_{1}$ or $\gamma_{2}$ are excessively high. Since $\gamma_{1}$ and $\gamma_{2}$ are exponents, large values quickly create unmanageable numbers.

We implemented the Silkswap algorithm inside Rust smart contracts on the Secret Network blockchain. Smart contract do not allow floating point calculations, and this required the use of a few workarounds. Since smart contracts only support integer arithmetics, we stored most numbers as $10^{18}$ larger than their actual value. This effectively allowed us to store 18 decimal places. Variables are stored as uint256 i.e. unsigned 256 bits integers. Additionally, we paired each variable with a boolean representing its sign, allowing us to calculate both positive and negative numbers. The square root in 2GM is computed by the Babilonian method.

In the previous sections we have seen how the parameters $A$, $\gamma_{1}$ and $\gamma_{2}$ modify the shape of the Silkswap invariant. Specifically, parameter $A$ controls the flatness of the curve in the area around the equilibrium point (see Fig 3), while the parameters $\gamma_{1}$ and $\gamma_{2}$ are used to control the curvature and asymmetry of the invariant (see Fig 3). The shape of the invariant is directly reflected in the price curve (see Fig. 5). Thus, the parameters of the model must be chosen in such a way as to make the price curve attractive for DEX users. Therefore, they must maintain a low price impact and protect against possible imbalances in the liquidity pool.

In this section we present a practical method for deciding the value of parameters. Since it is not very practical to consider a curve that depends on the quantity of tokens in the liquidity pool, we define the fraction of token $X$ inside the pool as $\frac{x}{x+py}$ and similarly we define the fraction of token $Y$ as $\frac{py}{x+py}$. As usual we multiply $y$ by $p$ to have the same units as $x$. It is very convenient to express the price curve as a function of the fraction, because the fraction is a percentage and does not depend on the size of the pool. Now we can choose a level of price impact that suits us, and choose the parameters in such a way that the price is impacted by this price impact at a certain level of fraction of the tokens.
In Fig. 9, we show examples of price curves for different parameter values. We choose a 5% level of price impact for $Y$. We are considering the case $x>py$ such that we are able to calculate $\gamma_{2}$. The same procedure should also be done for $x<py$ and considering a price impact for $X$ in order to find $\gamma_{1}$.
In this example, the price of the token $Y$ (SILK) at equilibrium is $p=1$. Its value after the 5% price impact increase is $1.05$ and it is represented by the points on the price curves in the figure. We can see how by changing $A$ and $\gamma_{2}$ it is possible to move the point to the right or to the left as desired.
A look at the curvature can also be helpful in choosing these parameters.