Miller’s Algorithm
Consider $E[l]$, the $l$torsion points of an elliptic curve $E$. Write $f_P$ for the rational function with divisor $l(P)l(O)$.
Let $g_{P,Q}$ be the line between the points $P$ and $Q$, and let $f_c$ be the function with divisor $(f_c) = c(P)  (c P)  (c1)(O)$. Then for all $a,b \in \mathbb{Z}$, we have $f_{a+b}(Q) = f_a(Q) f_b(Q) g_{aP,bP}(Q) / g_{(a+b)P, (a+b)P}(Q)$. Let the binary representation of $l$ be $l_t ,..., l_0$. Then Miller’s algorithm computes $f_P(Q)$ as follows.

set $f \leftarrow 1$ and $V \leftarrow P$

for $i \leftarrow t1$ to 0 do

set $f \leftarrow f^2 g_{V,V}(Q)/g_{2V, 2V}(Q)$ and $V \leftarrow 2V$

if $l_i = 1$ then
set $f \leftarrow f g_{V,P}(Q)/g_{V+P, (V+P)}(Q)$ and $V \leftarrow V + P$

At the end, $f = f_l(Q) = f_P(Q)$, and $V = l P = O$.
By adding extra logic in the above algorithm, we can avoid handling points of infinity when computing the $g$ functions. On the last iteration, we know $V = P$ if $l$ is odd, and $V$ has order 2 if $l$ is even.
The $g_{P,Q}$ functions
Let the curve be $Y = X^3 + a X + b$.

Tangents: At the point $(x,y)$, the line describing the tangent at that point is $\lambda X + Y + (y + \lambda x)$, where $\lambda = \frac{3x^2 + a}{2y}$.

Vertical lines: These are lines between $P$ and $P$. Let $P=(x,y)$. Then the vertical line through $P$ is $X + (x)$.

Other lines: The line between $P=(x_1, y_1)$ and $Q=(x_2, y_2)$ is given by $  \lambda X + Y + (y_1 + \lambda x_1)$ where $\lambda = \frac{y_2  y_1}{x_2  x_1}$.
It may be more efficient scale $\lambda$ appropriately and do all the divisions at once. If the embedding degree is at least two, and $P, Q$ are in the base field, the Tate exponentiation will get rid of any scaling factors. If the embedding degree is one, then since one has to compute $f(Q+R)/f(Q)$ (which should be computed with one iteration of Miller’s algorithm), the scaling factors in the numerator cancel those in the denominator. Thus the functions may be scaled as follows:

Tangents:

Verticals:

Lines:
(TODO: move to optimizations page since it assumes we’re using Tate. mention projective coords and mixin) If $P$ lies in the ground field and $k \gt 2$, then many divisions can be avoided since we may scale by quantities such as $2y$ and $x_2  x_1$ in the above equations. These quantities lie in the ground field hence become 1 during exponentiation.
Note
that is, $f_l = f_{l+1}$, hence if $l = 2^s  1$, we can easily compute $f_l(Q)$ using basically $s$ squarings. For Solinas primes of the form $l = 2^a + 2^b  1$ (where $a \gt b$), computing $f_{2^b}(Q)$ is an intermediate step in computing $f_{2^a}(Q)$, and then we may simply use one of the above formulas to compute $f_{2^a + 2^b}(Q) = f_{l+1}(Q) = f_l(Q)$.
Sometimes we need to compute $f_{ab}$. First note
hence
Thus
Under certain conditions(described later), then we have $f_{l1} = f_l$ so the last iteration of Miller’s algorithm can be simplified.