diff --git a/zksnark/bn128_curve.py b/zksnark/bn128_curve.py
index 8827193..9b8940f 100644
--- a/zksnark/bn128_curve.py
+++ b/zksnark/bn128_curve.py
@@ -7,7 +7,6 @@ b = FQ(3)
 b2 = FQ2([3, 0])
 b12 = FQ12([3] + [0] * 11) / FQ12([0] * 6 + [1] + [0] * 5)
 
-ate_loop_count = 29793968203157093288
 
 G1 = (FQ(1), FQ(2))
 # Second element corresponds to modsqrt(67) * i in our quadratic field representation
diff --git a/zksnark/bn128_pairing.py b/zksnark/bn128_pairing.py
new file mode 100644
index 0000000..41febd3
--- /dev/null
+++ b/zksnark/bn128_pairing.py
@@ -0,0 +1,64 @@
+# NOT YET FINISHED! Pairing code has bugs in it, is NOT bilinear!
+
+from bn128_curve import double, add, multiply, is_on_curve, twist, b, b2, b12, curve_order, G1, G2
+from bn128_field_elements import field_modulus, FQ, FQ2, FQ12
+
+ate_loop_count = 29793968203157093288
+log_ate_loop_count = 64
+
+# Create a function representing the line between P1 and P2,
+# and evaluate it at T
+def linefunc(P1, P2, T):
+    x1, y1 = P1
+    x2, y2 = P2
+    xt, yt = T
+    if x1 != x2:
+        m = (y2 - y1) / (x2 - x1)
+        return (yt - y1) - m * (xt - x1)
+    elif y1 == y2:
+        m = 3 * x1**2 / (2 * y1)
+        return (yt - y1) - m * (xt - x1)
+    else:
+        return xt - x1
+
+def cast_point_to_fq12(pt):
+    if pt is None:
+        return None
+    x, y = pt
+    return (FQ12([x.n] + [0] * 11), FQ12([y.n] + [0] * 11))
+
+# Check consistency of the "line function"
+one, two, three, negone = G1, double(G1), multiply(G1, 3), multiply(G1, curve_order - 1)
+
+assert linefunc(one, two, one) == FQ(0)
+assert linefunc(one, two, two) == FQ(0)
+assert linefunc(one, two, three) != FQ(0)
+assert linefunc(one, negone, one) == FQ(0)
+assert linefunc(one, negone, negone) == FQ(0)
+assert linefunc(one, negone, two) != FQ(0)
+assert linefunc(one, one, one) == FQ(0)
+assert linefunc(one, one, two) != FQ(0)
+
+# Main miller loop
+def miller_loop(Q, P):
+    R = Q
+    f = FQ12.one()
+    for i in range(log_ate_loop_count, -1, -1):
+        f = f * f / linefunc(R, R, P)
+        R = double(R)
+        if ate_loop_count & (2**i):
+            f = f * linefunc(R, Q, P)
+            R = add(R, Q)
+    Q1 = (Q[0] ** field_modulus, Q[1] ** field_modulus)
+    nQ2 = (Q[0] ** (field_modulus ** 2), -Q[1] ** (field_modulus ** 2))
+    f = f * linefunc(R, Q1, P)
+    R = add(R, Q1)
+    f = f * linefunc(R, nQ2, P)
+    R = add(R, nQ2)
+    return f ** ((field_modulus ** 12 - 1) / curve_order)
+
+# Pairing computation
+def pairing(Q, P):
+    assert is_on_curve(Q, b2)
+    assert is_on_curve(P, b)
+    return miller_loop(twist(Q), cast_point_to_fq12(P))