diff --git a/src/blst_util.c b/src/blst_util.c
index b1b6ed8..04c5fa9 100644
--- a/src/blst_util.c
+++ b/src/blst_util.c
@@ -52,6 +52,19 @@ void p1_sub(blst_p1 *out, const blst_p1 *a, const blst_p1 *b) {
     blst_p1_add_or_double(out, a, &bneg);
 }
 
+// TODO: Is there really no better way to do this?
+void p2_mul(blst_p2 *out, const blst_p2 *a, const blst_fr *b) {
+    blst_scalar s;
+    blst_scalar_from_fr(&s, b);
+    blst_p2_mult(out, a, s.b, 8 * sizeof(blst_scalar));
+}
+
+void p2_sub(blst_p2 *out, const blst_p2 *a, const blst_p2 *b) {
+    blst_p2 bneg = *b;
+    blst_p2_cneg(&bneg, true);
+    blst_p2_add_or_double(out, a, &bneg);
+}
+
 // TODO: would be good to have an optimised multiexp for this
 void linear_combination_g1(blst_p1 *out, const blst_p1 *p, const blst_fr *coeffs, const uint64_t len) {
     blst_p1 tmp;
@@ -61,3 +74,27 @@ void linear_combination_g1(blst_p1 *out, const blst_p1 *p, const blst_fr *coeffs
         blst_p1_add_or_double(out, out, &tmp);
     }
 }
+
+bool pairings_verify(const blst_p1 *aa1, const blst_p2 *aa2, const blst_p1 *bb1, const blst_p2 *bb2) {
+    blst_fp12 loop0, loop1, gt_point;
+    blst_p1_affine a1, b1;
+    blst_p2_affine a2, b2;
+
+    // As an optimisation, we want to invert one of the pairings,
+    // so we negate one of the points.
+    blst_p1 a1neg = *aa1;
+    blst_p1_cneg(&a1neg, true);
+
+    blst_p1_to_affine(&a1, &a1neg);
+    blst_p1_to_affine(&b1, bb1);
+    blst_p2_to_affine(&a2, aa2);
+    blst_p2_to_affine(&b2, bb2);
+
+    blst_miller_loop(&loop0, &a2, &a1);
+    blst_miller_loop(&loop1, &b2, &b1);
+
+    blst_fp12_mul(&gt_point, &loop0, &loop1);
+    blst_final_exp(&gt_point, &gt_point);
+
+    return blst_fp12_is_one(&gt_point);
+}
diff --git a/src/blst_util.h b/src/blst_util.h
index 2720573..a8a5641 100644
--- a/src/blst_util.h
+++ b/src/blst_util.h
@@ -29,4 +29,7 @@ bool fr_equal(const blst_fr *aa, const blst_fr *bb);
 void fr_negate(blst_fr *out, const blst_fr *in);
 void p1_mul(blst_p1 *out, const blst_p1 *a, const blst_fr *b);
 void p1_sub(blst_p1 *out, const blst_p1 *a, const blst_p1 *b);
+void p2_mul(blst_p2 *out, const blst_p2 *a, const blst_fr *b);
+void p2_sub(blst_p2 *out, const blst_p2 *a, const blst_p2 *b);
 void linear_combination_g1(blst_p1 *out, const blst_p1 *p, const blst_fr *coeffs, const uint64_t len);
+bool pairings_verify(const blst_p1 *a1, const blst_p2 *a2, const blst_p1 *b1, const blst_p2 *b2);
diff --git a/src/blst_util_test.c b/src/blst_util_test.c
index 1e7d2ed..4926405 100644
--- a/src/blst_util_test.c
+++ b/src/blst_util_test.c
@@ -78,6 +78,35 @@ void p1_sub_works(void) {
     TEST_CHECK(blst_p1_is_equal(&tmp, &res));
 }
 
+void p2_mul_works(void) {
+    blst_fr minus1;
+    blst_p2 g2_gen, g2_gen_neg, res;
+
+    // Multiply the generator by minus one (the second root of unity)
+    blst_p2_from_affine(&g2_gen, &BLS12_381_G2);
+    blst_fr_from_uint64(&minus1, m1);
+    p2_mul(&res, &g2_gen, &minus1);
+
+    // We should end up with negative the generator
+    blst_p2_from_affine(&g2_gen_neg, &BLS12_381_NEG_G2);
+
+    TEST_CHECK(blst_p2_is_equal(&res, &g2_gen_neg));
+}
+
+void p2_sub_works(void) {
+    blst_p2 g2_gen, g2_gen_neg;
+    blst_p2 tmp, res;
+
+    blst_p2_from_affine(&g2_gen, &BLS12_381_G2);
+    blst_p2_from_affine(&g2_gen_neg, &BLS12_381_NEG_G2);
+
+    // 2 * g2_gen = g2_gen - g2_gen_neg
+    blst_p2_double(&tmp, &g2_gen);
+    p2_sub(&res, &g2_gen, &g2_gen_neg);
+
+    TEST_CHECK(blst_p2_is_equal(&tmp, &res));
+}
+
 void identity_g1_is_infinity(void) {
     blst_p1 identity_g1;
     blst_p1_from_affine(&identity_g1, &identity_g1_affine);
@@ -103,6 +132,25 @@ void g1_linear_combination(void) {
     TEST_CHECK(blst_p1_is_equal(&exp, &res));
 }
 
+void pairings_work(void) {
+    // Verify that e([3]g1, [5]g2) = e([5]g1, [3]g2)
+    blst_fr three, five;
+    blst_p1 g1_3, g1_5;
+    blst_p2 g2_3, g2_5;
+
+    // Set up
+    fr_from_uint64(&three, 3);
+    fr_from_uint64(&five, 5);
+    p1_mul(&g1_3, blst_p1_generator(), &three);
+    p1_mul(&g1_5, blst_p1_generator(), &five);
+    p2_mul(&g2_3, blst_p2_generator(), &three);
+    p2_mul(&g2_5, blst_p2_generator(), &five);
+
+    // Verify the pairing
+    TEST_CHECK(true == pairings_verify(&g1_3, &g2_5, &g1_5, &g2_3));
+    TEST_CHECK(false == pairings_verify(&g1_3, &g2_3, &g1_5, &g2_5));
+}
+
 TEST_LIST =
     {
      {"fr_is_one_works", fr_is_one_works },
@@ -111,7 +159,10 @@ TEST_LIST =
      {"fr_negate_works", fr_negate_works},
      {"p1_mul_works", p1_mul_works},
      {"p1_sub_works", p1_sub_works},
+     {"p2_mul_works", p2_mul_works},
+     {"p2_sub_works", p2_sub_works},
      {"identity_g1_is_infinity", identity_g1_is_infinity},
      {"g1_linear_combination", g1_linear_combination},
+     {"pairings_work", pairings_work},
      { NULL, NULL }     /* zero record marks the end of the list */
     };