Initial commit

src/bls12_381/README.md

@@ -0,0 +1,57 @@
+# BLS12-381
+
+This is an implementation of the BLS12-381 pairing-friendly elliptic curve construction.
+
+## BLS12 Parameterization
+
+BLS12 curves are parameterized by a value *x* such that the base field modulus *q* and subgroup *r* can be computed by:
+
+* q = (x - 1)<sup>2</sup> ((x<sup>4</sup> - x<sup>2</sup> + 1) / 3) + x
+* r = (x<sup>4</sup> - x<sup>2</sup> + 1)
+
+Given primes *q* and *r* parameterized as above, we can easily construct an elliptic curve over the prime field F<sub>*q*</sub> which contains a subgroup of order *r* such that *r* | (*q*<sup>12</sup> - 1), giving it an embedding degree of 12. Instantiating its sextic twist over an extension field F<sub>q<sup>2</sup></sub> gives rise to an efficient bilinear pairing function between elements of the order *r* subgroups of either curves, into an order *r* multiplicative subgroup of F<sub>q<sup>12</sup></sub>.
+
+In zk-SNARK schemes, we require F<sub>r</sub> with large 2<sup>n</sup> roots of unity for performing efficient fast-fourier transforms. As such, guaranteeing that large 2<sup>n</sup> | (r - 1), or equivalently that *x* has a large 2<sup>n</sup> factor, gives rise to BLS12 curves suitable for zk-SNARKs.
+
+Due to recent research, it is estimated by many that *q* should be approximately 384 bits to target 128-bit security. Conveniently, *r* is approximately 256 bits when *q* is approximately 384 bits, making BLS12 curves ideal for 128-bit security. It also makes them ideal for many zk-SNARK applications, as the scalar field can be used for keying material such as embedded curve constructions.
+
+Many curves match our descriptions, but we require some extra properties for efficiency purposes:
+
+* *q* should be smaller than 2<sup>383</sup>, and *r* should be smaller than 2<sup>255</sup>, so that the most significant bit is unset when using 64-bit or 32-bit limbs. This allows for cheap reductions.
+* F<sub>q<sup>12</sup></sub> is typically constructed using towers of extension fields. As a byproduct of [research](https://eprint.iacr.org/2011/465.pdf) for BLS curves of embedding degree 24, we can identify subfamilies of BLS12 curves (for our purposes, where x mod 72 = {16, 64}) that produce efficient extension field towers and twisting isomorphisms.
+* We desire *x* of small Hamming weight, to increase the performance of the pairing function.
+
+## BLS12-381 Instantiation
+
+The BLS12-381 construction is instantiated by `x = -0xd201000000010000`, which produces the largest `q` and smallest Hamming weight of `x` that meets the above requirements. This produces:
+
+* q = `0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaab` (381 bits)
+* r = `0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001` (255 bits)
+
+Our extension field tower is constructed as follows:
+
+1. F<sub>q<sup>2</sup></sub> is constructed as F<sub>q</sub>(u) / (u<sup>2</sup> - β) where β = -1.
+2. F<sub>q<sup>6</sup></sub> is constructed as F<sub>q<sup>2</sup></sub>(v) / (v<sup>3</sup> - ξ) where ξ = u + 1
+3. F<sub>q<sup>12</sup></sub> is constructed as F<sub>q<sup>6</sup></sub>(w) / (w<sup>2</sup> - γ) where γ = v
+
+Now, we instantiate the elliptic curve E(F<sub>q</sub>) : y<sup>2</sup> = x<sup>3</sup> + 4, and the elliptic curve E'(F<sub>q<sup>2</sup></sub>) : y<sup>2</sup> = x<sup>3</sup> + 4v.
+
+The group G<sub>1</sub> is the *r* order subgroup of E, which has cofactor (x - 1)<sup>2</sup> / 3. The group G<sub>2</sub> is the *r* order subgroup of E', which has cofactor (x<sup>8</sup> - 4x<sup>7</sup> + 5x<sup>6</sup> - 4x<sup>4</sup> + 6x<sup>3</sup> - 4x<sup>2</sup> - 4x + 13) / 9.
+
+### Generators
+
+The generators of G<sub>1</sub> and G<sub>2</sub> are computed by finding the lexicographically smallest valid `x`-coordinate, and its lexicographically smallest `y`-coordinate and scaling it by the cofactor such that the result is not the point at infinity.
+
+#### G1
+
+```
+x = 3685416753713387016781088315183077757961620795782546409894578378688607592378376318836054947676345821548104185464507
+y = 1339506544944476473020471379941921221584933875938349620426543736416511423956333506472724655353366534992391756441569
+```
+
+#### G2
+
+```
+x = 3059144344244213709971259814753781636986470325476647558659373206291635324768958432433509563104347017837885763365758*u + 352701069587466618187139116011060144890029952792775240219908644239793785735715026873347600343865175952761926303160
+y = 927553665492332455747201965776037880757740193453592970025027978793976877002675564980949289727957565575433344219582*u + 1985150602287291935568054521177171638300868978215655730859378665066344726373823718423869104263333984641494340347905
+```

src/bls12_381/fq12.rs

@@ -0,0 +1,289 @@
+use rand::{Rng, Rand};
+use ::{Field};
+use super::fq6::Fq6;
+use super::fq2::Fq2;
+use super::fq::{FROBENIUS_COEFF_FQ12_C1};
+
+/// An element of F_{q^12}, represented by c0 + c1 * w.
+#[derive(Copy, Clone, Debug, Eq, PartialEq)]
+pub struct Fq12 {
+    pub c0: Fq6,
+    pub c1: Fq6
+}
+
+impl Rand for Fq12 {
+    fn rand<R: Rng>(rng: &mut R) -> Self {
+        Fq12 {
+            c0: rng.gen(),
+            c1: rng.gen()
+        }
+    }
+}
+
+impl Fq12 {
+    pub fn unitary_inverse(&mut self)
+    {
+        self.c1.negate();
+    }
+
+    pub fn mul_by_014(
+        &mut self,
+        c0: &Fq2,
+        c1: &Fq2,
+        c4: &Fq2
+    )
+    {
+        let mut aa = self.c0;
+        aa.mul_by_01(c0, c1);
+        let mut bb = self.c1;
+        bb.mul_by_1(c4);
+        let mut o = *c1;
+        o.add_assign(c4);
+        self.c1.add_assign(&self.c0);
+        self.c1.mul_by_01(c0, &o);
+        self.c1.sub_assign(&aa);
+        self.c1.sub_assign(&bb);
+        self.c0 = bb;
+        self.c0.mul_by_nonresidue();
+        self.c0.add_assign(&aa);
+    }
+}
+
+impl Field for Fq12
+{
+    fn zero() -> Self {
+        Fq12 {
+            c0: Fq6::zero(),
+            c1: Fq6::zero()
+        }
+    }
+
+    fn one() -> Self {
+        Fq12 {
+            c0: Fq6::one(),
+            c1: Fq6::zero()
+        }
+    }
+
+    fn is_zero(&self) -> bool {
+        self.c0.is_zero() && self.c1.is_zero()
+    }
+
+    fn double(&mut self) {
+        self.c0.double();
+        self.c1.double();
+    }
+
+    fn negate(&mut self) {
+        self.c0.negate();
+        self.c1.negate();
+    }
+
+    fn add_assign(&mut self, other: &Self) {
+        self.c0.add_assign(&other.c0);
+        self.c1.add_assign(&other.c1);
+    }
+
+    fn sub_assign(&mut self, other: &Self) {
+        self.c0.sub_assign(&other.c0);
+        self.c1.sub_assign(&other.c1);
+    }
+
+    fn frobenius_map(&mut self, power: usize)
+    {
+        self.c0.frobenius_map(power);
+        self.c1.frobenius_map(power);
+
+        self.c1.c0.mul_assign(&FROBENIUS_COEFF_FQ12_C1[power % 12]);
+        self.c1.c1.mul_assign(&FROBENIUS_COEFF_FQ12_C1[power % 12]);
+        self.c1.c2.mul_assign(&FROBENIUS_COEFF_FQ12_C1[power % 12]);
+    }
+
+    fn square(&mut self) {
+        let mut ab = self.c0;
+        ab.mul_assign(&self.c1);
+        let mut c0c1 = self.c0;
+        c0c1.add_assign(&self.c1);
+        let mut c0 = self.c1;
+        c0.mul_by_nonresidue();
+        c0.add_assign(&self.c0);
+        c0.mul_assign(&c0c1);
+        c0.sub_assign(&ab);
+        self.c1 = ab;
+        self.c1.add_assign(&ab);
+        ab.mul_by_nonresidue();
+        c0.sub_assign(&ab);
+        self.c0 = c0;
+    }
+
+    fn mul_assign(&mut self, other: &Self) {
+        let mut aa = self.c0;
+        aa.mul_assign(&other.c0);
+        let mut bb = self.c1;
+        bb.mul_assign(&other.c1);
+        let mut o = other.c0;
+        o.add_assign(&other.c1);
+        self.c1.add_assign(&self.c0);
+        self.c1.mul_assign(&o);
+        self.c1.sub_assign(&aa);
+        self.c1.sub_assign(&bb);
+        self.c0 = bb;
+        self.c0.mul_by_nonresidue();
+        self.c0.add_assign(&aa);
+    }
+
+    fn inverse(&self) -> Option<Self> {
+        let mut c0s = self.c0;
+        c0s.square();
+        let mut c1s = self.c1;
+        c1s.square();
+        c1s.mul_by_nonresidue();
+        c0s.sub_assign(&c1s);
+
+        c0s.inverse().map(|t| {
+            let mut tmp = Fq12 {
+                c0: t,
+                c1: t
+            };
+            tmp.c0.mul_assign(&self.c0);
+            tmp.c1.mul_assign(&self.c1);
+            tmp.c1.negate();
+
+            tmp
+        })
+    }
+}
+
+#[cfg(test)]
+use rand::{SeedableRng, XorShiftRng};
+
+#[test]
+fn test_fq12_mul_by_014() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..1000 {
+        let c0 = Fq2::rand(&mut rng);
+        let c1 = Fq2::rand(&mut rng);
+        let c5 = Fq2::rand(&mut rng);
+        let mut a = Fq12::rand(&mut rng);
+        let mut b = a;
+
+        a.mul_by_014(&c0, &c1, &c5);
+        b.mul_assign(&Fq12 {
+            c0: Fq6 {
+                c0: c0,
+                c1: c1,
+                c2: Fq2::zero()
+            },
+            c1: Fq6 {
+                c0: Fq2::zero(),
+                c1: c5,
+                c2: Fq2::zero()
+            }
+        });
+
+        assert_eq!(a, b);
+    }
+}
+
+#[bench]
+fn bench_fq12_add_assign(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<(Fq12, Fq12)> = (0..SAMPLES).map(|_| {
+        (Fq12::rand(&mut rng), Fq12::rand(&mut rng))
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let mut tmp = v[count].0;
+        tmp.add_assign(&v[count].1);
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_fq12_sub_assign(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<(Fq12, Fq12)> = (0..SAMPLES).map(|_| {
+        (Fq12::rand(&mut rng), Fq12::rand(&mut rng))
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let mut tmp = v[count].0;
+        tmp.sub_assign(&v[count].1);
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_fq12_mul_assign(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<(Fq12, Fq12)> = (0..SAMPLES).map(|_| {
+        (Fq12::rand(&mut rng), Fq12::rand(&mut rng))
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let mut tmp = v[count].0;
+        tmp.mul_assign(&v[count].1);
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_fq12_squaring(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<Fq12> = (0..SAMPLES).map(|_| {
+        Fq12::rand(&mut rng)
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let mut tmp = v[count];
+        tmp.square();
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_fq12_inverse(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<Fq12> = (0..SAMPLES).map(|_| {
+        Fq12::rand(&mut rng)
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let tmp = v[count].inverse();
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[test]
+fn fq12_field_tests() {
+    use ::PrimeField;
+
+    ::tests::field::random_field_tests::<Fq12>();
+    ::tests::field::random_frobenius_tests::<Fq12, _>(super::fq::Fq::char(), 13);
+}

src/bls12_381/fq2.rs

@@ -0,0 +1,504 @@
+use rand::{Rng, Rand};
+use ::{Field, SqrtField};
+use super::fq::{Fq, FROBENIUS_COEFF_FQ2_C1, NEGATIVE_ONE};
+
+/// An element of F_{q^2}, represented by c0 + c1 * u.
+#[derive(Copy, Clone, Debug, Eq, PartialEq)]
+pub struct Fq2 {
+    pub c0: Fq,
+    pub c1: Fq
+}
+
+impl Fq2 {
+    /// Multiply this element by the cubic and quadratic nonresidue 1 + u.
+    pub fn mul_by_nonresidue(&mut self) {
+        let t0 = self.c0;
+        self.c0.sub_assign(&self.c1);
+        self.c1.add_assign(&t0);
+    }
+}
+
+impl Rand for Fq2 {
+    fn rand<R: Rng>(rng: &mut R) -> Self {
+        Fq2 {
+            c0: rng.gen(),
+            c1: rng.gen()
+        }
+    }
+}
+
+impl Field for Fq2 {
+    fn zero() -> Self {
+        Fq2 { c0: Fq::zero(), c1: Fq::zero() }
+    }
+
+    fn one() -> Self {
+        Fq2 { c0: Fq::one(), c1: Fq::zero() }
+    }
+
+    fn is_zero(&self) -> bool {
+        self.c0.is_zero() && self.c1.is_zero()
+    }
+
+    fn square(&mut self) {
+        let mut ab = self.c0;
+        ab.mul_assign(&self.c1);
+        let mut c0c1 = self.c0;
+        c0c1.add_assign(&self.c1);
+        let mut c0 = self.c1;
+        c0.negate();
+        c0.add_assign(&self.c0);
+        c0.mul_assign(&c0c1);
+        c0.sub_assign(&ab);
+        self.c1 = ab;
+        self.c1.add_assign(&ab);
+        c0.add_assign(&ab);
+        self.c0 = c0;
+    }
+
+    fn double(&mut self) {
+        self.c0.double();
+        self.c1.double();
+    }
+
+    fn negate(&mut self) {
+        self.c0.negate();
+        self.c1.negate();
+    }
+
+    fn add_assign(&mut self, other: &Self) {
+        self.c0.add_assign(&other.c0);
+        self.c1.add_assign(&other.c1);
+    }
+
+    fn sub_assign(&mut self, other: &Self) {
+        self.c0.sub_assign(&other.c0);
+        self.c1.sub_assign(&other.c1);
+    }
+
+    fn mul_assign(&mut self, other: &Self) {
+        let mut aa = self.c0;
+        aa.mul_assign(&other.c0);
+        let mut bb = self.c1;
+        bb.mul_assign(&other.c1);
+        let mut o = other.c0;
+        o.add_assign(&other.c1);
+        self.c1.add_assign(&self.c0);
+        self.c1.mul_assign(&o);
+        self.c1.sub_assign(&aa);
+        self.c1.sub_assign(&bb);
+        self.c0 = aa;
+        self.c0.sub_assign(&bb);
+    }
+
+    fn inverse(&self) -> Option<Self> {
+        let mut t1 = self.c1;
+        t1.square();
+        let mut t0 = self.c0;
+        t0.square();
+        t0.add_assign(&t1);
+        t0.inverse().map(|t| {
+            let mut tmp = Fq2 {
+                c0: self.c0,
+                c1: self.c1
+            };
+            tmp.c0.mul_assign(&t);
+            tmp.c1.mul_assign(&t);
+            tmp.c1.negate();
+
+            tmp
+        })
+    }
+
+    fn frobenius_map(&mut self, power: usize)
+    {
+        self.c1.mul_assign(&FROBENIUS_COEFF_FQ2_C1[power % 2]);
+    }
+}
+
+impl SqrtField for Fq2 {
+    fn sqrt(&self) -> Option<Self> {
+        // Algorithm 9, https://eprint.iacr.org/2012/685.pdf
+
+        if self.is_zero() {
+            return Some(Self::zero());
+        } else {
+            // a1 = self^((q - 3) / 4)
+            let mut a1 = self.pow([0xee7fbfffffffeaaa, 0x7aaffffac54ffff, 0xd9cc34a83dac3d89, 0xd91dd2e13ce144af, 0x92c6e9ed90d2eb35, 0x680447a8e5ff9a6]);
+            let mut alpha = a1;
+            alpha.square();
+            alpha.mul_assign(self);
+            let mut a0 = alpha;
+            a0.frobenius_map(1);
+            a0.mul_assign(&alpha);
+
+            let neg1 = Fq2 {
+                c0: NEGATIVE_ONE,
+                c1: Fq::zero()
+            };
+
+            if a0 == neg1 {
+                None
+            } else {
+                a1.mul_assign(self);
+
+                if alpha == neg1 {
+                    a1.mul_assign(&Fq2{c0: Fq::zero(), c1: Fq::one()});
+                } else {
+                    alpha.add_assign(&Fq2::one());
+                    // alpha = alpha^((q - 1) / 2)
+                    alpha = alpha.pow([0xdcff7fffffffd555, 0xf55ffff58a9ffff, 0xb39869507b587b12, 0xb23ba5c279c2895f, 0x258dd3db21a5d66b, 0xd0088f51cbff34d]);
+                    a1.mul_assign(&alpha);
+                }
+
+                Some(a1)
+            }
+        }
+    }
+}
+
+#[test]
+fn test_fq2_basics() {
+    assert_eq!(Fq2 { c0: Fq::zero(), c1: Fq::zero() }, Fq2::zero());
+    assert_eq!(Fq2 { c0: Fq::one(), c1: Fq::zero() }, Fq2::one());
+    assert!(Fq2::zero().is_zero());
+    assert!(!Fq2::one().is_zero());
+    assert!(!Fq2{c0: Fq::zero(), c1: Fq::one()}.is_zero());
+}
+
+#[test]
+fn test_fq2_squaring() {
+    use ::PrimeField;
+    use super::fq::{FqRepr};
+
+    let mut a = Fq2 { c0: Fq::one(), c1: Fq::one() }; // u + 1
+    a.square();
+    assert_eq!(a, Fq2 { c0: Fq::zero(), c1: Fq::from_repr(FqRepr::from(2)).unwrap() }); // 2u
+
+    let mut a = Fq2 { c0: Fq::zero(), c1: Fq::one() }; // u
+    a.square();
+    assert_eq!(a, {
+        let mut neg1 = Fq::one();
+        neg1.negate();
+        Fq2 { c0: neg1, c1: Fq::zero() }
+    }); // -1
+
+    let mut a = Fq2 {
+        c0: Fq::from_repr(FqRepr([0x9c2c6309bbf8b598, 0x4eef5c946536f602, 0x90e34aab6fb6a6bd, 0xf7f295a94e58ae7c, 0x41b76dcc1c3fbe5e, 0x7080c5fa1d8e042])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x38f473b3c870a4ab, 0x6ad3291177c8c7e5, 0xdac5a4c911a4353e, 0xbfb99020604137a0, 0xfc58a7b7be815407, 0x10d1615e75250a21])).unwrap()
+    };
+    a.square();
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0xf262c28c538bcf68, 0xb9f2a66eae1073ba, 0xdc46ab8fad67ae0, 0xcb674157618da176, 0x4cf17b5893c3d327, 0x7eac81369c43361])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0xc1579cf58e980cf8, 0xa23eb7e12dd54d98, 0xe75138bce4cec7aa, 0x38d0d7275a9689e1, 0x739c983042779a65, 0x1542a61c8a8db994])).unwrap()
+    });
+}
+
+#[test]
+fn test_fq2_mul() {
+    use ::PrimeField;
+    use super::fq::{FqRepr};
+
+    let mut a = Fq2 {
+        c0: Fq::from_repr(FqRepr([0x85c9f989e1461f03, 0xa2e33c333449a1d6, 0x41e461154a7354a3, 0x9ee53e7e84d7532e, 0x1c202d8ed97afb45, 0x51d3f9253e2516f])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0xa7348a8b511aedcf, 0x143c215d8176b319, 0x4cc48081c09b8903, 0x9533e4a9a5158be, 0x7a5e1ecb676d65f9, 0x180c3ee46656b008])).unwrap()
+    };
+    a.mul_assign(&Fq2 {
+        c0: Fq::from_repr(FqRepr([0xe21f9169805f537e, 0xfc87e62e179c285d, 0x27ece175be07a531, 0xcd460f9f0c23e430, 0x6c9110292bfa409, 0x2c93a72eb8af83e])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x4b1c3f936d8992d4, 0x1d2a72916dba4c8a, 0x8871c508658d1e5f, 0x57a06d3135a752ae, 0x634cd3c6c565096d, 0x19e17334d4e93558])).unwrap()
+    });
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0x95b5127e6360c7e4, 0xde29c31a19a6937e, 0xf61a96dacf5a39bc, 0x5511fe4d84ee5f78, 0x5310a202d92f9963, 0x1751afbe166e5399])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x84af0e1bd630117a, 0x6c63cd4da2c2aa7, 0x5ba6e5430e883d40, 0xc975106579c275ee, 0x33a9ac82ce4c5083, 0x1ef1a36c201589d])).unwrap()
+    });
+}
+
+#[test]
+fn test_fq2_inverse() {
+    use ::PrimeField;
+    use super::fq::{FqRepr};
+
+    assert!(Fq2::zero().inverse().is_none());
+
+    let a = Fq2 {
+        c0: Fq::from_repr(FqRepr([0x85c9f989e1461f03, 0xa2e33c333449a1d6, 0x41e461154a7354a3, 0x9ee53e7e84d7532e, 0x1c202d8ed97afb45, 0x51d3f9253e2516f])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0xa7348a8b511aedcf, 0x143c215d8176b319, 0x4cc48081c09b8903, 0x9533e4a9a5158be, 0x7a5e1ecb676d65f9, 0x180c3ee46656b008])).unwrap()
+    };
+    let a = a.inverse().unwrap();
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0x70300f9bcb9e594, 0xe5ecda5fdafddbb2, 0x64bef617d2915a8f, 0xdfba703293941c30, 0xa6c3d8f9586f2636, 0x1351ef01941b70c4])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x8c39fd76a8312cb4, 0x15d7b6b95defbff0, 0x947143f89faedee9, 0xcbf651a0f367afb2, 0xdf4e54f0d3ef15a6, 0x103bdf241afb0019])).unwrap()
+    });
+}
+
+#[test]
+fn test_fq2_addition() {
+    use ::PrimeField;
+    use super::fq::{FqRepr};
+
+    let mut a = Fq2 {
+        c0: Fq::from_repr(FqRepr([0x2d0078036923ffc7, 0x11e59ea221a3b6d2, 0x8b1a52e0a90f59ed, 0xb966ce3bc2108b13, 0xccc649c4b9532bf3, 0xf8d295b2ded9dc])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x977df6efcdaee0db, 0x946ae52d684fa7ed, 0xbe203411c66fb3a5, 0xb3f8afc0ee248cad, 0x4e464dea5bcfd41e, 0x12d1137b8a6a837])).unwrap()
+    };
+    a.add_assign(&Fq2 {
+        c0: Fq::from_repr(FqRepr([0x619a02d78dc70ef2, 0xb93adfc9119e33e8, 0x4bf0b99a9f0dca12, 0x3b88899a42a6318f, 0x986a4a62fa82a49d, 0x13ce433fa26027f5])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x66323bf80b58b9b9, 0xa1379b6facf6e596, 0x402aef1fb797e32f, 0x2236f55246d0d44d, 0x4c8c1800eb104566, 0x11d6e20e986c2085])).unwrap()
+    });
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0x8e9a7adaf6eb0eb9, 0xcb207e6b3341eaba, 0xd70b0c7b481d23ff, 0xf4ef57d604b6bca2, 0x65309427b3d5d090, 0x14c715d5553f01d2])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0xfdb032e7d9079a94, 0x35a2809d15468d83, 0xfe4b23317e0796d5, 0xd62fa51334f560fa, 0x9ad265eb46e01984, 0x1303f3465112c8bc])).unwrap()
+    });
+}
+
+#[test]
+fn test_fq2_subtraction() {
+    use ::PrimeField;
+    use super::fq::{FqRepr};
+
+    let mut a = Fq2 {
+        c0: Fq::from_repr(FqRepr([0x2d0078036923ffc7, 0x11e59ea221a3b6d2, 0x8b1a52e0a90f59ed, 0xb966ce3bc2108b13, 0xccc649c4b9532bf3, 0xf8d295b2ded9dc])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x977df6efcdaee0db, 0x946ae52d684fa7ed, 0xbe203411c66fb3a5, 0xb3f8afc0ee248cad, 0x4e464dea5bcfd41e, 0x12d1137b8a6a837])).unwrap()
+    };
+    a.sub_assign(&Fq2 {
+        c0: Fq::from_repr(FqRepr([0x619a02d78dc70ef2, 0xb93adfc9119e33e8, 0x4bf0b99a9f0dca12, 0x3b88899a42a6318f, 0x986a4a62fa82a49d, 0x13ce433fa26027f5])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x66323bf80b58b9b9, 0xa1379b6facf6e596, 0x402aef1fb797e32f, 0x2236f55246d0d44d, 0x4c8c1800eb104566, 0x11d6e20e986c2085])).unwrap()
+    });
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0x8565752bdb5c9b80, 0x7756bed7c15982e9, 0xa65a6be700b285fe, 0xe255902672ef6c43, 0x7f77a718021c342d, 0x72ba14049fe9881])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0xeb4abaf7c255d1cd, 0x11df49bc6cacc256, 0xe52617930588c69a, 0xf63905f39ad8cb1f, 0x4cd5dd9fb40b3b8f, 0x957411359ba6e4c])).unwrap()
+    });
+}
+
+#[test]
+fn test_fq2_negation() {
+    use ::PrimeField;
+    use super::fq::{FqRepr};
+
+    let mut a = Fq2 {
+        c0: Fq::from_repr(FqRepr([0x2d0078036923ffc7, 0x11e59ea221a3b6d2, 0x8b1a52e0a90f59ed, 0xb966ce3bc2108b13, 0xccc649c4b9532bf3, 0xf8d295b2ded9dc])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x977df6efcdaee0db, 0x946ae52d684fa7ed, 0xbe203411c66fb3a5, 0xb3f8afc0ee248cad, 0x4e464dea5bcfd41e, 0x12d1137b8a6a837])).unwrap()
+    };
+    a.negate();
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0x8cfe87fc96dbaae4, 0xcc6615c8fb0492d, 0xdc167fc04da19c37, 0xab107d49317487ab, 0x7e555df189f880e3, 0x19083f5486a10cbd])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x228109103250c9d0, 0x8a411ad149045812, 0xa9109e8f3041427e, 0xb07e9bc405608611, 0xfcd559cbe77bd8b8, 0x18d400b280d93e62])).unwrap()
+    });
+}
+
+#[test]
+fn test_fq2_doubling() {
+    use ::PrimeField;
+    use super::fq::{FqRepr};
+
+    let mut a = Fq2 {
+        c0: Fq::from_repr(FqRepr([0x2d0078036923ffc7, 0x11e59ea221a3b6d2, 0x8b1a52e0a90f59ed, 0xb966ce3bc2108b13, 0xccc649c4b9532bf3, 0xf8d295b2ded9dc])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x977df6efcdaee0db, 0x946ae52d684fa7ed, 0xbe203411c66fb3a5, 0xb3f8afc0ee248cad, 0x4e464dea5bcfd41e, 0x12d1137b8a6a837])).unwrap()
+    };
+    a.double();
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0x5a00f006d247ff8e, 0x23cb3d4443476da4, 0x1634a5c1521eb3da, 0x72cd9c7784211627, 0x998c938972a657e7, 0x1f1a52b65bdb3b9])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x2efbeddf9b5dc1b6, 0x28d5ca5ad09f4fdb, 0x7c4068238cdf674b, 0x67f15f81dc49195b, 0x9c8c9bd4b79fa83d, 0x25a226f714d506e])).unwrap()
+    });
+}
+
+#[test]
+fn test_fq2_frobenius_map() {
+    use ::PrimeField;
+    use super::fq::{FqRepr};
+
+    let mut a = Fq2 {
+        c0: Fq::from_repr(FqRepr([0x2d0078036923ffc7, 0x11e59ea221a3b6d2, 0x8b1a52e0a90f59ed, 0xb966ce3bc2108b13, 0xccc649c4b9532bf3, 0xf8d295b2ded9dc])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x977df6efcdaee0db, 0x946ae52d684fa7ed, 0xbe203411c66fb3a5, 0xb3f8afc0ee248cad, 0x4e464dea5bcfd41e, 0x12d1137b8a6a837])).unwrap()
+    };
+    a.frobenius_map(0);
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0x2d0078036923ffc7, 0x11e59ea221a3b6d2, 0x8b1a52e0a90f59ed, 0xb966ce3bc2108b13, 0xccc649c4b9532bf3, 0xf8d295b2ded9dc])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x977df6efcdaee0db, 0x946ae52d684fa7ed, 0xbe203411c66fb3a5, 0xb3f8afc0ee248cad, 0x4e464dea5bcfd41e, 0x12d1137b8a6a837])).unwrap()
+    });
+    a.frobenius_map(1);
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0x2d0078036923ffc7, 0x11e59ea221a3b6d2, 0x8b1a52e0a90f59ed, 0xb966ce3bc2108b13, 0xccc649c4b9532bf3, 0xf8d295b2ded9dc])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x228109103250c9d0, 0x8a411ad149045812, 0xa9109e8f3041427e, 0xb07e9bc405608611, 0xfcd559cbe77bd8b8, 0x18d400b280d93e62])).unwrap()
+    });
+    a.frobenius_map(1);
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0x2d0078036923ffc7, 0x11e59ea221a3b6d2, 0x8b1a52e0a90f59ed, 0xb966ce3bc2108b13, 0xccc649c4b9532bf3, 0xf8d295b2ded9dc])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x977df6efcdaee0db, 0x946ae52d684fa7ed, 0xbe203411c66fb3a5, 0xb3f8afc0ee248cad, 0x4e464dea5bcfd41e, 0x12d1137b8a6a837])).unwrap()
+    });
+    a.frobenius_map(2);
+    assert_eq!(a, Fq2 {
+        c0: Fq::from_repr(FqRepr([0x2d0078036923ffc7, 0x11e59ea221a3b6d2, 0x8b1a52e0a90f59ed, 0xb966ce3bc2108b13, 0xccc649c4b9532bf3, 0xf8d295b2ded9dc])).unwrap(),
+        c1: Fq::from_repr(FqRepr([0x977df6efcdaee0db, 0x946ae52d684fa7ed, 0xbe203411c66fb3a5, 0xb3f8afc0ee248cad, 0x4e464dea5bcfd41e, 0x12d1137b8a6a837])).unwrap()
+    });
+}
+
+#[test]
+fn test_fq2_sqrt() {
+    use ::PrimeField;
+    use super::fq::{FqRepr};
+
+    assert_eq!(
+        Fq2 {
+            c0: Fq::from_repr(FqRepr([0x476b4c309720e227, 0x34c2d04faffdab6, 0xa57e6fc1bab51fd9, 0xdb4a116b5bf74aa1, 0x1e58b2159dfe10e2, 0x7ca7da1f13606ac])).unwrap(),
+            c1: Fq::from_repr(FqRepr([0xfa8de88b7516d2c3, 0x371a75ed14f41629, 0x4cec2dca577a3eb6, 0x212611bca4e99121, 0x8ee5394d77afb3d, 0xec92336650e49d5])).unwrap()
+        }.sqrt().unwrap(),
+        Fq2 {
+            c0: Fq::from_repr(FqRepr([0x40b299b2704258c5, 0x6ef7de92e8c68b63, 0x6d2ddbe552203e82, 0x8d7f1f723d02c1d3, 0x881b3e01b611c070, 0x10f6963bbad2ebc5])).unwrap(),
+            c1: Fq::from_repr(FqRepr([0xc099534fc209e752, 0x7670594665676447, 0x28a20faed211efe7, 0x6b852aeaf2afcb1b, 0xa4c93b08105d71a9, 0x8d7cfff94216330])).unwrap()
+        }
+    );
+
+    assert_eq!(
+        Fq2 {
+            c0: Fq::from_repr(FqRepr([0xb9f78429d1517a6b, 0x1eabfffeb153ffff, 0x6730d2a0f6b0f624, 0x64774b84f38512bf, 0x4b1ba7b6434bacd7, 0x1a0111ea397fe69a])).unwrap(),
+            c1: Fq::zero()
+        }.sqrt().unwrap(),
+        Fq2 {
+            c0: Fq::zero(),
+            c1: Fq::from_repr(FqRepr([0xb9fefffffd4357a3, 0x1eabfffeb153ffff, 0x6730d2a0f6b0f624, 0x64774b84f38512bf, 0x4b1ba7b6434bacd7, 0x1a0111ea397fe69a])).unwrap()
+        }
+    );
+}
+
+#[cfg(test)]
+use rand::{SeedableRng, XorShiftRng};
+
+#[test]
+fn test_fq2_mul_nonresidue() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let nqr = Fq2 {
+        c0: Fq::one(),
+        c1: Fq::one()
+    };
+
+    for _ in 0..1000 {
+        let mut a = Fq2::rand(&mut rng);
+        let mut b = a;
+        a.mul_by_nonresidue();
+        b.mul_assign(&nqr);
+
+        assert_eq!(a, b);
+    }
+}
+
+#[bench]
+fn bench_fq2_add_assign(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<(Fq2, Fq2)> = (0..SAMPLES).map(|_| {
+        (Fq2::rand(&mut rng), Fq2::rand(&mut rng))
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let mut tmp = v[count].0;
+        tmp.add_assign(&v[count].1);
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_fq2_sub_assign(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<(Fq2, Fq2)> = (0..SAMPLES).map(|_| {
+        (Fq2::rand(&mut rng), Fq2::rand(&mut rng))
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let mut tmp = v[count].0;
+        tmp.sub_assign(&v[count].1);
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_fq2_mul_assign(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<(Fq2, Fq2)> = (0..SAMPLES).map(|_| {
+        (Fq2::rand(&mut rng), Fq2::rand(&mut rng))
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let mut tmp = v[count].0;
+        tmp.mul_assign(&v[count].1);
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_fq2_squaring(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<Fq2> = (0..SAMPLES).map(|_| {
+        Fq2::rand(&mut rng)
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let mut tmp = v[count];
+        tmp.square();
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_fq2_inverse(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<Fq2> = (0..SAMPLES).map(|_| {
+        Fq2::rand(&mut rng)
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let tmp = v[count].inverse();
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_fq2_sqrt(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<Fq2> = (0..SAMPLES).map(|_| {
+        Fq2::rand(&mut rng)
+    }).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let tmp = v[count].sqrt();
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[test]
+fn fq2_field_tests() {
+    use ::PrimeField;
+    
+    ::tests::field::random_field_tests::<Fq2>();
+    ::tests::field::random_sqrt_tests::<Fq2>();
+    ::tests::field::random_frobenius_tests::<Fq2, _>(super::fq::Fq::char(), 13);
+}

src/bls12_381/fq6.rs

@@ -0,0 +1,372 @@
+use rand::{Rng, Rand};
+use ::{Field};
+use super::fq2::Fq2;
+use super::fq::{FROBENIUS_COEFF_FQ6_C1, FROBENIUS_COEFF_FQ6_C2};
+
+/// An element of F_{q^6}, represented by c0 + c1 * v + c2 * v^2.
+#[derive(Copy, Clone, Debug, Eq, PartialEq)]
+pub struct Fq6 {
+    pub c0: Fq2,
+    pub c1: Fq2,
+    pub c2: Fq2
+}
+
+impl Rand for Fq6 {
+    fn rand<R: Rng>(rng: &mut R) -> Self {
+        Fq6 {
+            c0: rng.gen(),
+            c1: rng.gen(),
+            c2: rng.gen()
+        }
+    }
+}
+
+impl Fq6 {
+    /// Multiply by quadratic nonresidue v.
+    pub fn mul_by_nonresidue(&mut self) {
+        use std::mem::swap;
+        swap(&mut self.c0, &mut self.c1);
+        swap(&mut self.c0, &mut self.c2);
+
+        self.c0.mul_by_nonresidue();
+    }
+
+    pub fn mul_by_1(&mut self, c1: &Fq2)
+    {
+        let mut b_b = self.c1;
+        b_b.mul_assign(c1);
+
+        let mut t1 = *c1;
+        {
+            let mut tmp = self.c1;
+            tmp.add_assign(&self.c2);
+
+            t1.mul_assign(&tmp);
+            t1.sub_assign(&b_b);
+            t1.mul_by_nonresidue();
+        }
+
+        let mut t2 = *c1;
+        {
+            let mut tmp = self.c0;
+            tmp.add_assign(&self.c1);
+
+            t2.mul_assign(&tmp);
+            t2.sub_assign(&b_b);
+        }
+
+        self.c0 = t1;
+        self.c1 = t2;
+        self.c2 = b_b;
+    }
+
+    pub fn mul_by_01(&mut self, c0: &Fq2, c1: &Fq2)
+    {
+        let mut a_a = self.c0;
+        let mut b_b = self.c1;
+        a_a.mul_assign(c0);
+        b_b.mul_assign(c1);
+
+        let mut t1 = *c1;
+        {
+            let mut tmp = self.c1;
+            tmp.add_assign(&self.c2);
+
+            t1.mul_assign(&tmp);
+            t1.sub_assign(&b_b);
+            t1.mul_by_nonresidue();
+            t1.add_assign(&a_a);
+        }
+
+        let mut t3 = *c0;
+        {
+            let mut tmp = self.c0;
+            tmp.add_assign(&self.c2);
+
+            t3.mul_assign(&tmp);
+            t3.sub_assign(&a_a);
+            t3.add_assign(&b_b);
+        }
+
+        let mut t2 = *c0;
+        t2.add_assign(c1);
+        {
+            let mut tmp = self.c0;
+            tmp.add_assign(&self.c1);
+
+            t2.mul_assign(&tmp);
+            t2.sub_assign(&a_a);
+            t2.sub_assign(&b_b);
+        }
+
+        self.c0 = t1;
+        self.c1 = t2;
+        self.c2 = t3;
+    }
+}
+
+impl Field for Fq6
+{
+    fn zero() -> Self {
+        Fq6 {
+            c0: Fq2::zero(),
+            c1: Fq2::zero(),
+            c2: Fq2::zero()
+        }
+    }
+
+    fn one() -> Self {
+        Fq6 {
+            c0: Fq2::one(),
+            c1: Fq2::zero(),
+            c2: Fq2::zero()
+        }
+    }
+
+    fn is_zero(&self) -> bool {
+        self.c0.is_zero() && self.c1.is_zero() && self.c2.is_zero()
+    }
+
+    fn double(&mut self) {
+        self.c0.double();
+        self.c1.double();
+        self.c2.double();
+    }
+
+    fn negate(&mut self) {
+        self.c0.negate();
+        self.c1.negate();
+        self.c2.negate();
+    }
+
+    fn add_assign(&mut self, other: &Self) {
+        self.c0.add_assign(&other.c0);
+        self.c1.add_assign(&other.c1);
+        self.c2.add_assign(&other.c2);
+    }
+
+    fn sub_assign(&mut self, other: &Self) {
+        self.c0.sub_assign(&other.c0);
+        self.c1.sub_assign(&other.c1);
+        self.c2.sub_assign(&other.c2);
+    }
+
+    fn frobenius_map(&mut self, power: usize)
+    {
+        self.c0.frobenius_map(power);
+        self.c1.frobenius_map(power);
+        self.c2.frobenius_map(power);
+
+        self.c1.mul_assign(&FROBENIUS_COEFF_FQ6_C1[power % 6]);
+        self.c2.mul_assign(&FROBENIUS_COEFF_FQ6_C2[power % 6]);
+    }
+
+    fn square(&mut self) {
+        let mut s0 = self.c0;
+        s0.square();
+        let mut ab = self.c0;
+        ab.mul_assign(&self.c1);
+        let mut s1 = ab;
+        s1.double();
+        let mut s2 = self.c0;
+        s2.sub_assign(&self.c1);
+        s2.add_assign(&self.c2);
+        s2.square();
+        let mut bc = self.c1;
+        bc.mul_assign(&self.c2);
+        let mut s3 = bc;
+        s3.double();
+        let mut s4 = self.c2;
+        s4.square();
+
+        self.c0 = s3;
+        self.c0.mul_by_nonresidue();
+        self.c0.add_assign(&s0);
+
+        self.c1 = s4;
+        self.c1.mul_by_nonresidue();
+        self.c1.add_assign(&s1);
+
+        self.c2 = s1;
+        self.c2.add_assign(&s2);
+        self.c2.add_assign(&s3);
+        self.c2.sub_assign(&s0);
+        self.c2.sub_assign(&s4);
+    }
+
+    fn mul_assign(&mut self, other: &Self) {
+        let mut a_a = self.c0;
+        let mut b_b = self.c1;
+        let mut c_c = self.c2;
+        a_a.mul_assign(&other.c0);
+        b_b.mul_assign(&other.c1);
+        c_c.mul_assign(&other.c2);
+
+        let mut t1 = other.c1;
+        t1.add_assign(&other.c2);
+        {
+            let mut tmp = self.c1;
+            tmp.add_assign(&self.c2);
+
+            t1.mul_assign(&tmp);
+            t1.sub_assign(&b_b);
+            t1.sub_assign(&c_c);
+            t1.mul_by_nonresidue();
+            t1.add_assign(&a_a);
+        }
+
+        let mut t3 = other.c0;
+        t3.add_assign(&other.c2);
+        {
+            let mut tmp = self.c0;
+            tmp.add_assign(&self.c2);
+
+            t3.mul_assign(&tmp);
+            t3.sub_assign(&a_a);
+            t3.add_assign(&b_b);
+            t3.sub_assign(&c_c);
+        }
+
+        let mut t2 = other.c0;
+        t2.add_assign(&other.c1);
+        {
+            let mut tmp = self.c0;
+            tmp.add_assign(&self.c1);
+
+            t2.mul_assign(&tmp);
+            t2.sub_assign(&a_a);
+            t2.sub_assign(&b_b);
+            c_c.mul_by_nonresidue();
+            t2.add_assign(&c_c);
+        }
+
+        self.c0 = t1;
+        self.c1 = t2;
+        self.c2 = t3;
+    }
+
+    fn inverse(&self) -> Option<Self> {
+        let mut c0 = self.c2;
+        c0.mul_by_nonresidue();
+        c0.mul_assign(&self.c1);
+        c0.negate();
+        {
+            let mut c0s = self.c0;
+            c0s.square();
+            c0.add_assign(&c0s);
+        }
+        let mut c1 = self.c2;
+        c1.square();
+        c1.mul_by_nonresidue();
+        {
+            let mut c01 = self.c0;
+            c01.mul_assign(&self.c1);
+            c1.sub_assign(&c01);
+        }
+        let mut c2 = self.c1;
+        c2.square();
+        {
+            let mut c02 = self.c0;
+            c02.mul_assign(&self.c2);
+            c2.sub_assign(&c02);
+        }
+
+        let mut tmp1 = self.c2;
+        tmp1.mul_assign(&c1);
+        let mut tmp2 = self.c1;
+        tmp2.mul_assign(&c2);
+        tmp1.add_assign(&tmp2);
+        tmp1.mul_by_nonresidue();
+        tmp2 = self.c0;
+        tmp2.mul_assign(&c0);
+        tmp1.add_assign(&tmp2);
+
+        match tmp1.inverse() {
+            Some(t) => {
+                let mut tmp = Fq6 {
+                    c0: t,
+                    c1: t,
+                    c2: t
+                };
+                tmp.c0.mul_assign(&c0);
+                tmp.c1.mul_assign(&c1);
+                tmp.c2.mul_assign(&c2);
+
+                Some(tmp)
+            },
+            None => None
+        }
+    }
+}
+
+#[cfg(test)]
+use rand::{SeedableRng, XorShiftRng};
+
+#[test]
+fn test_fq6_mul_nonresidue() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let nqr = Fq6 {
+        c0: Fq2::zero(),
+        c1: Fq2::one(),
+        c2: Fq2::zero()
+    };
+
+    for _ in 0..1000 {
+        let mut a = Fq6::rand(&mut rng);
+        let mut b = a;
+        a.mul_by_nonresidue();
+        b.mul_assign(&nqr);
+
+        assert_eq!(a, b);
+    }
+}
+
+#[test]
+fn test_fq6_mul_by_1() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..1000 {
+        let c1 = Fq2::rand(&mut rng);
+        let mut a = Fq6::rand(&mut rng);
+        let mut b = a;
+
+        a.mul_by_1(&c1);
+        b.mul_assign(&Fq6 {
+            c0: Fq2::zero(),
+            c1: c1,
+            c2: Fq2::zero()
+        });
+
+        assert_eq!(a, b);
+    }
+}
+
+#[test]
+fn test_fq6_mul_by_01() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..1000 {
+        let c0 = Fq2::rand(&mut rng);
+        let c1 = Fq2::rand(&mut rng);
+        let mut a = Fq6::rand(&mut rng);
+        let mut b = a;
+
+        a.mul_by_01(&c0, &c1);
+        b.mul_assign(&Fq6 {
+            c0: c0,
+            c1: c1,
+            c2: Fq2::zero()
+        });
+
+        assert_eq!(a, b);
+    }
+}
+
+#[test]
+fn fq6_field_tests() {
+    use ::PrimeField;
+    
+    ::tests::field::random_field_tests::<Fq6>();
+    ::tests::field::random_frobenius_tests::<Fq6, _>(super::fq::Fq::char(), 13);
+}

src/bls12_381/mod.rs

@@ -0,0 +1,459 @@
+mod fq;
+mod fr;
+mod fq2;
+mod fq6;
+mod fq12;
+mod ec;
+
+pub use self::fr::{Fr, FrRepr};
+pub use self::fq::{Fq, FqRepr};
+pub use self::fq2::Fq2;
+pub use self::fq12::Fq12;
+pub use self::ec::{G1, G2, G1Affine, G2Affine, G1Prepared, G2Prepared};
+
+use super::{Engine, CurveAffine, Field, BitIterator};
+
+// The BLS parameter x for BLS12-381 is -0xd201000000010000
+const BLS_X: u64 = 0xd201000000010000;
+const BLS_X_IS_NEGATIVE: bool = true;
+
+pub struct Bls12;
+
+impl Engine for Bls12 {
+    type Fr = Fr;
+    type G1 = G1;
+    type G1Affine = G1Affine;
+    type G2 = G2;
+    type G2Affine = G2Affine;
+    type Fq = Fq;
+    type Fqe = Fq2;
+    type Fqk = Fq12;
+
+    fn miller_loop<'a, I>(i: I) -> Self::Fqk
+        where I: IntoIterator<Item=&'a (
+                                    &'a <Self::G1Affine as CurveAffine>::Prepared,
+                                    &'a <Self::G2Affine as CurveAffine>::Prepared
+                               )>
+    {
+        let mut pairs = vec![];
+        for &(p, q) in i {
+            if !p.is_zero() && !q.is_zero() {
+                pairs.push((p, q.coeffs.iter()));
+            }
+        }
+
+        // Twisting isomorphism from E to E'
+        fn ell(
+            f: &mut Fq12,
+            coeffs: &(Fq2, Fq2, Fq2),
+            p: &G1Affine
+        )
+        {
+            let mut c0 = coeffs.0;
+            let mut c1 = coeffs.1;
+
+            c0.c0.mul_assign(&p.y);
+            c0.c1.mul_assign(&p.y);
+
+            c1.c0.mul_assign(&p.x);
+            c1.c1.mul_assign(&p.x);
+
+            // Sparse multiplication in Fq12
+            f.mul_by_014(&coeffs.2, &c1, &c0);
+        }
+
+        let mut f = Fq12::one();
+
+        let mut found_one = false;
+        for i in BitIterator::new(&[BLS_X >> 1]) {
+            if !found_one {
+                found_one = i;
+                continue;
+            }
+
+            for &mut (p, ref mut coeffs) in &mut pairs {
+                ell(&mut f, coeffs.next().unwrap(), &p.0);
+            }
+
+            if i {
+                for &mut (p, ref mut coeffs) in &mut pairs {
+                    ell(&mut f, coeffs.next().unwrap(), &p.0);
+                }
+            }
+
+            f.square();
+        }
+
+        for &mut (p, ref mut coeffs) in &mut pairs {
+            ell(&mut f, coeffs.next().unwrap(), &p.0);
+        }
+
+        f
+    }
+
+    fn final_exponentiation(r: &Fq12) -> Option<Fq12> {
+        let mut f1 = *r;
+        f1.unitary_inverse();
+
+        match r.inverse() {
+            Some(mut f2) => {
+                let mut r = f1;
+                r.mul_assign(&f2);
+                f2 = r;
+                r.frobenius_map(2);
+                r.mul_assign(&f2);
+
+                fn exp_by_x(f: &mut Fq12, x: u64)
+                {
+                    *f = f.pow(&[x]);
+                    if BLS_X_IS_NEGATIVE {
+                        f.unitary_inverse();
+                    }
+                }
+
+                let mut x = BLS_X;
+                let mut y0 = r;
+                y0.square();
+                let mut y1 = y0;
+                exp_by_x(&mut y1, x);
+                x >>= 1;
+                let mut y2 = y1;
+                exp_by_x(&mut y2, x);
+                x <<= 1;
+                let mut y3 = r;
+                y3.unitary_inverse();
+                y1.mul_assign(&y3);
+                y1.unitary_inverse();
+                y1.mul_assign(&y2);
+                y2 = y1;
+                exp_by_x(&mut y2, x);
+                y3 = y2;
+                exp_by_x(&mut y3, x);
+                y1.unitary_inverse();
+                y3.mul_assign(&y1);
+                y1.unitary_inverse();
+                y1.frobenius_map(3);
+                y2.frobenius_map(2);
+                y1.mul_assign(&y2);
+                y2 = y3;
+                exp_by_x(&mut y2, x);
+                y2.mul_assign(&y0);
+                y2.mul_assign(&r);
+                y1.mul_assign(&y2);
+                y2 = y3;
+                y2.frobenius_map(1);
+                y1.mul_assign(&y2);
+                Some(y1)
+            },
+            None => None
+        }
+    }
+}
+
+impl G2Prepared {
+    pub fn is_zero(&self) -> bool {
+        self.infinity
+    }
+
+    pub fn from_affine(q: G2Affine) -> Self {
+        if q.is_zero() {
+            return G2Prepared {
+                coeffs: vec![],
+                infinity: true
+            }
+        }
+
+        fn doubling_step(
+            r: &mut G2
+        ) -> (Fq2, Fq2, Fq2)
+        {
+            // Adaptation of Algorithm 26, https://eprint.iacr.org/2010/354.pdf
+            let mut tmp0 = r.x;
+            tmp0.square();
+
+            let mut tmp1 = r.y;
+            tmp1.square();
+
+            let mut tmp2 = tmp1;
+            tmp2.square();
+
+            let mut tmp3 = tmp1;
+            tmp3.add_assign(&r.x);
+            tmp3.square();
+            tmp3.sub_assign(&tmp0);
+            tmp3.sub_assign(&tmp2);
+            tmp3.double();
+
+            let mut tmp4 = tmp0;
+            tmp4.double();
+            tmp4.add_assign(&tmp0);
+
+            let mut tmp6 = r.x;
+            tmp6.add_assign(&tmp4);
+
+            let mut tmp5 = tmp4;
+            tmp5.square();
+
+            let mut zsquared = r.z;
+            zsquared.square();
+
+            r.x = tmp5;
+            r.x.sub_assign(&tmp3);
+            r.x.sub_assign(&tmp3);
+
+            r.z.add_assign(&r.y);
+            r.z.square();
+            r.z.sub_assign(&tmp1);
+            r.z.sub_assign(&zsquared);
+
+            r.y = tmp3;
+            r.y.sub_assign(&r.x);
+            r.y.mul_assign(&tmp4);
+
+            tmp2.double();
+            tmp2.double();
+            tmp2.double();
+
+            r.y.sub_assign(&tmp2);
+
+            tmp3 = tmp4;
+            tmp3.mul_assign(&zsquared);
+            tmp3.double();
+            tmp3.negate();
+
+            tmp6.square();
+            tmp6.sub_assign(&tmp0);
+            tmp6.sub_assign(&tmp5);
+
+            tmp1.double();
+            tmp1.double();
+
+            tmp6.sub_assign(&tmp1);
+
+            tmp0 = r.z;
+            tmp0.mul_assign(&zsquared);
+            tmp0.double();
+
+            (tmp0, tmp3, tmp6)
+        }
+
+        fn addition_step(
+            r: &mut G2,
+            q: &G2Affine
+        ) -> (Fq2, Fq2, Fq2)
+        {
+            // Adaptation of Algorithm 27, https://eprint.iacr.org/2010/354.pdf
+            let mut zsquared = r.z;
+            zsquared.square();
+
+            let mut ysquared = q.y;
+            ysquared.square();
+
+            let mut t0 = zsquared;
+            t0.mul_assign(&q.x);
+
+            let mut t1 = q.y;
+            t1.add_assign(&r.z);
+            t1.square();
+            t1.sub_assign(&ysquared);
+            t1.sub_assign(&zsquared);
+            t1.mul_assign(&zsquared);
+
+            let mut t2 = t0;
+            t2.sub_assign(&r.x);
+
+            let mut t3 = t2;
+            t3.square();
+
+            let mut t4 = t3;
+            t4.double();
+            t4.double();
+
+            let mut t5 = t4;
+            t5.mul_assign(&t2);
+
+            let mut t6 = t1;
+            t6.sub_assign(&r.y);
+            t6.sub_assign(&r.y);
+
+            let mut t9 = t6;
+            t9.mul_assign(&q.x);
+
+            let mut t7 = t4;
+            t7.mul_assign(&r.x);
+
+            r.x = t6;
+            r.x.square();
+            r.x.sub_assign(&t5);
+            r.x.sub_assign(&t7);
+            r.x.sub_assign(&t7);
+
+            r.z.add_assign(&t2);
+            r.z.square();
+            r.z.sub_assign(&zsquared);
+            r.z.sub_assign(&t3);
+
+            let mut t10 = q.y;
+            t10.add_assign(&r.z);
+
+            let mut t8 = t7;
+            t8.sub_assign(&r.x);
+            t8.mul_assign(&t6);
+
+            t0 = r.y;
+            t0.mul_assign(&t5);
+            t0.double();
+
+            r.y = t8;
+            r.y.sub_assign(&t0);
+
+            t10.square();
+            t10.sub_assign(&ysquared);
+
+            let mut ztsquared = r.z;
+            ztsquared.square();
+
+            t10.sub_assign(&ztsquared);
+
+            t9.double();
+            t9.sub_assign(&t10);
+
+            t10 = r.z;
+            t10.double();
+
+            t6.negate();
+
+            t1 = t6;
+            t1.double();
+
+            (t10, t1, t9)
+        }
+
+        let mut coeffs = vec![];
+        let mut r: G2 = q.into();
+
+        let mut found_one = false;
+        for i in BitIterator::new([BLS_X >> 1]) {
+            if !found_one {
+                found_one = i;
+                continue;
+            }
+
+            coeffs.push(doubling_step(&mut r));
+
+            if i {
+                coeffs.push(addition_step(&mut r, &q));
+            }
+        }
+
+        coeffs.push(doubling_step(&mut r));
+
+        G2Prepared {
+            coeffs: coeffs,
+            infinity: false
+        }
+    }
+}
+
+#[test]
+fn bls12_engine_tests() {
+    ::tests::engine::engine_tests::<Bls12>();
+}
+
+#[cfg(test)]
+use rand::{Rand, SeedableRng, XorShiftRng};
+
+#[bench]
+fn bench_pairing_g1_preparation(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<G1> = (0..SAMPLES).map(|_| G1::rand(&mut rng)).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let tmp = G1Affine::from(v[count]).prepare();
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_pairing_g2_preparation(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<G2> = (0..SAMPLES).map(|_| G2::rand(&mut rng)).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let tmp = G2Affine::from(v[count]).prepare();
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_pairing_miller_loop(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<(G1Prepared, G2Prepared)> = (0..SAMPLES).map(|_|
+        (
+            G1Affine::from(G1::rand(&mut rng)).prepare(),
+            G2Affine::from(G2::rand(&mut rng)).prepare()
+        )
+    ).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let tmp = Bls12::miller_loop(&[(&v[count].0, &v[count].1)]);
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_pairing_final_exponentiation(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<Fq12> = (0..SAMPLES).map(|_|
+        (
+            G1Affine::from(G1::rand(&mut rng)).prepare(),
+            G2Affine::from(G2::rand(&mut rng)).prepare()
+        )
+    ).map(|(ref p, ref q)| Bls12::miller_loop(&[(p, q)])).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let tmp = Bls12::final_exponentiation(&v[count]);
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}
+
+#[bench]
+fn bench_pairing_full(b: &mut ::test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    let v: Vec<(G1, G2)> = (0..SAMPLES).map(|_|
+        (
+            G1::rand(&mut rng),
+            G2::rand(&mut rng)
+        )
+    ).collect();
+
+    let mut count = 0;
+    b.iter(|| {
+        let tmp = Bls12::pairing(v[count].0, v[count].1);
+        count = (count + 1) % SAMPLES;
+        tmp
+    });
+}

src/lib.rs

@@ -0,0 +1,409 @@
+#![feature(i128_type)]
+
+#![cfg_attr(test, feature(test))]
+#[cfg(test)]
+extern crate test;
+
+extern crate rand;
+
+#[cfg(test)]
+pub mod tests;
+
+pub mod bls12_381;
+
+use std::fmt;
+
+/// An "engine" is a collection of types (fields, elliptic curve groups, etc.)
+/// with well-defined relationships. In particular, the G1/G2 curve groups are
+/// of prime order `r`, and are equipped with a bilinear pairing function.
+pub trait Engine {
+    /// This is the scalar field of the G1/G2 groups.
+    type Fr: PrimeField;
+
+    /// The projective representation of an element in G1.
+    type G1: CurveProjective<Base=Self::Fq, Scalar=Self::Fr, Affine=Self::G1Affine> + From<Self::G1Affine>;
+
+    /// The affine representation of an element in G1.
+    type G1Affine: CurveAffine<Base=Self::Fq, Scalar=Self::Fr, Projective=Self::G1> + From<Self::G1>;
+
+    /// The projective representation of an element in G2.
+    type G2: CurveProjective<Base=Self::Fqe, Scalar=Self::Fr, Affine=Self::G2Affine> + From<Self::G2Affine>;
+
+    /// The affine representation of an element in G2.
+    type G2Affine: CurveAffine<Base=Self::Fqe, Scalar=Self::Fr, Projective=Self::G2> + From<Self::G2>;
+
+    /// The base field that hosts G1.
+    type Fq: PrimeField + SqrtField;
+
+    /// The extension field that hosts G2.
+    type Fqe: SqrtField;
+
+    /// The extension field that hosts the target group of the pairing.
+    type Fqk: Field;
+
+    /// Perform a miller loop with some number of (G1, G2) pairs.
+    fn miller_loop<'a, I>(i: I) -> Self::Fqk
+        where I: IntoIterator<Item=&'a (
+                                    &'a <Self::G1Affine as CurveAffine>::Prepared,
+                                    &'a <Self::G2Affine as CurveAffine>::Prepared
+                               )>;
+
+    /// Perform final exponentiation of the result of a miller loop.
+    fn final_exponentiation(&Self::Fqk) -> Option<Self::Fqk>;
+
+    /// Performs a complete pairing operation `(p, q)`.
+    fn pairing<G1, G2>(p: G1, q: G2) -> Self::Fqk
+        where G1: Into<Self::G1Affine>,
+              G2: Into<Self::G2Affine>
+    {
+        Self::final_exponentiation(&Self::miller_loop(
+            [(
+                &(p.into().prepare()),
+                &(q.into().prepare())
+            )].into_iter()
+        )).unwrap()
+    }
+}
+
+/// Projective representation of an elliptic curve point guaranteed to be
+/// in the correct prime order subgroup.
+pub trait CurveProjective: PartialEq +
+                           Eq +
+                           Sized +
+                           Copy +
+                           Clone +
+                           Send +
+                           Sync +
+                           fmt::Debug +
+                           rand::Rand +
+                           'static
+{
+    type Scalar: PrimeField;
+    type Base: SqrtField;
+    type Affine: CurveAffine<Projective=Self, Scalar=Self::Scalar>;
+
+    /// Returns the additive identity.
+    fn zero() -> Self;
+
+    /// Returns a fixed generator of unknown exponent.
+    fn one() -> Self;
+
+    /// Determines if this point is the point at infinity.
+    fn is_zero(&self) -> bool;
+
+    /// Normalizes a slice of projective elements so that
+    /// conversion to affine is cheap.
+    fn batch_normalization(v: &mut [Self]);
+
+    /// Checks if the point is already "normalized" so that
+    /// cheap affine conversion is possible.
+    fn is_normalized(&self) -> bool;
+
+    /// Doubles this element.
+    fn double(&mut self);
+
+    /// Adds another element to this element.
+    fn add_assign(&mut self, other: &Self);
+
+    /// Subtracts another element from this element.
+    fn sub_assign(&mut self, other: &Self) {
+        let mut tmp = *other;
+        tmp.negate();
+        self.add_assign(&tmp);
+    }
+
+    /// Adds an affine element to this element.
+    fn add_assign_mixed(&mut self, other: &Self::Affine);
+
+    /// Negates this element.
+    fn negate(&mut self);
+
+    /// Performs scalar multiplication of this element.
+    fn mul_assign<S: Into<<Self::Scalar as PrimeField>::Repr>>(&mut self, other: S);
+
+    /// Converts this element into its affine representation.
+    fn to_affine(&self) -> Self::Affine;
+}
+
+/// Affine representation of an elliptic curve point guaranteed to be
+/// in the correct prime order subgroup.
+pub trait CurveAffine: Copy +
+                       Clone +
+                       Sized +
+                       Send +
+                       Sync +
+                       fmt::Debug +
+                       PartialEq +
+                       Eq +
+                       'static
+{
+    type Scalar: PrimeField;
+    type Base: SqrtField;
+    type Projective: CurveProjective<Affine=Self, Scalar=Self::Scalar>;
+    type Prepared: Clone + Send + Sync + 'static;
+
+    /// Returns the additive identity.
+    fn zero() -> Self;
+
+    /// Returns a fixed generator of unknown exponent.
+    fn one() -> Self;
+
+    /// Determines if this point represents the point at infinity; the
+    /// additive identity.
+    fn is_zero(&self) -> bool;
+
+    /// Determines if this point is on the curve and in the correct subgroup.
+    fn is_valid(&self) -> bool;
+
+    /// Negates this element.
+    fn negate(&mut self);
+
+    /// Performs scalar multiplication of this element with mixed addition.
+    fn mul<S: Into<<Self::Scalar as PrimeField>::Repr>>(&self, other: S) -> Self::Projective;
+
+    /// Prepares this element for pairing purposes.
+    fn prepare(&self) -> Self::Prepared;
+
+    /// Converts this element into its affine representation.
+    fn to_projective(&self) -> Self::Projective;
+}
+
+/// This trait represents an element of a field.
+pub trait Field: Sized +
+                 Eq +
+                 Copy +
+                 Clone +
+                 Send +
+                 Sync +
+                 fmt::Debug +
+                 'static +
+                 rand::Rand
+{
+    /// Returns the zero element of the field, the additive identity.
+    fn zero() -> Self;
+
+    /// Returns the one element of the field, the multiplicative identity.
+    fn one() -> Self;
+
+    /// Returns true iff this element is zero.
+    fn is_zero(&self) -> bool;
+
+    /// Squares this element.
+    fn square(&mut self);
+
+    /// Doubles this element.
+    fn double(&mut self);
+
+    /// Negates this element.
+    fn negate(&mut self);
+
+    /// Adds another element to this element.
+    fn add_assign(&mut self, other: &Self);
+
+    /// Subtracts another element from this element.
+    fn sub_assign(&mut self, other: &Self);
+
+    /// Multiplies another element by this element.
+    fn mul_assign(&mut self, other: &Self);
+
+    /// Computes the multiplicative inverse of this element, if nonzero.
+    fn inverse(&self) -> Option<Self>;
+
+    /// Exponentiates this element by a power of the base prime modulus via
+    /// the Frobenius automorphism.
+    fn frobenius_map(&mut self, power: usize);
+
+    /// Exponentiates this element by a number represented with `u64` limbs,
+    /// least significant digit first.
+    fn pow<S: AsRef<[u64]>>(&self, exp: S) -> Self
+    {
+        let mut res = Self::one();
+
+        for i in BitIterator::new(exp) {
+            res.square();
+            if i {
+                res.mul_assign(self);
+            }
+        }
+
+        res
+    }
+}
+
+/// This trait represents an element of a field that has a square root operation described for it.
+pub trait SqrtField: Field
+{
+    /// Returns the square root of the field element, if it is
+    /// quadratic residue.
+    fn sqrt(&self) -> Option<Self>;
+}
+
+/// This trait represents a wrapper around a biginteger which can encode any element of a particular
+/// prime field. It is a smart wrapper around a sequence of `u64` limbs, least-significant digit
+/// first.
+pub trait PrimeFieldRepr: Sized +
+                          Copy +
+                          Clone +
+                          Eq +
+                          Ord +
+                          Send +
+                          Sync +
+                          fmt::Debug +
+                          'static +
+                          rand::Rand +
+                          AsRef<[u64]> +
+                          From<u64>
+{
+    /// Subtract another reprensetation from this one, returning the borrow bit.
+    fn sub_noborrow(&mut self, other: &Self) -> bool;
+
+    /// Add another representation to this one, returning the carry bit.
+    fn add_nocarry(&mut self, other: &Self) -> bool;
+
+    /// Compute the number of bits needed to encode this number.
+    fn num_bits(&self) -> u32;
+
+    /// Returns true iff this number is zero.
+    fn is_zero(&self) -> bool;
+
+    /// Returns true iff this number is odd.
+    fn is_odd(&self) -> bool;
+
+    /// Returns true iff this number is even.
+    fn is_even(&self) -> bool;
+
+    /// Performs a rightwise bitshift of this number, effectively dividing
+    /// it by 2.
+    fn div2(&mut self);
+
+    /// Performs a rightwise bitshift of this number by some amount.
+    fn divn(&mut self, amt: usize);
+
+    /// Performs a leftwise bitshift of this number, effectively multiplying
+    /// it by 2. Overflow is ignored.
+    fn mul2(&mut self);
+}
+
+/// This represents an element of a prime field.
+pub trait PrimeField: Field
+{
+    /// The prime field can be converted back and forth into this biginteger
+    /// representation.
+    type Repr: PrimeFieldRepr + From<Self>;
+
+    /// Convert this prime field element into a biginteger representation.
+    fn from_repr(Self::Repr) -> Result<Self, ()>;
+
+    /// Convert a biginteger reprensentation into a prime field element, if
+    /// the number is an element of the field.
+    fn into_repr(&self) -> Self::Repr;
+
+    /// Returns the field characteristic; the modulus.
+    fn char() -> Self::Repr;
+
+    /// Returns how many bits are needed to represent an element of this
+    /// field.
+    fn num_bits() -> u32;
+
+    /// Returns how many bits of information can be reliably stored in the
+    /// field element.
+    fn capacity() -> u32;
+
+    /// Returns the multiplicative generator of `char()` - 1 order. This element
+    /// must also be quadratic nonresidue.
+    fn multiplicative_generator() -> Self;
+
+    /// Returns s such that 2^s * t = `char()` - 1 with t odd.
+    fn s() -> usize;
+
+    /// Returns the 2^s root of unity computed by exponentiating the `multiplicative_generator()`
+    /// by t.
+    fn root_of_unity() -> Self;
+}
+
+pub struct BitIterator<E> {
+    t: E,
+    n: usize
+}
+
+impl<E: AsRef<[u64]>> BitIterator<E> {
+    pub fn new(t: E) -> Self {
+        let n = t.as_ref().len() * 64;
+
+        BitIterator {
+            t: t,
+            n: n
+        }
+    }
+}
+
+impl<E: AsRef<[u64]>> Iterator for BitIterator<E> {
+    type Item = bool;
+
+    fn next(&mut self) -> Option<bool> {
+        if self.n == 0 {
+            None
+        } else {
+            self.n -= 1;
+            let part = self.n / 64;
+            let bit = self.n - (64 * part);
+
+            Some(self.t.as_ref()[part] & (1 << bit) > 0)
+        }
+    }
+}
+
+#[test]
+fn test_bit_iterator() {
+    let mut a = BitIterator::new([0xa953d79b83f6ab59, 0x6dea2059e200bd39]);
+    let expected = "01101101111010100010000001011001111000100000000010111101001110011010100101010011110101111001101110000011111101101010101101011001";
+
+    for e in expected.chars() {
+        assert!(a.next().unwrap() == (e == '1'));
+    }
+
+    assert!(a.next().is_none());
+
+    let expected = "1010010101111110101010000101101011101000011101110101001000011001100100100011011010001011011011010001011011101100110100111011010010110001000011110100110001100110011101101000101100011100100100100100001010011101010111110011101011000011101000111011011101011001";
+
+    let mut a = BitIterator::new([0x429d5f3ac3a3b759, 0xb10f4c66768b1c92, 0x92368b6d16ecd3b4, 0xa57ea85ae8775219]);
+
+    for e in expected.chars() {
+        assert!(a.next().unwrap() == (e == '1'));
+    }
+
+    assert!(a.next().is_none());
+}
+
+/// Calculate a - b - borrow, returning the result and modifying
+/// the borrow value.
+#[inline(always)]
+pub(crate) fn sbb(a: u64, b: u64, borrow: &mut u64) -> u64 {
+    let tmp = (1u128 << 64) + (a as u128) - (b as u128) - (*borrow as u128);
+
+    *borrow = if tmp >> 64 == 0 { 1 } else { 0 };
+
+    tmp as u64
+}
+
+/// Calculate a + b + carry, returning the sum and modifying the
+/// carry value.
+#[inline(always)]
+pub(crate) fn adc(a: u64, b: u64, carry: &mut u64) -> u64 {
+    let tmp = (a as u128) + (b as u128) + (*carry as u128);
+
+    *carry = (tmp >> 64) as u64;
+
+    tmp as u64
+}
+
+/// Calculate a + (b * c) + carry, returning the least significant digit
+/// and setting carry to the most significant digit.
+#[inline(always)]
+pub(crate) fn mac_with_carry(a: u64, b: u64, c: u64, carry: &mut u64) -> u64 {
+    let tmp = (a as u128) + (b as u128) * (c as u128) + (*carry as u128);
+
+    *carry = (tmp >> 64) as u64;
+
+    tmp as u64
+}

src/tests/curve.rs

@@ -0,0 +1,267 @@
+use rand::{SeedableRng, XorShiftRng, Rand};
+
+use ::{CurveProjective, CurveAffine, Field};
+
+pub fn curve_tests<G: CurveProjective>()
+{
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    // Negation edge case with zero.
+    {
+        let mut z = G::zero();
+        z.negate();
+        assert!(z.is_zero());
+    }
+
+    // Doubling edge case with zero.
+    {
+        let mut z = G::zero();
+        z.double();
+        assert!(z.is_zero());
+    }
+
+    // Addition edge cases with zero
+    {
+        let mut r = G::rand(&mut rng);
+        let rcopy = r;
+        r.add_assign(&G::zero());
+        assert_eq!(r, rcopy);
+        r.add_assign_mixed(&G::Affine::zero());
+        assert_eq!(r, rcopy);
+
+        let mut z = G::zero();
+        z.add_assign(&G::zero());
+        assert!(z.is_zero());
+        z.add_assign_mixed(&G::Affine::zero());
+        assert!(z.is_zero());
+
+        let mut z2 = z;
+        z2.add_assign(&r);
+
+        z.add_assign_mixed(&r.to_affine());
+
+        assert_eq!(z, z2);
+        assert_eq!(z, r);
+    }
+
+    // Transformations
+    {
+        let a = G::rand(&mut rng);
+        let b = a.to_affine().to_projective();
+        let c = a.to_affine().to_projective().to_affine().to_projective();
+        assert_eq!(a, b);
+        assert_eq!(b, c);
+    }
+
+    random_addition_tests::<G>();
+    random_multiplication_tests::<G>();
+    random_doubling_tests::<G>();
+    random_negation_tests::<G>();
+    random_transformation_tests::<G>();
+}
+
+fn random_negation_tests<G: CurveProjective>() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..1000 {
+        let r = G::rand(&mut rng);
+
+        let s = G::Scalar::rand(&mut rng);
+        let mut sneg = s;
+        sneg.negate();
+
+        let mut t1 = r;
+        t1.mul_assign(s);
+
+        let mut t2 = r;
+        t2.mul_assign(sneg);
+
+        let mut t3 = t1;
+        t3.add_assign(&t2);
+        assert!(t3.is_zero());
+
+        let mut t4 = t1;
+        t4.add_assign_mixed(&t2.to_affine());
+        assert!(t4.is_zero());
+
+        t1.negate();
+        assert_eq!(t1, t2);
+    }
+}
+
+fn random_doubling_tests<G: CurveProjective>() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..1000 {
+        let mut a = G::rand(&mut rng);
+        let mut b = G::rand(&mut rng);
+
+        // 2(a + b)
+        let mut tmp1 = a;
+        tmp1.add_assign(&b);
+        tmp1.double();
+
+        // 2a + 2b
+        a.double();
+        b.double();
+
+        let mut tmp2 = a;
+        tmp2.add_assign(&b);
+
+        let mut tmp3 = a;
+        tmp3.add_assign_mixed(&b.to_affine());
+
+        assert_eq!(tmp1, tmp2);
+        assert_eq!(tmp1, tmp3);
+    }
+}
+
+fn random_multiplication_tests<G: CurveProjective>() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..1000 {
+        let mut a = G::rand(&mut rng);
+        let mut b = G::rand(&mut rng);
+        let a_affine = a.to_affine();
+        let b_affine = b.to_affine();
+
+        let s = G::Scalar::rand(&mut rng);
+
+        // s ( a + b )
+        let mut tmp1 = a;
+        tmp1.add_assign(&b);
+        tmp1.mul_assign(s);
+
+        // sa + sb
+        a.mul_assign(s);
+        b.mul_assign(s);
+
+        let mut tmp2 = a;
+        tmp2.add_assign(&b);
+
+        // Affine multiplication
+        let mut tmp3 = a_affine.mul(s);
+        tmp3.add_assign(&b_affine.mul(s));
+
+        assert_eq!(tmp1, tmp2);
+        assert_eq!(tmp1, tmp3);
+    }
+}
+
+fn random_addition_tests<G: CurveProjective>() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..1000 {
+        let a = G::rand(&mut rng);
+        let b = G::rand(&mut rng);
+        let c = G::rand(&mut rng);
+        let a_affine = a.to_affine();
+        let b_affine = b.to_affine();
+        let c_affine = c.to_affine();
+
+        // a + a should equal the doubling
+        {
+            let mut aplusa = a;
+            aplusa.add_assign(&a);
+
+            let mut aplusamixed = a;
+            aplusamixed.add_assign_mixed(&a.to_affine());
+
+            let mut adouble = a;
+            adouble.double();
+
+            assert_eq!(aplusa, adouble);
+            assert_eq!(aplusa, aplusamixed);
+        }
+
+        let mut tmp = vec![G::zero(); 6];
+
+        // (a + b) + c
+        tmp[0] = a;
+        tmp[0].add_assign(&b);
+        tmp[0].add_assign(&c);
+
+        // a + (b + c)
+        tmp[1] = b;
+        tmp[1].add_assign(&c);
+        tmp[1].add_assign(&a);
+
+        // (a + c) + b
+        tmp[2] = a;
+        tmp[2].add_assign(&c);
+        tmp[2].add_assign(&b);
+
+        // Mixed addition
+
+        // (a + b) + c
+        tmp[3] = a_affine.to_projective();
+        tmp[3].add_assign_mixed(&b_affine);
+        tmp[3].add_assign_mixed(&c_affine);
+
+        // a + (b + c)
+        tmp[4] = b_affine.to_projective();
+        tmp[4].add_assign_mixed(&c_affine);
+        tmp[4].add_assign_mixed(&a_affine);
+
+        // (a + c) + b
+        tmp[5] = a_affine.to_projective();
+        tmp[5].add_assign_mixed(&c_affine);
+        tmp[5].add_assign_mixed(&b_affine);
+
+        // Comparisons
+        for i in 0..6 {
+            for j in 0..6 {
+                assert_eq!(tmp[i], tmp[j]);
+                assert_eq!(tmp[i].to_affine(), tmp[j].to_affine());
+            }
+
+            assert!(tmp[i] != a);
+            assert!(tmp[i] != b);
+            assert!(tmp[i] != c);
+
+            assert!(a != tmp[i]);
+            assert!(b != tmp[i]);
+            assert!(c != tmp[i]);
+        }
+    }
+}
+
+fn random_transformation_tests<G: CurveProjective>() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..1000 {
+        let g = G::rand(&mut rng);
+        let g_affine = g.to_affine();
+        let g_projective = g_affine.to_projective();
+        assert_eq!(g, g_projective);
+    }
+
+    // Batch normalization
+    for _ in 0..10 {
+        let mut v = (0..1000).map(|_| G::rand(&mut rng)).collect::<Vec<_>>();
+
+        for i in &v {
+            assert!(!i.is_normalized());
+        }
+
+        use rand::distributions::{IndependentSample, Range};
+        let between = Range::new(0, 1000);
+        // Sprinkle in some normalized points
+        for _ in 0..5 {
+            v[between.ind_sample(&mut rng)] = G::zero();
+        }
+        for _ in 0..5 {
+            let s = between.ind_sample(&mut rng);
+            v[s] = v[s].to_affine().to_projective();
+        }
+
+        let expected_v = v.iter().map(|v| v.to_affine().to_projective()).collect::<Vec<_>>();
+        G::batch_normalization(&mut v);
+
+        for i in &v {
+            assert!(i.is_normalized());
+        }
+
+        assert_eq!(v, expected_v);
+    }
+}

src/tests/engine.rs

@@ -0,0 +1,120 @@
+use rand::{SeedableRng, XorShiftRng, Rand};
+
+use ::{Engine, CurveProjective, CurveAffine, Field, PrimeField};
+
+pub fn engine_tests<E: Engine>()
+{
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+    
+    for _ in 0..1000 {
+        let z1 = E::G1Affine::zero().prepare();
+        let z2 = E::G2Affine::zero().prepare();
+
+        let a = E::G1::rand(&mut rng).to_affine().prepare();
+        let b = E::G2::rand(&mut rng).to_affine().prepare();
+        let c = E::G1::rand(&mut rng).to_affine().prepare();
+        let d = E::G2::rand(&mut rng).to_affine().prepare();
+
+        assert_eq!(
+            E::Fqk::one(),
+            E::final_exponentiation(&E::miller_loop(&[(&z1, &b)])).unwrap()
+        );
+
+        assert_eq!(
+            E::Fqk::one(),
+            E::final_exponentiation(&E::miller_loop(&[(&a, &z2)])).unwrap()
+        );
+
+        assert_eq!(
+            E::final_exponentiation(&E::miller_loop(&[(&z1, &b), (&c, &d)])).unwrap(),
+            E::final_exponentiation(&E::miller_loop(&[(&a, &z2), (&c, &d)])).unwrap()
+        );
+
+        assert_eq!(
+            E::final_exponentiation(&E::miller_loop(&[(&a, &b), (&z1, &d)])).unwrap(),
+            E::final_exponentiation(&E::miller_loop(&[(&a, &b), (&c, &z2)])).unwrap()
+        );
+    }
+
+    random_bilinearity_tests::<E>();
+    random_miller_loop_tests::<E>();
+}
+
+fn random_miller_loop_tests<E: Engine>() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    // Exercise the miller loop for a reduced pairing
+    for _ in 0..1000 {
+        let a = E::G1::rand(&mut rng);
+        let b = E::G2::rand(&mut rng);
+
+        let p2 = E::pairing(a, b);
+
+        let a = a.to_affine().prepare();
+        let b = b.to_affine().prepare();
+
+        let p1 = E::final_exponentiation(&E::miller_loop(&[(&a, &b)])).unwrap();
+
+        assert_eq!(p1, p2);
+    }
+
+    // Exercise a double miller loop
+    for _ in 0..1000 {
+        let a = E::G1::rand(&mut rng);
+        let b = E::G2::rand(&mut rng);
+        let c = E::G1::rand(&mut rng);
+        let d = E::G2::rand(&mut rng);
+
+        let ab = E::pairing(a, b);
+        let cd = E::pairing(c, d);
+
+        let mut abcd = ab;
+        abcd.mul_assign(&cd);
+
+        let a = a.to_affine().prepare();
+        let b = b.to_affine().prepare();
+        let c = c.to_affine().prepare();
+        let d = d.to_affine().prepare();
+
+        let abcd_with_double_loop = E::final_exponentiation(
+            &E::miller_loop(&[(&a, &b), (&c, &d)])
+        ).unwrap();
+
+        assert_eq!(abcd, abcd_with_double_loop);
+    }
+}
+
+fn random_bilinearity_tests<E: Engine>() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..1000 {
+        let a = E::G1::rand(&mut rng);
+        let b = E::G2::rand(&mut rng);
+
+        let c = E::Fr::rand(&mut rng);
+        let d = E::Fr::rand(&mut rng);
+
+        let mut ac = a;
+        ac.mul_assign(c);
+
+        let mut ad = a;
+        ad.mul_assign(d);
+
+        let mut bc = b;
+        bc.mul_assign(c);
+
+        let mut bd = b;
+        bd.mul_assign(d);
+
+        let acbd = E::pairing(ac, bd);
+        let adbc = E::pairing(ad, bc);
+
+        let mut cd = c;
+        cd.mul_assign(&d);
+
+        let abcd = E::pairing(a, b).pow(cd.into_repr());
+
+        assert_eq!(acbd, adbc);
+        assert_eq!(acbd, abcd);
+    }
+}

src/tests/field.rs

@@ -0,0 +1,229 @@
+use rand::{Rng, SeedableRng, XorShiftRng};
+use ::{SqrtField, Field};
+
+pub fn random_frobenius_tests<F: Field, C: AsRef<[u64]>>(characteristic: C, maxpower: usize) {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..100 {
+        for i in 0..(maxpower+1) {
+            let mut a = F::rand(&mut rng);
+            let mut b = a;
+
+            for _ in 0..i {
+                a = a.pow(&characteristic);
+            }
+            b.frobenius_map(i);
+
+            assert_eq!(a, b);
+        }
+    }
+}
+
+pub fn random_sqrt_tests<F: SqrtField>() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    for _ in 0..10000 {
+        let a = F::rand(&mut rng);
+        let mut b = a;
+        b.square();
+
+        let b = b.sqrt().unwrap();
+        let mut negb = b;
+        negb.negate();
+
+        assert!(a == b || a == negb);
+    }
+
+    let mut c = F::one();
+    for _ in 0..10000 {
+        let mut b = c;
+        b.square();
+        b = b.sqrt().unwrap();
+
+        if b != c {
+            b.negate();
+        }
+
+        assert_eq!(b, c);
+
+        c.add_assign(&F::one());
+    }
+}
+
+pub fn random_field_tests<F: Field>() {
+    let mut rng = XorShiftRng::from_seed([0x5dbe6259, 0x8d313d76, 0x3237db17, 0xe5bc0654]);
+
+    random_multiplication_tests::<F, _>(&mut rng);
+    random_addition_tests::<F, _>(&mut rng);
+    random_subtraction_tests::<F, _>(&mut rng);
+    random_negation_tests::<F, _>(&mut rng);
+    random_doubling_tests::<F, _>(&mut rng);
+    random_squaring_tests::<F, _>(&mut rng);
+    random_inversion_tests::<F, _>(&mut rng);
+    random_expansion_tests::<F, _>(&mut rng);
+
+    assert!(F::zero().is_zero());
+    {
+        let mut z = F::zero();
+        z.negate();
+        assert!(z.is_zero());
+    }
+
+    assert!(F::zero().inverse().is_none());
+
+    // Multiplication by zero
+    {
+        let mut a = F::rand(&mut rng);
+        a.mul_assign(&F::zero());
+        assert!(a.is_zero());
+    }
+
+    // Addition by zero
+    {
+        let mut a = F::rand(&mut rng);
+        let copy = a;
+        a.add_assign(&F::zero());
+        assert_eq!(a, copy);
+    }
+}
+
+fn random_multiplication_tests<F: Field, R: Rng>(rng: &mut R) {
+    for _ in 0..10000 {
+        let a = F::rand(rng);
+        let b = F::rand(rng);
+        let c = F::rand(rng);
+
+        let mut t0 = a; // (a * b) * c
+        t0.mul_assign(&b);
+        t0.mul_assign(&c);
+
+        let mut t1 = a; // (a * c) * b
+        t1.mul_assign(&c);
+        t1.mul_assign(&b);
+
+        let mut t2 = b; // (b * c) * a
+        t2.mul_assign(&c);
+        t2.mul_assign(&a);
+
+        assert_eq!(t0, t1);
+        assert_eq!(t1, t2);
+    }
+}
+
+fn random_addition_tests<F: Field, R: Rng>(rng: &mut R) {
+    for _ in 0..10000 {
+        let a = F::rand(rng);
+        let b = F::rand(rng);
+        let c = F::rand(rng);
+
+        let mut t0 = a; // (a + b) + c
+        t0.add_assign(&b);
+        t0.add_assign(&c);
+
+        let mut t1 = a; // (a + c) + b
+        t1.add_assign(&c);
+        t1.add_assign(&b);
+
+        let mut t2 = b; // (b + c) + a
+        t2.add_assign(&c);
+        t2.add_assign(&a);
+
+        assert_eq!(t0, t1);
+        assert_eq!(t1, t2);
+    }
+}
+
+fn random_subtraction_tests<F: Field, R: Rng>(rng: &mut R) {
+    for _ in 0..10000 {
+        let a = F::rand(rng);
+        let b = F::rand(rng);
+
+        let mut t0 = a; // (a - b)
+        t0.sub_assign(&b);
+
+        let mut t1 = b; // (b - a)
+        t1.sub_assign(&a);
+
+        let mut t2 = t0; // (a - b) + (b - a) = 0
+        t2.add_assign(&t1);
+
+        assert!(t2.is_zero());
+    }
+}
+
+fn random_negation_tests<F: Field, R: Rng>(rng: &mut R) {
+    for _ in 0..10000 {
+        let a = F::rand(rng);
+        let mut b = a;
+        b.negate();
+        b.add_assign(&a);
+
+        assert!(b.is_zero());
+    }
+}
+
+fn random_doubling_tests<F: Field, R: Rng>(rng: &mut R) {
+    for _ in 0..10000 {
+        let mut a = F::rand(rng);
+        let mut b = a;
+        a.add_assign(&b);
+        b.double();
+
+        assert_eq!(a, b);
+    }
+}
+
+fn random_squaring_tests<F: Field, R: Rng>(rng: &mut R) {
+    for _ in 0..10000 {
+        let mut a = F::rand(rng);
+        let mut b = a;
+        a.mul_assign(&b);
+        b.square();
+
+        assert_eq!(a, b);
+    }
+}
+
+fn random_inversion_tests<F: Field, R: Rng>(rng: &mut R) {
+    assert!(F::zero().inverse().is_none());
+
+    for _ in 0..10000 {
+        let mut a = F::rand(rng);
+        let b = a.inverse().unwrap(); // probablistically nonzero
+        a.mul_assign(&b);
+
+        assert_eq!(a, F::one());
+    }
+}
+
+fn random_expansion_tests<F: Field, R: Rng>(rng: &mut R) {
+    for _ in 0..10000 {
+        // Compare (a + b)(c + d) and (a*c + b*c + a*d + b*d)
+
+        let a = F::rand(rng);
+        let b = F::rand(rng);
+        let c = F::rand(rng);
+        let d = F::rand(rng);
+
+        let mut t0 = a;
+        t0.add_assign(&b);
+        let mut t1 = c;
+        t1.add_assign(&d);
+        t0.mul_assign(&t1);
+
+        let mut t2 = a;
+        t2.mul_assign(&c);
+        let mut t3 = b;
+        t3.mul_assign(&c);
+        let mut t4 = a;
+        t4.mul_assign(&d);
+        let mut t5 = b;
+        t5.mul_assign(&d);
+
+        t2.add_assign(&t3);
+        t2.add_assign(&t4);
+        t2.add_assign(&t5);
+
+        assert_eq!(t0, t2);
+    }
+}

src/tests/mod.rs

@@ -0,0 +1,3 @@
+pub mod curve;
+pub mod field;
+pub mod engine;