How to write any function in Solidity

11 February 2020 (updated 7 March 2024)

Rational approximations

$$ \frac{P(x_i)}{Q(x_i)} = f(x_i) $$

$$ P(x) = p_0 + p_1 x + p_2 x^2 + \dots $$

$$ Q(x) = 1 + q_1 x + q_2 x^2 + \dots $$

$$ P(x_i) = Q(x_i) f(x_i) $$

$$ \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots \\ 1 & x_1 & x_1^2 & \cdots \\ 1 & x_2 & x_2^2 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix} p_0 \\ p_1 \\ p_2 \\ \vdots \end{bmatrix} = \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots \\ 1 & x_1 & x_1^2 & \cdots \\ 1 & x_2 & x_2^2 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix} 1 \\ q_1 \\ q_2 \\ \vdots \end{bmatrix} \odot \begin{bmatrix} f(x_0) \\ f(x_1) \\ f(x_2) \\ \vdots \end{bmatrix} $$

$$ \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots \\ 1 & x_1 & x_1^2 & \cdots \\ 1 & x_2 & x_2^2 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix} p_0 \\ p_1 \\ p_2 \\ \vdots \end{bmatrix} = \begin{bmatrix} f(x_0) & f(x_0) x_0 & f(x_0) x_0^2 & \cdots \\ f(x_1) & f(x_1)x_1 & f(x_1)x_1^2 & \cdots \\ f(x_2) & f(x_2)x_2 & f(x_2)x_2^2 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix} 1 \\ q_1 \\ q_2 \\ \vdots \end{bmatrix} $$

$$ \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots \\ 1 & x_1 & x_1^2 & \cdots \\ 1 & x_2 & x_2^2 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix} p_0 \\ p_1 \\ p_2 \\ \vdots \end{bmatrix} = \begin{bmatrix} f(x_0) x_0 & f(x_0) x_0^2 & \cdots \\ f(x_1)x_1 & f(x_1)x_1^2 & \cdots \\ f(x_2)x_2 & f(x_2)x_2^2 & \cdots \\ \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix} q_1 \\ q_2 \\ \vdots \end{bmatrix} + \begin{bmatrix} f(x_0) \\ f(x_1) \\ f(x_2) \\ \vdots \end{bmatrix} $$

$$ \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & -f(x_0)x_0 & -f(x_0)x_0^2 & \cdots \\ 1 & x_1 & x_1^2 & \cdots & -f(x_1)x_1 & -f(x_1)x_1^2 & \cdots \\ 1 & x_2 & x_2^2 & \cdots & -f(x_2)x_2 & -f(x_2)x_2^2 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix} p_0 \\ p_1 \\ p_2 \\ \vdots \\ q_1 \\ q_2 \\ \vdots \end{bmatrix} = \begin{bmatrix} f(x_0) \\ f(x_1) \\ f(x_2) \\ \vdots \end{bmatrix} $$

Remez algorithm

By the Equioscilation theorem we know that when $P$ and $Q$ are optimal, the errors is oscilating between $\pm E$. Given the points of extremal error, $x_i$, this means

$$ \frac{P(x_i)}{Q(x_i)} - f(x_i) = (-1)^i E $$

This formula is the basis of the Remez algorithm. We itteratively solve for this equation, and use the result to find better approximations of the extremal points $x_i$.

$$ \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & -f(x_0)x_0 & -f(x_0)x_0^2 & \cdots & -1 \\ 1 & x_1 & x_1^2 & \cdots & -f(x_1)x_1 & -f(x_1)x_1^2 & \cdots & +1 \\ 1 & x_2 & x_2^2 & \cdots & -f(x_2)x_2 & -f(x_2)x_2^2 & \cdots & -1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \end{bmatrix} \begin{bmatrix} p_0 \\ p_1 \\ p_2 \\ \vdots \\ q_1 \\ q_2 \\ \vdots \\ E \end{bmatrix} = \begin{bmatrix} f(x_0) \\ f(x_1) \\ f(x_2) \\ \vdots \end{bmatrix} $$

TODO. In the complex domain, the ideal errors seem to be $e^{\pi i \frac{i}{n}} E$, this suggests that in the real domain we can use $\cos \frac{i}{n} E$. In a way, this is a generalization of the 'Modified Remez Algorithm' that also adds zero-error points $\frac{P(x_i)}{Q(x_i)} - f(x_i) = 0$.

Clenshaw algorithm and generalization of basis

We started out with polynomials of the form

$$ P(x) = p_0 + p_1 x + p_2 x^2 + \dots $$

If $\vert x \vert > 1$, the $x^i$ terms will grow exponentially, which can lead to numericall instabilities. A simple trick to avoid this is to apply a linear function to the domain to map it into $[-1, 1]$.

A more sophisticated way of avoiding it is by changing our basis polynomials. Let's say that till now we have used the basis

$$ \phi_i(x) = x^i $$

such that

$$ P(x) = p_0 \phi_0(x) + p_1 \phi_1(x) + p_2 \phi_2(x) + \dots $$

If these basis functions satisfy a recurrence

$$ \phi_{i+1}(x) = \alpha_i(x) \phi_i(x) + \beta_i(x) \phi_{i-1}(x) $$

then a generalization of Horner's evaluation mechanism applies:

The regular basis $\phi_i(x) = x^i$ satisfies with $\alpha_i(x) = x$ and $\beta_i(x) = 0$.

The Chebyshev basis satisfies $\alpha_i(x) = 2x$ and $\beta_i(x) = -1$.

$$ r_k = p_k + \alpha_k(x) r_{k+1} + \beta_{k+1} r_{k+2} $$

Assuming $\alpha$ and $\beta$ do not depend on the index:

Horner case using one variable

a = p_k + x * a

Chebyshev case using two variables

a = p_4 + 2x * b - a
b = p_3 + 2x * a - b
a = p_2 + 2x * b - a
b = p_1 + 2x * a - b

Chebyshev uses one more variable and has an additional subtraction, but same number of multiplications.

TODO Barycentric basis

https://web.archive.org/web/20200211232711/https://www.mn.uio.no/math/english/people/aca/michaelf/papers/rational.pdf