Pairing Programming
Needs maintenance; refer to my thesis instead.
Minimal pairingbased cryptography requires:

Arithmetic in \(\mathbb{Z}_p\): I built mine on top of the GMP library, which conveniently provides number theoretic functions such as inversion modulo a prime and the Jacobi symbol.

Elliptic curve groups: mostly routines for solving \(y^2 = x^3 + a x + b\) over \(\mathbb{Z}_p\), point addition and multiplication.

Bilinear pairing: Miller’s algorithm.
The above suffices for what I call Type A and B (and A1) pairings: very fast, but elements take a lot of space to represent.
Field extensions must be implemented for other pairing types.
Finding certain pairingfriendly curves requires more work. The the MNT curve construction method requires routines for finding roots modulo a given prime, testing polynomial irreducibility, computing Hilbert polynomials. These in turn depend on high precision complex floating point arithmetic and also an algorithm to solve a Pelltype equation.
Arithmetic on Elliptic Curves
Curves over fields of characteristic greater than 3 can be written in the form \(y^2 = x^3 + a x + b\). For such curves, point addition and multiplication routines are straightforward. (See Blake, Seroussi and Smart.) We give the explicit formulae here.
Suppose we wish to compute \(R = P + Q\). If \(P = O\), then \(R = Q\) and similarly \(Q = O\) is easily handled. If \(P_x = Q_x\) then either \(P_y = Q_y\) in which case \(R = O\) (for we have \(P = Q\)), or \(P = Q\) in which case the point doubling algorithm is invoked. Note that the point doubling algorithm assumes the input does not have order 2, i.e. the \(y\)coordinate is nonzero. (If \(P\) has order 2, the result is \(O\).)
Point addition for \(P \ne Q\):
Point doubling:
Generating Random Points
The standard method to generate a random point on an elliptic curve is to choose a random \(x\)coordinate and solve a quadratic equation for \(y\). (If no solution exists, a new \(x\)coordinate is chosen.) For odd characteristics, this can be done once one is able to find square roots of elements. In a finite field of prime order, square roots can be obtained via the TonelliShanks algorithm. In general fields, see D.J. Bernstein’s draft of a paper entitled "Faster square roots in annoying finite fields". At the moment I just solve \(x^2  a\) using the CantorZassenhaus method (which is referred to as Legendre’s method in the aforementioned reference).
Easier ways to generate random points exist for particular cases. For example if \(p\) is a prime with \(p = 2 \bmod 3\) then the cube root of \(x \in \mathbb{F}_q\) is simply \(x^{(2p1)/3}\). Thus if \(a = 0\) in the equation of the elliptic curve, then it is a simple matter to pick \(y\) first and then solve for \(x\).