# How to write any function in Solidity

## 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

Remco Bloemen
Math & Engineering
https://2π.com