In with the new.

Cargo.toml

@@ -9,4 +9,10 @@ repository = "https://github.com/ebfull/bellman"
 version = "0.0.1"
 
 [dependencies]
-rand = "0.3.15"
+rand = "0.3.*"
+rayon = "0.6.*"
+byteorder = "1.*"
+serde = "0.9.*"
+
+[dev-dependencies.bincode]
+git = "https://github.com/TyOverby/bincode.git"

README.md

@@ -1,7 +1,3 @@
 # bellman [![Build status](https://api.travis-ci.org/ebfull/bellman.svg)](https://travis-ci.org/ebfull/bellman) [![Crates.io](https://img.shields.io/crates/v/bellman.svg)](https://crates.io/crates/bellman) #
 
-This is a Rust language zk-SNARK crate. (You should not use this in production for anything yet.)
-
-------
-
 This is a research project being built for [Zcash](https://z.cash/).

benches/curve.rs

@@ -0,0 +1,235 @@
+#![feature(test)]
+#![allow(unused_mut)]
+
+extern crate rand;
+extern crate bellman;
+extern crate test;
+
+use bellman::curves::*;
+use bellman::curves::bls381::*;
+
+const SAMPLES: usize = 30;
+
+macro_rules! benchmark(
+    ($name:ident, $engine:ident, $input:ident($rng:ident) = $pre:expr; $post:expr) => (
+        #[bench]
+        fn $name(b: &mut test::Bencher) {
+            let $rng = &mut rand::thread_rng();
+            let $engine = &Bls381::new();
+            let $input: Vec<_> = (0..SAMPLES).map(|_| $pre).collect();
+
+            let mut c = 0;
+
+            b.iter(|| {
+                c += 1;
+
+                let mut $input = $input[c % SAMPLES].clone();
+
+                $post
+            })
+        }
+    )
+);
+
+benchmark!(g1_multiexp, e,
+    input(rng) = {
+        let mut g = G1::random(e, rng);
+        let mut a = G1::random(e, rng);
+        (
+            (0..1000).map(|_| {
+                g.add_assign(e, &a);
+                a.double(e);
+                g.to_affine(e)
+            }).collect::<Vec<_>>(),
+            (0..1000).map(|_| Fr::random(e, rng)).collect::<Vec<_>>(),
+        )
+    };
+
+    e.multiexp::<G1>(&input.0, &input.1)
+);
+
+benchmark!(g2_multiexp, e,
+    input(rng) = {
+        let mut g = G2::random(e, rng);
+        let mut a = G2::random(e, rng);
+        (
+            (0..1000).map(|_| {
+                g.add_assign(e, &a);
+                a.double(e);
+                g.to_affine(e)
+            }).collect::<Vec<_>>(),
+            (0..1000).map(|_| Fr::random(e, rng)).collect::<Vec<_>>(),
+        )
+    };
+
+    e.multiexp::<G2>(&input.0, &input.1)
+);
+
+benchmark!(full_pairing, e,
+    input(rng) = (G1::random(e, rng), G2::random(e, rng));
+
+    e.pairing(&input.0, &input.1)
+);
+
+benchmark!(g1_pairing_preparation, e,
+    input(rng) = G1::random(e, rng);
+
+    input.prepare(e)
+);
+
+benchmark!(g2_pairing_preparation, e,
+    input(rng) = G2::random(e, rng);
+
+    input.prepare(e)
+);
+
+benchmark!(miller_loop, e,
+    input(rng) = (G1::random(e, rng).prepare(e), G2::random(e, rng).prepare(e));
+
+    e.miller_loop([(&input.0, &input.1)].into_iter())
+);
+
+benchmark!(double_miller_loop, e,
+    input(rng) = (G1::random(e, rng).prepare(e), G2::random(e, rng).prepare(e), G1::random(e, rng).prepare(e), G2::random(e, rng).prepare(e));
+
+    e.miller_loop([
+        (&input.0, &input.1),
+        (&input.2, &input.3),
+    ].into_iter())
+);
+
+benchmark!(final_exponentiation, e,
+    input(rng) = e.miller_loop([
+        (&G1::random(e, rng).prepare(e), &G2::random(e, rng).prepare(e)),
+    ].into_iter());
+
+    e.final_exponentiation(&input)
+);
+
+macro_rules! group_tests(
+    (
+        $name:ident,
+        $mul:ident,
+        $mul_mixed:ident,
+        $add:ident
+    ) => {
+        benchmark!($mul, e,
+            input(rng) = ($name::random(e, rng), Fr::random(e, rng));
+
+            {input.0.mul_assign(e, &input.1); input.0}
+        );
+
+        benchmark!($mul_mixed, e,
+            input(rng) = ($name::random(e, rng).to_affine(e), Fr::random(e, rng));
+
+            {input.0.mul(e, &input.1)}
+        );
+
+        benchmark!($add, e,
+            input(rng) = ($name::random(e, rng), $name::random(e, rng));
+
+            {input.0.add_assign(e, &input.1); input.0}
+        );
+    };
+);
+
+macro_rules! field_tests(
+    (
+        @nosqrt,
+        $name:ident,
+        $mul:ident,
+        $square:ident,
+        $add:ident,
+        $inverse:ident
+    ) => {
+        benchmark!($mul, e,
+            input(rng) = ($name::random(e, rng), $name::random(e, rng));
+
+            {input.0.mul_assign(e, &input.1); input.0}
+        );
+
+        benchmark!($square, e,
+            input(rng) = $name::random(e, rng);
+
+            {input.square(e); input}
+        );
+
+        benchmark!($add, e,
+            input(rng) = ($name::random(e, rng), $name::random(e, rng));
+
+            {input.0.add_assign(e, &input.1); input.0}
+        );
+
+        benchmark!($inverse, e,
+            input(rng) = $name::random(e, rng);
+
+            {input.inverse(e)}
+        );
+    };
+    (
+        $name:ident,
+        $mul:ident,
+        $square:ident,
+        $add:ident,
+        $inverse:ident,
+        $sqrt:ident
+    ) => {
+        field_tests!(@nosqrt, $name, $mul, $square, $add, $inverse);
+
+        benchmark!($sqrt, e,
+            input(rng) = {let mut tmp = $name::random(e, rng); tmp.square(e); tmp};
+
+            {input.sqrt(e)}
+        );
+    };
+);
+
+field_tests!(
+    Fr,
+    fr_multiplication,
+    fr_squaring,
+    fr_addition,
+    fr_inverse,
+    fr_sqrt
+);
+
+field_tests!(
+    Fq,
+    fq_multiplication,
+    fq_squaring,
+    fq_addition,
+    fq_inverse,
+    fq_sqrt
+);
+
+field_tests!(
+    Fq2,
+    fq2_multiplication,
+    fq2_squaring,
+    fq2_addition,
+    fq2_inverse,
+    fq2_sqrt
+);
+
+field_tests!(
+    @nosqrt,
+    Fq12,
+    fq12_multiplication,
+    fq12_squaring,
+    fq12_addition,
+    fq12_inverse
+);
+
+group_tests!(
+    G1,
+    g1_multiplication,
+    g1_multiplication_mixed,
+    g1_addition
+);
+
+group_tests!(
+    G2,
+    g2_multiplication,
+    g2_multiplication_mixed,
+    g2_addition
+);

src/curves/bls381/ec.rs

@@ -0,0 +1,398 @@
+macro_rules! curve_impl {
+    (
+        $engine:ident,
+        $name:ident,
+        $name_affine:ident,
+        $name_prepared:ident,
+        $name_uncompressed:ident,
+        $params_name:ident,
+        $params_field:ident,
+        $basefield:ident,
+        $scalarfield:ident
+    ) => {
+        #[repr(C)]
+        #[derive(Copy, Clone, PartialEq, Eq, Debug)]
+        pub struct $name_affine {
+            x: $basefield,
+            y: $basefield,
+            infinity: bool
+        }
+
+        #[repr(C)]
+        #[derive(Copy, Clone, Debug)]
+        pub struct $name {
+            x: $basefield,
+            y: $basefield,
+            z: $basefield
+        }
+
+        struct $params_name {
+            zero: $name,
+            one: $name,
+            coeff_b: $basefield,
+            windows: Vec<usize>
+        }
+
+        impl Convert<$name_affine, $engine> for $name {
+            type Target = $name_affine;
+
+            fn convert(&self, engine: &$engine) -> Cow<$name_affine> {
+                Cow::Owned(self.to_affine(engine))
+            }
+        }
+
+        impl GroupAffine<$engine, $name> for $name_affine {
+            type Uncompressed = $name_uncompressed;
+
+            fn is_valid(&self, e: &$engine) -> bool {
+                if self.is_zero() {
+                    true
+                } else {
+                    // Check that the point is on the curve
+                    let mut y2 = self.y;
+                    y2.square(e);
+
+                    let mut x3b = self.x;
+                    x3b.square(e);
+                    x3b.mul_assign(e, &self.x);
+                    x3b.add_assign(e, &e.$params_field.coeff_b);
+
+                    if y2 == x3b {
+                        // Check that the point is in the correct subgroup
+                        if self.mul(e, &$scalarfield::char(e)).is_zero() {
+                            true
+                        } else {
+                            false
+                        }
+                    } else {
+                        false
+                    }
+                }
+            }
+
+            fn to_uncompressed(&self, engine: &$engine) -> Self::Uncompressed {
+                $name_uncompressed::from_affine(self, engine)
+            }
+
+            fn to_jacobian(&self, engine: &$engine) -> $name {
+                if self.infinity {
+                    $name::zero(engine)
+                } else {
+                    $name {
+                        x: self.x,
+                        y: self.y,
+                        z: $basefield::one(engine)
+                    }
+                }
+            }
+
+            fn prepare(self, e: &$engine) -> $name_prepared {
+                $name_prepared::from_engine(e, self)
+            }
+
+            fn is_zero(&self) -> bool {
+                self.infinity
+            }
+
+            fn mul<S: Convert<[u64], $engine>>(&self, e: &$engine, other: &S) -> $name {
+                let mut res = $name::zero(e);
+
+                for i in BitIterator::from((*other.convert(e)).borrow())
+                {
+                    res.double(e);
+
+                    if i {
+                        res.add_assign_mixed(e, self);
+                    }
+                }
+
+                res
+            }
+
+            fn negate(&mut self, e: &$engine) {
+                if !self.is_zero() {
+                    self.y.negate(e);
+                }
+            }
+        }
+
+        impl Group<$engine> for $name {
+            type Affine = $name_affine;
+            type Prepared = $name_prepared;
+
+            fn optimal_window(engine: &$engine, scalar_bits: usize) -> Option<usize> {
+                for (i, bits) in engine.$params_field.windows.iter().enumerate().rev() {
+                    if &scalar_bits >= bits {
+                        return Some(i + 2);
+                    }
+                }
+
+                None
+            }
+
+            fn zero(engine: &$engine) -> Self {
+                engine.$params_field.zero
+            }
+
+            fn one(engine: &$engine) -> Self {
+                engine.$params_field.one
+            }
+
+            fn random<R: rand::Rng>(engine: &$engine, rng: &mut R) -> Self {
+                let mut tmp = Self::one(engine);
+                tmp.mul_assign(engine, &$scalarfield::random(engine, rng));
+
+                tmp
+            }
+
+            fn is_zero(&self) -> bool {
+                self.z.is_zero()
+            }
+
+            fn is_equal(&self, engine: &$engine, other: &Self) -> bool {
+                if self.is_zero() {
+                    return other.is_zero();
+                }
+
+                if other.is_zero() {
+                    return false;
+                }
+
+                let mut z1 = self.z;
+                z1.square(engine);
+                let mut z2 = other.z;
+                z2.square(engine);
+
+                let mut tmp1 = self.x;
+                tmp1.mul_assign(engine, &z2);
+
+                let mut tmp2 = other.x;
+                tmp2.mul_assign(engine, &z1);
+
+                if tmp1 != tmp2 {
+                    return false;
+                }
+
+                z1.mul_assign(engine, &self.z);
+                z2.mul_assign(engine, &other.z);
+                z2.mul_assign(engine, &self.y);
+                z1.mul_assign(engine, &other.y);
+
+                if z1 != z2 {
+                    return false;
+                }
+
+                true
+            }
+
+            fn to_affine(&self, engine: &$engine) -> Self::Affine {
+                if self.is_zero() {
+                    $name_affine {
+                        x: $basefield::zero(),
+                        y: $basefield::one(engine),
+                        infinity: true
+                    }
+                } else {
+                    let zinv = self.z.inverse(engine).unwrap();
+                    let mut zinv_powered = zinv;
+                    zinv_powered.square(engine);
+
+                    let mut x = self.x;
+                    x.mul_assign(engine, &zinv_powered);
+
+                    let mut y = self.y;
+                    zinv_powered.mul_assign(engine, &zinv);
+                    y.mul_assign(engine, &zinv_powered);
+
+                    $name_affine {
+                        x: x,
+                        y: y,
+                        infinity: false
+                    }
+                }
+            }
+
+            fn prepare(&self, e: &$engine) -> $name_prepared {
+                self.to_affine(e).prepare(e)
+            }
+
+            fn double(&mut self, engine: &$engine) {
+                if self.is_zero() {
+                    return;
+                }
+
+                let mut a = self.x;
+                a.square(engine);
+                let mut c = self.y;
+                c.square(engine);
+                let mut d = c;
+                c.square(engine);
+                d.add_assign(engine, &self.x);
+                d.square(engine);
+                d.sub_assign(engine, &a);
+                d.sub_assign(engine, &c);
+                d.double(engine);
+                let mut e = a;
+                e.add_assign(engine, &a);
+                e.add_assign(engine, &a);
+                self.x = e;
+                self.x.square(engine);
+                self.x.sub_assign(engine, &d);
+                self.x.sub_assign(engine, &d);
+                c.double(engine);
+                c.double(engine);
+                c.double(engine);
+                self.z.mul_assign(engine, &self.y);
+                self.z.double(engine);
+                self.y = d;
+                self.y.sub_assign(engine, &self.x);
+                self.y.mul_assign(engine, &e);
+                self.y.sub_assign(engine, &c);
+            }
+
+            fn negate(&mut self, engine: &$engine) {
+                if !self.is_zero() {
+                    self.y.negate(engine)
+                }
+            }
+
+            fn mul_assign<S: Convert<[u64], $engine>>(&mut self, engine: &$engine, other: &S) {
+                let mut res = Self::zero(engine);
+
+                for i in BitIterator::from((*other.convert(engine)).borrow())
+                {
+                    res.double(engine);
+
+                    if i {
+                        res.add_assign(engine, self);
+                    }
+                }
+
+                *self = res;
+            }
+
+            fn sub_assign(&mut self, engine: &$engine, other: &Self) {
+                let mut tmp = *other;
+                tmp.negate(engine);
+                self.add_assign(engine, &tmp);
+            }
+
+            fn add_assign_mixed(&mut self, e: &$engine, other: &$name_affine) {
+                if other.is_zero() {
+                    return;
+                }
+
+                if self.is_zero() {
+                    self.x = other.x;
+                    self.y = other.y;
+                    self.z = $basefield::one(e);
+                    return;
+                }
+                
+                let mut z1z1 = self.z;
+                z1z1.square(e);
+                let mut u2 = other.x;
+                u2.mul_assign(e, &z1z1);
+                let mut z1cubed = self.z;
+                z1cubed.mul_assign(e, &z1z1);
+                let mut s2 = other.y;
+                s2.mul_assign(e, &z1cubed);
+
+                if self.x == u2 && self.y == s2 {
+                    self.double(e);
+                    return;
+                }
+
+                let mut h = u2;
+                h.sub_assign(e, &self.x);
+                let mut hh = h;
+                hh.square(e);
+                let mut i = hh;
+                i.double(e);
+                i.double(e);
+                let mut j = h;
+                j.mul_assign(e, &i);
+                let mut r = s2;
+                r.sub_assign(e, &self.y);
+                r.double(e);
+                let mut v = self.x;
+                v.mul_assign(e, &i);
+
+                self.x = r;
+                self.x.square(e);
+                self.x.sub_assign(e, &j);
+                self.x.sub_assign(e, &v);
+                self.x.sub_assign(e, &v);
+
+                self.y.mul_assign(e, &j);
+                let mut tmp = v;
+                tmp.sub_assign(e, &self.x);
+                tmp.mul_assign(e, &r);
+                tmp.sub_assign(e, &self.y);
+                tmp.sub_assign(e, &self.y);
+                self.y = tmp;
+
+                self.z.add_assign(e, &h);
+                self.z.square(e);
+                self.z.sub_assign(e, &z1z1);
+                self.z.sub_assign(e, &hh);
+            }
+
+            fn add_assign(&mut self, engine: &$engine, other: &Self) {
+                if self.is_zero() {
+                    *self = *other;
+                    return;
+                }
+
+                if other.is_zero() {
+                    return;
+                }
+
+                let mut z1_squared = self.z;
+                z1_squared.square(engine);
+                let mut z2_squared = other.z;
+                z2_squared.square(engine);
+                let mut u1 = self.x;
+                u1.mul_assign(engine, &z2_squared);
+                let mut u2 = other.x;
+                u2.mul_assign(engine, &z1_squared);
+                let mut s1 = other.z;
+                s1.mul_assign(engine, &z2_squared);
+                s1.mul_assign(engine, &self.y);
+                let mut s2 = self.z;
+                s2.mul_assign(engine, &z1_squared);
+                s2.mul_assign(engine, &other.y);
+
+                if u1 == u2 && s1 == s2 {
+                    self.double(engine);
+                } else {
+                    u2.sub_assign(engine, &u1);
+                    s2.sub_assign(engine, &s1);
+                    s2.double(engine);
+                    let mut i = u2;
+                    i.double(engine);
+                    i.square(engine);
+                    let mut v = i;
+                    v.mul_assign(engine, &u1);
+                    i.mul_assign(engine, &u2);
+                    s1.mul_assign(engine, &i);
+                    s1.double(engine);
+                    self.x = s2;
+                    self.x.square(engine);
+                    self.x.sub_assign(engine, &i);
+                    self.x.sub_assign(engine, &v);
+                    self.x.sub_assign(engine, &v);
+                    self.y = v;
+                    self.y.sub_assign(engine, &self.x);
+                    self.y.mul_assign(engine, &s2);
+                    self.y.sub_assign(engine, &s1);
+                    self.z.add_assign(engine, &other.z);
+                    self.z.square(engine);
+                    self.z.sub_assign(engine, &z1_squared);
+                    self.z.sub_assign(engine, &z2_squared);
+                    self.z.mul_assign(engine, &u2);
+                }
+            }
+        }
+    }
+}

src/curves/bls381/fp.rs

@@ -0,0 +1,714 @@
+macro_rules! fp_params_impl {
+    (
+        $name:ident = (3 mod 4),
+        engine = $engine:ident,
+        params = $params_field:ident : $params_name:ident,
+        limbs = $limbs:expr,
+        modulus = $modulus:expr,
+        r1 = $r1:expr,
+        r2 = $r2:expr,
+        modulus_minus_3_over_4 = $modulus_minus_3_over_4:expr,
+        modulus_minus_1_over_2 = $modulus_minus_1_over_2:expr,
+        inv = $inv:expr
+    ) => {
+        struct $params_name {
+            modulus: [u64; $limbs],
+            r1: $name,
+            r2: $name,
+            inv: u64,
+            one: $name,
+            num_bits: usize,
+            modulus_minus_3_over_4: [u64; $limbs],
+            modulus_minus_1_over_2: [u64; $limbs],
+            base10: [$name; 11]
+        }
+
+        impl $params_name {
+            fn partial_init() -> $params_name {
+                let mut tmp = $params_name {
+                    modulus: $modulus,
+                    r1: $name($r1),
+                    r2: $name($r2),
+                    inv: $inv,
+                    one: $name::zero(),
+                    num_bits: 0,
+                    modulus_minus_3_over_4: $modulus_minus_3_over_4,
+                    modulus_minus_1_over_2: $modulus_minus_1_over_2,
+                    base10: [$name::zero(); 11]
+                };
+
+                tmp.one.0[0] = 1;
+
+                tmp
+            }
+        }
+    };
+    (
+        $name:ident = (1 mod 16),
+        engine = $engine:ident,
+        params = $params_field:ident : $params_name:ident,
+        limbs = $limbs:expr,
+        modulus = $modulus:expr,
+        r1 = $r1:expr,
+        r2 = $r2:expr,
+        modulus_minus_1_over_2 = $modulus_minus_1_over_2:expr,
+        s = $s:expr,
+        t = $t:expr,
+        t_plus_1_over_2 = $t_plus_1_over_2:expr,
+        inv = $inv:expr
+    ) => {
+        struct $params_name {
+            modulus: [u64; $limbs],
+            r1: $name,
+            r2: $name,
+            inv: u64,
+            one: $name,
+            num_bits: usize,
+            modulus_minus_1_over_2: [u64; $limbs],
+            s: u64,
+            t: [u64; $limbs],
+            t_plus_1_over_2: [u64; $limbs],
+            root_of_unity: $name,
+            multiplicative_generator: $name,
+            base10: [$name; 11]
+        }
+
+        impl $params_name {
+            fn partial_init() -> $params_name {
+                let mut tmp = $params_name {
+                    modulus: $modulus,
+                    r1: $name($r1),
+                    r2: $name($r2),
+                    inv: $inv,
+                    one: $name::zero(),
+                    num_bits: 0,
+                    modulus_minus_1_over_2: $modulus_minus_1_over_2,
+                    s: $s,
+                    t: $t,
+                    t_plus_1_over_2: $t_plus_1_over_2,
+                    root_of_unity: $name::zero(),
+                    multiplicative_generator: $name::zero(),
+                    base10: [$name::zero(); 11]
+                };
+
+                tmp.one.0[0] = 1;
+
+                tmp
+            }
+        }
+    };
+}
+
+macro_rules! fp_sqrt_impl {
+    (
+        $name:ident = (3 mod 4),
+        engine = $engine:ident,
+        params = $params_field:ident : $params_name:ident
+    ) => {
+        impl SqrtField<$engine> for $name {
+            fn sqrt(&self, engine: &$engine) -> Option<Self> {
+                let mut a1 = self.pow(engine, &engine.$params_field.modulus_minus_3_over_4);
+
+                let mut a0 = a1;
+                a0.square(engine);
+                a0.mul_assign(engine, self);
+
+                let mut neg1 = Self::one(engine);
+                neg1.negate(engine);
+
+                if a0 == neg1 {
+                    None
+                } else {
+                    a1.mul_assign(engine, self);
+                    Some(a1)
+                }
+            }
+        }
+    };
+    (
+        $name:ident = (1 mod 16),
+        engine = $engine:ident,
+        params = $params_field:ident : $params_name:ident
+    ) => {
+        impl SqrtField<$engine> for $name {
+            fn sqrt(&self, engine: &$engine) -> Option<Self> {
+                if self.is_zero() {
+                    return Some(*self);
+                }
+
+                if self.pow(engine, &engine.$params_field.modulus_minus_1_over_2) != $name::one(engine) {
+                    None
+                } else {
+                    let mut c = engine.$params_field.root_of_unity;
+
+                    let mut r = self.pow(engine, &engine.$params_field.t_plus_1_over_2);
+                    let mut t = self.pow(engine, &engine.$params_field.t);
+                    let mut m = engine.$params_field.s;
+                    
+                    while t != Self::one(engine) {
+                        let mut i = 1;
+                        {
+                            let mut t2i = t;
+                            t2i.square(engine);
+                            loop {
+                                if t2i == Self::one(engine) {
+                                    break;
+                                }
+                                t2i.square(engine);
+                                i += 1;
+                            }
+                        }
+
+                        for _ in 0..(m - i - 1) {
+                            c.square(engine);
+                        }
+                        r.mul_assign(engine, &c);
+                        c.square(engine);
+                        t.mul_assign(engine, &c);
+                        m = i;
+                    }
+
+                    Some(r)
+                }
+            }
+        }
+    };
+}
+
+macro_rules! fp_impl {
+    (
+        $name:ident = ($($congruency:tt)*),
+        engine = $engine:ident,
+        params = $params_field:ident : $params_name:ident,
+        arith = $arith_mod:ident,
+        limbs = $limbs:expr,
+        $($params:tt)*
+    ) => {
+        fp_params_impl!(
+            $name = ($($congruency)*),
+            engine = $engine,
+            params = $params_field : $params_name,
+            limbs = $limbs,
+            $($params)*
+        );
+
+        impl $params_name {
+            fn base10(e: &$engine) -> [$name; 11] {
+                let mut ret = [$name::zero(); 11];
+
+                let mut acc = $name::zero();
+                for i in 0..11 {
+                    ret[i] = acc;
+                    acc.add_assign(e, &$name::one(e));
+                }
+
+                ret
+            }
+        }
+
+        fp_sqrt_impl!(
+            $name = ($($congruency)*),
+            engine = $engine,
+            params = $params_field : $params_name
+        );
+
+        #[derive(Copy, Clone, PartialEq, Eq)]
+        #[repr(C)]
+        pub struct $name([u64; $limbs]);
+
+        impl fmt::Debug for $name
+        {
+            fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+                ENGINE.with(|e| {
+                    let mut repr = self.into_repr(&e);
+                    repr.reverse();
+
+                    try!(write!(f, "Fp(0x"));
+                    for i in &repr {
+                        try!(write!(f, "{:016x}", *i));
+                    }
+                    write!(f, ")")
+                })
+            }
+        }
+
+        impl $name
+        {
+            #[inline]
+            fn mont_reduce(&mut self, engine: &$engine, res: &mut [u64; $limbs*2]) {
+                // The Montgomery reduction here is based on Algorithm 14.32 in
+                // Handbook of Applied Cryptography
+                // <http://cacr.uwaterloo.ca/hac/about/chap14.pdf>.
+
+                for i in 0..$limbs {
+                    let k = res[i].wrapping_mul(engine.$params_field.inv);
+                    $arith_mod::mac_digit(&mut res[i..], &engine.$params_field.modulus, k);
+                }
+
+                self.0.copy_from_slice(&res[$limbs..]);
+
+                self.reduce(engine);
+            }
+
+            #[inline]
+            fn reduce(&mut self, engine: &$engine) {
+                if !$arith_mod::lt(&self.0, &engine.$params_field.modulus) {
+                    $arith_mod::sub_noborrow(&mut self.0, &engine.$params_field.modulus);
+                }
+            }
+        }
+
+        impl Convert<[u64], $engine> for $name
+        {
+            type Target = [u64; $limbs];
+
+            fn convert(&self, engine: &$engine) -> Cow<[u64; $limbs]> {
+                Cow::Owned(self.into_repr(engine))
+            }
+        }
+
+        impl PrimeField<$engine> for $name
+        {
+            type Repr = [u64; $limbs];
+
+            fn from_repr(engine: &$engine, repr: Self::Repr) -> Result<Self, ()> {
+                let mut tmp = $name(repr);
+                if $arith_mod::lt(&tmp.0, &engine.$params_field.modulus) {
+                    tmp.mul_assign(engine, &engine.$params_field.r2);
+                    Ok(tmp)
+                } else {
+                    Err(())
+                }
+            }
+
+            fn into_repr(&self, engine: &$engine) -> Self::Repr {
+                let mut tmp = *self;
+                tmp.mul_assign(engine, &engine.$params_field.one);
+                tmp.0
+            }
+
+            fn from_u64(engine: &$engine, n: u64) -> Self {
+                let mut r = [0; $limbs];
+                r[0] = n;
+
+                Self::from_repr(engine, r).unwrap()
+            }
+
+            fn from_str(engine: &$engine, s: &str) -> Result<Self, ()> {
+                let mut res = Self::zero();
+                for c in s.chars() {
+                    match c.to_digit(10) {
+                        Some(d) => {
+                            res.mul_assign(engine, &engine.$params_field.base10[10]);
+                            res.add_assign(engine, &engine.$params_field.base10[d as usize]);
+                        },
+                        None => {
+                            return Err(());
+                        }
+                    }
+                }
+
+                Ok(res)
+            }
+
+            fn bits(&self, engine: &$engine) -> BitIterator<Self::Repr> {
+                self.into_repr(engine).into()
+            }
+
+            fn char(engine: &$engine) -> Self::Repr {
+                engine.$params_field.modulus
+            }
+
+            fn num_bits(engine: &$engine) -> usize {
+                engine.$params_field.num_bits
+            }
+
+            fn capacity(engine: &$engine) -> usize {
+                Self::num_bits(engine) - 1
+            }
+        }
+
+        impl Field<$engine> for $name
+        {
+            fn zero() -> Self {
+                $name([0; $limbs])
+            }
+
+            fn one(engine: &$engine) -> Self {
+                engine.$params_field.r1
+            }
+
+            fn random<R: rand::Rng>(engine: &$engine, rng: &mut R) -> Self {
+                let mut tmp = [0; $limbs*2];
+                for i in &mut tmp {
+                    *i = rng.gen();
+                }
+
+                $name($arith_mod::divrem(&tmp, &engine.$params_field.modulus).1)
+            }
+
+            fn is_zero(&self) -> bool {
+                self.0.iter().all(|&e| e==0)
+            }
+
+            fn double(&mut self, engine: &$engine) {
+                $arith_mod::mul2(&mut self.0);
+
+                self.reduce(engine);
+            }
+
+            fn frobenius_map(&mut self, _: &$engine, _: usize)
+            {
+                // This is the identity function for a prime field.
+            }
+
+            fn negate(&mut self, engine: &$engine) {
+                if !self.is_zero() {
+                    let mut tmp = engine.$params_field.modulus;
+                    $arith_mod::sub_noborrow(&mut tmp, &self.0);
+                    
+                    self.0 = tmp;
+                }
+            }
+
+            fn add_assign(&mut self, engine: &$engine, other: &Self) {
+                $arith_mod::add_nocarry(&mut self.0, &other.0);
+
+                self.reduce(engine);
+            }
+
+            fn sub_assign(&mut self, engine: &$engine, other: &Self) {
+                if $arith_mod::lt(&self.0, &other.0) {
+                    $arith_mod::add_nocarry(&mut self.0, &engine.$params_field.modulus);
+                }
+
+                $arith_mod::sub_noborrow(&mut self.0, &other.0);
+            }
+
+            fn square(&mut self, engine: &$engine)
+            {
+                let mut res = [0; $limbs*2];
+                $arith_mod::mac3(&mut res, &self.0, &self.0);
+
+                self.mont_reduce(engine, &mut res);
+            }
+
+            fn mul_assign(&mut self, engine: &$engine, other: &Self) {
+                let mut res = [0; $limbs*2];
+                $arith_mod::mac3(&mut res, &self.0, &other.0);
+
+                self.mont_reduce(engine, &mut res);
+            }
+
+            fn inverse(&self, engine: &$engine) -> Option<Self> {
+                if self.is_zero() {
+                    None
+                } else {
+                    // Guajardo Kumar Paar Pelzl
+                    // Efficient Software-Implementation of Finite Fields with Applications to Cryptography
+                    // Algorithm 16 (BEA for Inversion in Fp)
+
+                    let mut u = self.0;
+                    let mut v = engine.$params_field.modulus;
+                    let mut b = engine.$params_field.r2; // Avoids unnecessary reduction step.
+                    let mut c = Self::zero();
+
+                    while u != engine.$params_field.one.0 && v != engine.$params_field.one.0 {
+                        while $arith_mod::even(&u) {
+                            $arith_mod::div2(&mut u);
+
+                            if $arith_mod::even(&b.0) {
+                                $arith_mod::div2(&mut b.0);
+                            } else {
+                                $arith_mod::add_nocarry(&mut b.0, &engine.$params_field.modulus);
+                                $arith_mod::div2(&mut b.0);
+                            }
+                        }
+
+                        while $arith_mod::even(&v) {
+                            $arith_mod::div2(&mut v);
+
+                            if $arith_mod::even(&c.0) {
+                                $arith_mod::div2(&mut c.0);
+                            } else {
+                                $arith_mod::add_nocarry(&mut c.0, &engine.$params_field.modulus);
+                                $arith_mod::div2(&mut c.0);
+                            }
+                        }
+
+                        if $arith_mod::lt(&v, &u) {
+                            $arith_mod::sub_noborrow(&mut u, &v);
+                            b.sub_assign(engine, &c);
+                        } else {
+                            $arith_mod::sub_noborrow(&mut v, &u);
+                            c.sub_assign(engine, &b);
+                        }
+                    }
+
+                    if u == engine.$params_field.one.0 {
+                        Some(b)
+                    } else {
+                        Some(c)
+                    }
+                }
+            }
+        }
+
+        mod $arith_mod {
+            use super::BitIterator;
+            // Arithmetic
+            #[allow(dead_code)]
+            pub fn num_bits(v: &[u64; $limbs]) -> usize
+            {
+                // TODO: optimize
+                for (i, b) in BitIterator::from(&v[..]).enumerate() {
+                    if b {
+                        return ($limbs*64) - i;
+                    }
+                }
+
+                0
+            }
+
+            #[inline]
+            pub fn mac_digit(acc: &mut [u64], b: &[u64], c: u64)
+            {
+                #[inline]
+                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
+                }
+
+                let mut b_iter = b.iter();
+                let mut carry = 0;
+
+                for ai in acc.iter_mut() {
+                    if let Some(bi) = b_iter.next() {
+                        *ai = mac_with_carry(*ai, *bi, c, &mut carry);
+                    } else {
+                        *ai = mac_with_carry(*ai, 0, c, &mut carry);
+                    }
+                }
+
+                debug_assert!(carry == 0);
+            }
+
+            #[inline]
+            pub fn mac3_long(acc: &mut [u64], b: &[u64], c: &[u64]) {
+                for (i, xi) in b.iter().enumerate() {
+                    mac_digit(&mut acc[i..], c, *xi);
+                }
+            }
+
+            #[inline]
+            pub fn mac3(acc: &mut [u64; $limbs*2], b: &[u64; $limbs], c: &[u64; $limbs]) {
+                if $limbs > 4 {
+                    let (x0, x1) = b.split_at($limbs / 2);
+                    let (y0, y1) = c.split_at($limbs / 2);
+
+                    let mut p = [0; $limbs+1];
+                    mac3_long(&mut p, x1, y1);
+                    add_nocarry(&mut acc[$limbs/2..], &p);
+                    add_nocarry(&mut acc[$limbs..], &p);
+                    p = [0; $limbs+1];
+                    mac3_long(&mut p, x0, y0);
+                    add_nocarry(&mut acc[..], &p);
+                    add_nocarry(&mut acc[$limbs/2..], &p);
+
+                    let mut sign;
+                    let mut j0 = [0; $limbs / 2];
+                    let mut j1 = [0; $limbs / 2];
+                    if lt(x1, x0) {
+                        sign = false;
+                        j0.copy_from_slice(x0);
+                        sub_noborrow(&mut j0, x1);
+                    } else {
+                        sign = true;
+                        j0.copy_from_slice(x1);
+                        sub_noborrow(&mut j0, x0);
+                    }
+                    if lt(&y1, &y0) {
+                        sign = sign == false;
+                        j1.copy_from_slice(y0);
+                        sub_noborrow(&mut j1, y1);
+                    } else {
+                        sign = sign == true;
+                        j1.copy_from_slice(y1);
+                        sub_noborrow(&mut j1, y0);
+                    }
+
+                    if sign {
+                        p = [0; $limbs+1];
+                        mac3_long(&mut p, &j0, &j1);
+                        sub_noborrow(&mut acc[$limbs/2..], &p);
+                    } else {
+                        mac3_long(&mut acc[$limbs/2..], &j0, &j1);
+                    }
+                } else {
+                    mac3_long(acc, b, c);
+                }
+            }
+
+            #[inline]
+            pub 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
+            }
+
+            #[inline]
+            #[allow(dead_code)]
+            pub fn add_carry(a: &mut [u64], b: &[u64]) {
+                use std::iter;
+
+                let mut carry = 0;
+
+                for (a, b) in a.into_iter().zip(b.iter().chain(iter::repeat(&0))) {
+                    *a = adc(*a, *b, &mut carry);
+                }
+
+                debug_assert!(0 == carry);
+            }
+
+            #[inline]
+            pub fn add_nocarry(a: &mut [u64], b: &[u64]) {
+                let mut carry = 0;
+
+                for (a, b) in a.into_iter().zip(b.iter()) {
+                    *a = adc(*a, *b, &mut carry);
+                }
+
+                debug_assert!(0 == carry);
+            }
+
+            /// Returns true if a < b.
+            #[inline]
+            pub fn lt(a: &[u64], b: &[u64]) -> bool {
+                for (a, b) in a.iter().zip(b.iter()).rev() {
+                    if *a > *b {
+                        return false;
+                    } else if *a < *b {
+                        return true;
+                    }
+                }
+
+                false
+            }
+
+            #[inline]
+            pub fn sub_noborrow(a: &mut [u64], b: &[u64]) {
+                #[inline]
+                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
+                }
+
+                let mut borrow = 0;
+
+                for (a, b) in a.into_iter().zip(b.iter()) {
+                    *a = sbb(*a, *b, &mut borrow);
+                }
+
+                debug_assert!(0 == borrow);
+            }
+
+            /// Returns if number is even.
+            #[inline(always)]
+            #[allow(dead_code)]
+            pub fn even(a: &[u64; $limbs]) -> bool {
+                a[0] & 1 == 0
+            }
+
+            #[inline(always)]
+            #[allow(dead_code)]
+            pub fn odd(a: &[u64; $limbs]) -> bool {
+                a[0] & 1 == 1
+            }
+
+            /// Divide by two
+            #[inline]
+            pub fn div2(a: &mut [u64; $limbs]) {
+                let mut t = 0;
+                for i in a.iter_mut().rev() {
+                    let t2 = *i << 63;
+                    *i >>= 1;
+                    *i |= t;
+                    t = t2;
+                }
+            }
+
+            #[inline]
+            pub fn mul2(a: &mut [u64; $limbs]) {
+                let mut last = 0;
+                for i in a {
+                    let tmp = *i >> 63;
+                    *i <<= 1;
+                    *i |= last;
+                    last = tmp;
+                }
+            }
+
+            fn get_bit(this: &[u64; $limbs*2], n: usize) -> bool {
+                let part = n / 64;
+                let bit = n - (64 * part);
+
+                this[part] & (1 << bit) > 0
+            }
+
+            fn set_bit(this: &mut [u64; $limbs], n: usize, to: bool) -> bool
+            {
+                let part = n / 64;
+                let bit = n - (64 * part);
+
+                match this.get_mut(part) {
+                    Some(e) => {
+                        if to {
+                            *e |= 1 << bit;
+                        } else {
+                            *e &= !(1 << bit);
+                        }
+
+                        true
+                    },
+                    None => false
+                }
+            }
+
+            pub fn divrem(
+                this: &[u64; $limbs*2],
+                modulo: &[u64; $limbs]
+            ) -> (Option<[u64; $limbs]>, [u64; $limbs])
+            {
+                let mut q = Some([0; $limbs]);
+                let mut r = [0; $limbs];
+
+                for i in (0..($limbs*2*64)).rev() {
+                    // NB: modulo's first bit is unset so this will never
+                    // destroy information
+                    mul2(&mut r);
+                    assert!(set_bit(&mut r, 0, get_bit(this, i)));
+                    if !lt(&r, modulo) {
+                        sub_noborrow(&mut r, modulo);
+                        if q.is_some() && !set_bit(q.as_mut().unwrap(), i, true) {
+                            q = None
+                        }
+                    }
+                }
+
+                if q.is_some() && !lt(q.as_ref().unwrap(), modulo) {
+                    (None, r)
+                } else {
+                    (q, r)
+                }
+            }
+        }
+    }
+}

src/curves/bls381/fq12.rs

@@ -0,0 +1,151 @@
+use super::{Bls381, Fq2, Fq6};
+use rand;
+use super::Field;
+
+#[derive(Debug, Copy, Clone, PartialEq, Eq)]
+pub struct Fq12 {
+    pub c0: Fq6,
+    pub c1: Fq6
+}
+
+impl Fq12 {
+    pub fn unitary_inverse(&mut self, e: &Bls381)
+    {
+        self.c1.negate(e);
+    }
+
+    pub fn mul_by_015(
+        &mut self,
+        e: &Bls381,
+        a: &Fq2,
+        b: &Fq2,
+        c: &Fq2
+    )
+    {
+        let mut aa = self.c0;
+        aa.mul_by_01(e, a, b);
+        let mut bb = self.c1;
+        bb.mul_by_1(e, c);
+        let mut o = *b;
+        o.add_assign(e, &c);
+        self.c1.add_assign(e, &self.c0);
+        self.c1.mul_by_01(e, a, &o);
+        self.c1.sub_assign(e, &aa);
+        self.c1.sub_assign(e, &bb);
+        self.c0 = bb;
+        self.c0.mul_by_nonresidue(e);
+        self.c0.add_assign(e, &aa);
+    }
+}
+
+impl Field<Bls381> for Fq12
+{
+    fn zero() -> Self {
+        Fq12 {
+            c0: Fq6::zero(),
+            c1: Fq6::zero()
+        }
+    }
+
+    fn one(engine: &Bls381) -> Self {
+        Fq12 {
+            c0: Fq6::one(engine),
+            c1: Fq6::zero()
+        }
+    }
+
+    fn random<R: rand::Rng>(engine: &Bls381, rng: &mut R) -> Self {
+        Fq12 {
+            c0: Fq6::random(engine, rng),
+            c1: Fq6::random(engine, rng)
+        }
+    }
+
+    fn is_zero(&self) -> bool {
+        self.c0.is_zero() && self.c1.is_zero()
+    }
+
+    fn double(&mut self, engine: &Bls381) {
+        self.c0.double(engine);
+        self.c1.double(engine);
+    }
+
+    fn negate(&mut self, engine: &Bls381) {
+        self.c0.negate(engine);
+        self.c1.negate(engine);
+    }
+
+    fn add_assign(&mut self, engine: &Bls381, other: &Self) {
+        self.c0.add_assign(engine, &other.c0);
+        self.c1.add_assign(engine, &other.c1);
+    }
+
+    fn sub_assign(&mut self, engine: &Bls381, other: &Self) {
+        self.c0.sub_assign(engine, &other.c0);
+        self.c1.sub_assign(engine, &other.c1);
+    }
+
+    fn frobenius_map(&mut self, e: &Bls381, power: usize)
+    {
+        self.c0.frobenius_map(e, power);
+        self.c1.frobenius_map(e, power);
+
+        self.c1.c0.mul_assign(e, &e.frobenius_coeff_fq12[power % 12]);
+        self.c1.c1.mul_assign(e, &e.frobenius_coeff_fq12[power % 12]);
+        self.c1.c2.mul_assign(e, &e.frobenius_coeff_fq12[power % 12]);
+    }
+
+    fn square(&mut self, e: &Bls381) {
+        let mut ab = self.c0;
+        ab.mul_assign(e, &self.c1);
+        let mut c0c1 = self.c0;
+        c0c1.add_assign(e, &self.c1);
+        let mut c0 = self.c1;
+        c0.mul_by_nonresidue(e);
+        c0.add_assign(e, &self.c0);
+        c0.mul_assign(e, &c0c1);
+        c0.sub_assign(e, &ab);
+        self.c1 = ab;
+        self.c1.add_assign(e, &ab);
+        ab.mul_by_nonresidue(e);
+        c0.sub_assign(e, &ab);
+        self.c0 = c0;
+    }
+
+    fn mul_assign(&mut self, e: &Bls381, other: &Self) {
+        let mut aa = self.c0;
+        aa.mul_assign(e, &other.c0);
+        let mut bb = self.c1;
+        bb.mul_assign(e, &other.c1);
+        let mut o = other.c0;
+        o.add_assign(e, &other.c1);
+        self.c1.add_assign(e, &self.c0);
+        self.c1.mul_assign(e, &o);
+        self.c1.sub_assign(e, &aa);
+        self.c1.sub_assign(e, &bb);
+        self.c0 = bb;
+        self.c0.mul_by_nonresidue(e);
+        self.c0.add_assign(e, &aa);
+    }
+
+    fn inverse(&self, e: &Bls381) -> Option<Self> {
+        let mut c0s = self.c0;
+        c0s.square(e);
+        let mut c1s = self.c1;
+        c1s.square(e);
+        c1s.mul_by_nonresidue(e);
+        c0s.sub_assign(e, &c1s);
+
+        c0s.inverse(e).map(|t| {
+            let mut tmp = Fq12 {
+                c0: t,
+                c1: t
+            };
+            tmp.c0.mul_assign(e, &self.c0);
+            tmp.c1.mul_assign(e, &self.c1);
+            tmp.c1.negate(e);
+
+            tmp
+        })
+    }
+}

src/curves/bls381/fq2.rs

@@ -0,0 +1,156 @@
+use super::{Bls381, Fq};
+use rand;
+use super::{Field, SqrtField};
+
+#[derive(Debug, Copy, Clone, PartialEq, Eq)]
+pub struct Fq2 {
+    pub c0: Fq,
+    pub c1: Fq
+}
+
+impl Fq2 {
+    pub fn mul_by_nonresidue(&mut self, e: &Bls381) {
+        let t0 = self.c0;
+        self.c0.sub_assign(e, &self.c1);
+        self.c1.add_assign(e, &t0);
+    }
+}
+
+impl SqrtField<Bls381> for Fq2 {
+    fn sqrt(&self, engine: &Bls381) -> Option<Self> {
+        // Algorithm 9, https://eprint.iacr.org/2012/685.pdf
+
+        if self.is_zero() {
+            return Some(Self::zero());
+        } else {
+            let mut a1 = self.pow(engine, &engine.fqparams.modulus_minus_3_over_4);
+            let mut alpha = a1;
+            alpha.square(engine);
+            alpha.mul_assign(engine, self);
+            let mut a0 = alpha.pow(engine, &engine.fqparams.modulus);
+            a0.mul_assign(engine, &alpha);
+
+            let mut neg1 = Self::one(engine);
+            neg1.negate(engine);
+
+            if a0 == neg1 {
+                None
+            } else {
+                a1.mul_assign(engine, self);
+
+                if alpha == neg1 {
+                    a1.mul_assign(engine, &Fq2{c0: Fq::zero(), c1: Fq::one(engine)});
+                } else {
+                    alpha.add_assign(engine, &Fq2::one(engine));
+                    alpha = alpha.pow(engine, &engine.fqparams.modulus_minus_1_over_2);
+                    a1.mul_assign(engine, &alpha);
+                }
+
+                Some(a1)
+            }
+        }
+    }
+}
+
+impl Field<Bls381> for Fq2
+{
+    fn zero() -> Self {
+        Fq2 {
+            c0: Fq::zero(),
+            c1: Fq::zero()
+        }
+    }
+
+    fn one(engine: &Bls381) -> Self {
+        Fq2 {
+            c0: Fq::one(engine),
+            c1: Fq::zero()
+        }
+    }
+
+    fn random<R: rand::Rng>(engine: &Bls381, rng: &mut R) -> Self {
+        Fq2 {
+            c0: Fq::random(engine, rng),
+            c1: Fq::random(engine, rng)
+        }
+    }
+
+    fn is_zero(&self) -> bool {
+        self.c0.is_zero() && self.c1.is_zero()
+    }
+
+    fn double(&mut self, engine: &Bls381) {
+        self.c0.double(engine);
+        self.c1.double(engine);
+    }
+
+    fn negate(&mut self, engine: &Bls381) {
+        self.c0.negate(engine);
+        self.c1.negate(engine);
+    }
+
+    fn add_assign(&mut self, engine: &Bls381, other: &Self) {
+        self.c0.add_assign(engine, &other.c0);
+        self.c1.add_assign(engine, &other.c1);
+    }
+
+    fn sub_assign(&mut self, engine: &Bls381, other: &Self) {
+        self.c0.sub_assign(engine, &other.c0);
+        self.c1.sub_assign(engine, &other.c1);
+    }
+
+    fn frobenius_map(&mut self, e: &Bls381, power: usize)
+    {
+        self.c1.mul_assign(e, &e.frobenius_coeff_fq2[power % 2]);
+    }
+
+    fn square(&mut self, engine: &Bls381) {
+        let mut ab = self.c0;
+        ab.mul_assign(engine, &self.c1);
+        let mut c0c1 = self.c0;
+        c0c1.add_assign(engine, &self.c1);
+        let mut c0 = self.c1;
+        c0.negate(engine);
+        c0.add_assign(engine, &self.c0);
+        c0.mul_assign(engine, &c0c1);
+        c0.sub_assign(engine, &ab);
+        self.c1 = ab;
+        self.c1.add_assign(engine, &ab);
+        c0.add_assign(engine, &ab);
+        self.c0 = c0;
+    }
+
+    fn mul_assign(&mut self, engine: &Bls381, other: &Self) {
+        let mut aa = self.c0;
+        aa.mul_assign(engine, &other.c0);
+        let mut bb = self.c1;
+        bb.mul_assign(engine, &other.c1);
+        let mut o = other.c0;
+        o.add_assign(engine, &other.c1);
+        self.c1.add_assign(engine, &self.c0);
+        self.c1.mul_assign(engine, &o);
+        self.c1.sub_assign(engine, &aa);
+        self.c1.sub_assign(engine, &bb);
+        self.c0 = aa;
+        self.c0.sub_assign(engine, &bb);
+    }
+
+    fn inverse(&self, engine: &Bls381) -> Option<Self> {
+        let mut t1 = self.c1;
+        t1.square(engine);
+        let mut t0 = self.c0;
+        t0.square(engine);
+        t0.add_assign(engine, &t1);
+        t0.inverse(engine).map(|t| {
+            let mut tmp = Fq2 {
+                c0: self.c0,
+                c1: self.c1
+            };
+            tmp.c0.mul_assign(engine, &t);
+            tmp.c1.mul_assign(engine, &t);
+            tmp.c1.negate(engine);
+
+            tmp
+        })
+    }
+}

src/curves/bls381/fq6.rs

@@ -0,0 +1,295 @@
+use super::{Bls381, Fq2};
+use rand;
+use super::Field;
+
+#[derive(Debug, Copy, Clone, PartialEq, Eq)]
+pub struct Fq6 {
+    pub c0: Fq2,
+    pub c1: Fq2,
+    pub c2: Fq2
+}
+
+impl Fq6 {
+    pub fn mul_by_nonresidue(&mut self, e: &Bls381) {
+        use std::mem::swap;
+        swap(&mut self.c0, &mut self.c1);
+        swap(&mut self.c0, &mut self.c2);
+
+        self.c0.mul_by_nonresidue(e);
+    }
+
+    pub fn mul_by_1(&mut self, e: &Bls381, c1: &Fq2)
+    {
+        let mut b_b = self.c1;
+        b_b.mul_assign(e, c1);
+
+        let mut t1 = *c1;
+        {
+            let mut tmp = self.c1;
+            tmp.add_assign(e, &self.c2);
+
+            t1.mul_assign(e, &tmp);
+            t1.sub_assign(e, &b_b);
+            t1.mul_by_nonresidue(e);
+        }
+
+        let mut t2 = *c1;
+        {
+            let mut tmp = self.c0;
+            tmp.add_assign(e, &self.c1);
+
+            t2.mul_assign(e, &tmp);
+            t2.sub_assign(e, &b_b);
+        }
+
+        self.c0 = t1;
+        self.c1 = t2;
+        self.c2 = b_b;
+    }
+
+    pub fn mul_by_01(&mut self, e: &Bls381, c0: &Fq2, c1: &Fq2)
+    {
+        let mut a_a = self.c0;
+        let mut b_b = self.c1;
+        a_a.mul_assign(e, c0);
+        b_b.mul_assign(e, c1);
+
+        let mut t1 = *c1;
+        {
+            let mut tmp = self.c1;
+            tmp.add_assign(e, &self.c2);
+
+            t1.mul_assign(e, &tmp);
+            t1.sub_assign(e, &b_b);
+            t1.mul_by_nonresidue(e);
+            t1.add_assign(e, &a_a);
+        }
+
+        let mut t3 = *c0;
+        {
+            let mut tmp = self.c0;
+            tmp.add_assign(e, &self.c2);
+
+            t3.mul_assign(e, &tmp);
+            t3.sub_assign(e, &a_a);
+            t3.add_assign(e, &b_b);
+        }
+
+        let mut t2 = *c0;
+        t2.add_assign(e, c1);
+        {
+            let mut tmp = self.c0;
+            tmp.add_assign(e, &self.c1);
+
+            t2.mul_assign(e, &tmp);
+            t2.sub_assign(e, &a_a);
+            t2.sub_assign(e, &b_b);
+        }
+
+        self.c0 = t1;
+        self.c1 = t2;
+        self.c2 = t3;
+    }
+}
+
+impl Field<Bls381> for Fq6
+{
+    fn zero() -> Self {
+        Fq6 {
+            c0: Fq2::zero(),
+            c1: Fq2::zero(),
+            c2: Fq2::zero()
+        }
+    }
+
+    fn one(engine: &Bls381) -> Self {
+        Fq6 {
+            c0: Fq2::one(engine),
+            c1: Fq2::zero(),
+            c2: Fq2::zero()
+        }
+    }
+
+    fn random<R: rand::Rng>(engine: &Bls381, rng: &mut R) -> Self {
+        Fq6 {
+            c0: Fq2::random(engine, rng),
+            c1: Fq2::random(engine, rng),
+            c2: Fq2::random(engine, rng)
+        }
+    }
+
+    fn is_zero(&self) -> bool {
+        self.c0.is_zero() && self.c1.is_zero() && self.c2.is_zero()
+    }
+
+    fn double(&mut self, engine: &Bls381) {
+        self.c0.double(engine);
+        self.c1.double(engine);
+        self.c2.double(engine);
+    }
+
+    fn negate(&mut self, engine: &Bls381) {
+        self.c0.negate(engine);
+        self.c1.negate(engine);
+        self.c2.negate(engine);
+    }
+
+    fn add_assign(&mut self, engine: &Bls381, other: &Self) {
+        self.c0.add_assign(engine, &other.c0);
+        self.c1.add_assign(engine, &other.c1);
+        self.c2.add_assign(engine, &other.c2);
+    }
+
+    fn sub_assign(&mut self, engine: &Bls381, other: &Self) {
+        self.c0.sub_assign(engine, &other.c0);
+        self.c1.sub_assign(engine, &other.c1);
+        self.c2.sub_assign(engine, &other.c2);
+    }
+
+    fn frobenius_map(&mut self, e: &Bls381, power: usize)
+    {
+        self.c0.frobenius_map(e, power);
+        self.c1.frobenius_map(e, power);
+        self.c2.frobenius_map(e, power);
+
+        self.c1.mul_assign(e, &e.frobenius_coeff_fq6_c1[power % 6]);
+        self.c2.mul_assign(e, &e.frobenius_coeff_fq6_c2[power % 6]);
+    }
+
+    fn square(&mut self, e: &Bls381) {
+        let mut s0 = self.c0;
+        s0.square(e);
+        let mut ab = self.c0;
+        ab.mul_assign(e, &self.c1);
+        let mut s1 = ab;
+        s1.double(e);
+        let mut s2 = self.c0;
+        s2.sub_assign(e, &self.c1);
+        s2.add_assign(e, &self.c2);
+        s2.square(e);
+        let mut bc = self.c1;
+        bc.mul_assign(e, &self.c2);
+        let mut s3 = bc;
+        s3.double(e);
+        let mut s4 = self.c2;
+        s4.square(e);
+
+        self.c0 = s3;
+        self.c0.mul_by_nonresidue(e);
+        self.c0.add_assign(e, &s0);
+
+        self.c1 = s4;
+        self.c1.mul_by_nonresidue(e);
+        self.c1.add_assign(e, &s1);
+
+        self.c2 = s1;
+        self.c2.add_assign(e, &s2);
+        self.c2.add_assign(e, &s3);
+        self.c2.sub_assign(e, &s0);
+        self.c2.sub_assign(e, &s4);
+    }
+
+    fn mul_assign(&mut self, e: &Bls381, other: &Self) {
+        let mut a_a = self.c0;
+        let mut b_b = self.c1;
+        let mut c_c = self.c2;
+        a_a.mul_assign(e, &other.c0);
+        b_b.mul_assign(e, &other.c1);
+        c_c.mul_assign(e, &other.c2);
+
+        let mut t1 = other.c1;
+        t1.add_assign(e, &other.c2);
+        {
+            let mut tmp = self.c1;
+            tmp.add_assign(e, &self.c2);
+
+            t1.mul_assign(e, &tmp);
+            t1.sub_assign(e, &b_b);
+            t1.sub_assign(e, &c_c);
+            t1.mul_by_nonresidue(e);
+            t1.add_assign(e, &a_a);
+        }
+
+        let mut t3 = other.c0;
+        t3.add_assign(e, &other.c2);
+        {
+            let mut tmp = self.c0;
+            tmp.add_assign(e, &self.c2);
+
+            t3.mul_assign(e, &tmp);
+            t3.sub_assign(e, &a_a);
+            t3.add_assign(e, &b_b);
+            t3.sub_assign(e, &c_c);
+        }
+
+        let mut t2 = other.c0;
+        t2.add_assign(e, &other.c1);
+        {
+            let mut tmp = self.c0;
+            tmp.add_assign(e, &self.c1);
+
+            t2.mul_assign(e, &tmp);
+            t2.sub_assign(e, &a_a);
+            t2.sub_assign(e, &b_b);
+            c_c.mul_by_nonresidue(e);
+            t2.add_assign(e, &c_c);
+        }
+
+        self.c0 = t1;
+        self.c1 = t2;
+        self.c2 = t3;
+    }
+
+    fn inverse(&self, e: &Bls381) -> Option<Self> {
+        let mut c0 = self.c2;
+        c0.mul_by_nonresidue(e);
+        c0.mul_assign(e, &self.c1);
+        c0.negate(e);
+        {
+            let mut c0s = self.c0;
+            c0s.square(e);
+            c0.add_assign(e, &c0s);
+        }
+        let mut c1 = self.c2;
+        c1.square(e);
+        c1.mul_by_nonresidue(e);
+        {
+            let mut c01 = self.c0;
+            c01.mul_assign(e, &self.c1);
+            c1.sub_assign(e, &c01);
+        }
+        let mut c2 = self.c1;
+        c2.square(e);
+        {
+            let mut c02 = self.c0;
+            c02.mul_assign(e, &self.c2);
+            c2.sub_assign(e, &c02);
+        }
+
+        let mut tmp1 = self.c2;
+        tmp1.mul_assign(e, &c1);
+        let mut tmp2 = self.c1;
+        tmp2.mul_assign(e, &c2);
+        tmp1.add_assign(e, &tmp2);
+        tmp1.mul_by_nonresidue(e);
+        tmp2 = self.c0;
+        tmp2.mul_assign(e, &c0);
+        tmp1.add_assign(e, &tmp2);
+
+        match tmp1.inverse(e) {
+            Some(t) => {
+                let mut tmp = Fq6 {
+                    c0: t,
+                    c1: t,
+                    c2: t
+                };
+                tmp.c0.mul_assign(e, &c0);
+                tmp.c1.mul_assign(e, &c1);
+                tmp.c2.mul_assign(e, &c2);
+
+                Some(tmp)
+            },
+            None => None
+        }
+    }
+}

src/curves/bls381/tests/mod.rs

@@ -0,0 +1,41 @@
+extern crate bincode;
+
+use curves::*;
+use super::*;
+
+fn test_vectors<E: Engine, G: Group<E>>(e: &E, expected: &[u8]) {
+    let mut bytes = vec![];
+    let mut acc = G::zero(e);
+    let mut expected_reader = expected;
+
+    for _ in 0..10000 {
+        {
+            let acc = acc.to_affine(e);
+            let exp: <G::Affine as GroupAffine<E, G>>::Uncompressed =
+                bincode::deserialize_from(&mut expected_reader, bincode::SizeLimit::Infinite).unwrap();
+
+            assert!(acc == exp.to_affine(e).unwrap());
+
+            let acc = acc.to_uncompressed(e);
+            bincode::serialize_into(&mut bytes, &acc, bincode::SizeLimit::Infinite).unwrap();
+        }
+        acc.double(e);
+        acc.add_assign(e, &G::one(e));
+    }
+
+    assert_eq!(&bytes[..], expected);
+}
+
+#[test]
+fn g1_serialization_test_vectors() {
+    let engine = Bls381::new();
+
+    test_vectors::<Bls381, G1>(&engine, include_bytes!("g1_serialized.bin"));
+}
+
+#[test]
+fn g2_serialization_test_vectors() {
+    let engine = Bls381::new();
+
+    test_vectors::<Bls381, G2>(&engine, include_bytes!("g2_serialized.bin"));
+}

src/curves/mod.rs

@@ -0,0 +1,294 @@
+use rand;
+use std::fmt;
+
+use serde::{Serialize, Deserialize};
+
+pub mod bls381;
+
+pub trait Engine: Sized
+{
+    type Fq: PrimeField<Self>;
+    type Fr: SnarkField<Self>;
+    type Fqe: SqrtField<Self>;
+    type Fqk: Field<Self>;
+    type G1: Group<Self> + Convert<<Self::G1 as Group<Self>>::Affine, Self>;
+    type G2: Group<Self> + Convert<<Self::G2 as Group<Self>>::Affine, Self>;
+
+    fn new() -> Self;
+
+    fn pairing<G1, G2>(&self, p: &G1, q: &G2) -> Self::Fqk
+        where G1: Convert<<Self::G1 as Group<Self>>::Affine, Self>,
+              G2: Convert<<Self::G2 as Group<Self>>::Affine, Self>
+    {
+        self.final_exponentiation(&self.miller_loop(
+            [(
+                &(*p.convert(self)).borrow().prepare(self),
+                &(*q.convert(self)).borrow().prepare(self)
+            )].into_iter()
+        ))
+    }
+    fn miller_loop<'a, I>(&self, I) -> Self::Fqk
+        where I: IntoIterator<Item=&'a (
+                                    &'a <Self::G1 as Group<Self>>::Prepared,
+                                    &'a <Self::G2 as Group<Self>>::Prepared
+                               )>;
+    fn final_exponentiation(&self, &Self::Fqk) -> Self::Fqk;
+
+    fn multiexp<G: Group<Self>>(&self, &[G::Affine], &[Self::Fr]) -> G;
+    fn batch_baseexp<G: Group<Self>>(&self, base: &G, scalars: &[Self::Fr]) -> Vec<G::Affine>;
+}
+
+pub trait Group<E: Engine>: Sized +
+                            Copy +
+                            Clone +
+                            Send +
+                            Sync +
+                            fmt::Debug +
+                            'static
+{
+    type Affine: GroupAffine<E, Self>;
+    type Prepared: Clone + Send + Sync + 'static;
+
+    fn zero(&E) -> Self;
+    fn one(&E) -> Self;
+    fn random<R: rand::Rng>(&E, &mut R) -> Self;
+
+    fn is_zero(&self) -> bool;
+    fn is_equal(&self, &E, other: &Self) -> bool;
+
+    fn to_affine(&self, &E) -> Self::Affine;
+    fn prepare(&self, &E) -> Self::Prepared;
+
+    fn double(&mut self, &E);
+    fn negate(&mut self, engine: &E);
+    fn add_assign(&mut self, &E, other: &Self);
+    fn sub_assign(&mut self, &E, other: &Self);
+    fn add_assign_mixed(&mut self, &E, other: &Self::Affine);
+    fn mul_assign<S: Convert<[u64], E>>(&mut self, &E, other: &S);
+
+    fn optimal_window(&E, scalar_bits: usize) -> Option<usize>;
+}
+
+pub trait GroupAffine<E: Engine, G: Group<E>>: Copy +
+                                               Clone +
+                                               Sized +
+                                               Send +
+                                               Sync +
+                                               fmt::Debug +
+                                               PartialEq +
+                                               Eq +
+                                               'static
+{
+    type Uncompressed: GroupRepresentation<E, G>;
+
+    fn to_jacobian(&self, &E) -> G;
+    fn prepare(self, &E) -> G::Prepared;
+    fn is_zero(&self) -> bool;
+    fn mul<S: Convert<[u64], E>>(&self, &E, other: &S) -> G;
+    fn negate(&mut self, &E);
+
+    /// Returns true iff the point is on the curve and in the correct
+    /// subgroup. This is guaranteed to return true unless the user
+    /// invokes `to_affine_unchecked`.
+    fn is_valid(&self, &E) -> bool;
+
+    /// Produces an "uncompressed" representation of the curve point according
+    /// to IEEE standards.
+    fn to_uncompressed(&self, &E) -> Self::Uncompressed;
+}
+
+pub trait GroupRepresentation<E: Engine, G: Group<E>>: Serialize + Deserialize
+{
+    /// If the point representation is valid (lies on the curve, correct
+    /// subgroup) this function will return it.
+    fn to_affine(&self, e: &E) -> Result<G::Affine, ()> {
+        let p = try!(self.to_affine_unchecked(e));
+
+        if p.is_valid(e) {
+            Ok(p)
+        } else {
+            Err(())
+        }
+    }
+
+    /// Returns the point under the assumption that it is valid. Undefined
+    /// behavior if `to_affine` would have rejected the point.
+    fn to_affine_unchecked(&self, &E) -> Result<G::Affine, ()>;
+}
+
+pub trait Field<E: Engine>: Sized +
+                            Eq +
+                            PartialEq +
+                            Copy +
+                            Clone +
+                            Send +
+                            Sync +
+                            fmt::Debug +
+                            'static
+{
+    fn zero() -> Self;
+    fn one(&E) -> Self;
+    fn random<R: rand::Rng>(&E, &mut R) -> Self;
+    fn is_zero(&self) -> bool;
+    fn square(&mut self, engine: &E);
+    fn double(&mut self, engine: &E);
+    fn negate(&mut self, &E);
+    fn add_assign(&mut self, &E, other: &Self);
+    fn sub_assign(&mut self, &E, other: &Self);
+    fn mul_assign(&mut self, &E, other: &Self);
+    fn inverse(&self, &E) -> Option<Self>;
+    fn frobenius_map(&mut self, &E, power: usize);
+    fn pow<S: Convert<[u64], E>>(&self, engine: &E, exp: &S) -> Self
+    {
+        let mut res = Self::one(engine);
+
+        for i in BitIterator::from((*exp.convert(engine)).borrow()) {
+            res.square(engine);
+            if i {
+                res.mul_assign(engine, self);
+            }
+        }
+
+        res
+    }
+}
+
+pub trait SqrtField<E: Engine>: Field<E>
+{
+    /// Returns the square root of the field element, if it is
+    /// quadratic residue.
+    fn sqrt(&self, engine: &E) -> Option<Self>;
+}
+
+pub trait PrimeField<E: Engine>: SqrtField<E> + Convert<[u64], E>
+{
+    /// Little endian representation of a field element.
+    type Repr: Convert<[u64], E>;
+    fn from_u64(&E, u64) -> Self;
+    fn from_str(&E, s: &str) -> Result<Self, ()>;
+    fn from_repr(&E, Self::Repr) -> Result<Self, ()>;
+    fn into_repr(&self, &E) -> Self::Repr;
+
+    /// Returns an interator over all bits, most significant bit first.
+    fn bits(&self, &E) -> BitIterator<Self::Repr>;
+
+    /// Returns the field characteristic; the modulus.
+    fn char(&E) -> Self::Repr;
+
+    /// Returns how many bits are needed to represent an element of this
+    /// field.
+    fn num_bits(&E) -> usize;
+
+    /// Returns how many bits of information can be reliably stored in the
+    /// field element.
+    fn capacity(&E) -> usize;
+}
+
+pub trait SnarkField<E: Engine>: PrimeField<E>
+{
+    fn s(&E) -> u64;
+    fn multiplicative_generator(&E) -> Self;
+    fn root_of_unity(&E) -> Self;
+}
+
+pub struct BitIterator<T> {
+    t: T,
+    n: usize
+}
+
+impl<T: AsRef<[u64]>> Iterator for BitIterator<T> {
+    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)
+        }
+    }
+}
+
+impl<'a> From<&'a [u64]> for BitIterator<&'a [u64]>
+{
+    fn from(v: &'a [u64]) -> Self {
+        assert!(v.len() < 100);
+
+        BitIterator {
+            t: v,
+            n: v.len() * 64
+        }
+    }
+}
+
+use std::ops::Deref;
+use std::borrow::Borrow;
+
+pub enum Cow<'a, T: 'a> {
+    Owned(T),
+    Borrowed(&'a T)
+}
+
+impl<'a, T: 'a> Deref for Cow<'a, T> {
+    type Target = T;
+
+    fn deref(&self) -> &T {
+        match *self {
+            Cow::Owned(ref v) => v,
+            Cow::Borrowed(v) => v
+        }
+    }
+}
+
+pub trait Convert<T: ?Sized, E> {
+    type Target: Borrow<T>;
+
+    fn convert(&self, &E) -> Cow<Self::Target>;
+}
+
+impl<T, E> Convert<T, E> for T {
+    type Target = T;
+
+    fn convert(&self, _: &E) -> Cow<T> {
+        Cow::Borrowed(self)
+    }
+}
+
+macro_rules! bit_iter_impl(
+    ($n:expr) => {
+        impl From<[u64; $n]> for BitIterator<[u64; $n]> {
+            fn from(v: [u64; $n]) -> Self {
+                BitIterator {
+                    t: v,
+                    n: $n * 64
+                }
+            }
+        }
+
+        impl<E> Convert<[u64], E> for [u64; $n] {
+            type Target = [u64; $n];
+
+            fn convert(&self, _: &E) -> Cow<[u64; $n]> {
+                Cow::Borrowed(self)
+            }
+        }
+    };
+);
+
+bit_iter_impl!(1);
+bit_iter_impl!(2);
+bit_iter_impl!(3);
+bit_iter_impl!(4);
+bit_iter_impl!(5);
+bit_iter_impl!(6);
+
+#[cfg(test)]
+mod tests;
+
+#[test]
+fn bls381_test_suite() {
+    tests::test_engine::<bls381::Bls381>();
+}

src/curves/tests/fields.rs

@@ -0,0 +1,219 @@
+use rand::{self, Rng};
+use super::super::{Engine, Field, SqrtField, PrimeField};
+
+fn inversion_tests<E: Engine, F: Field<E>, R: Rng>(e: &E, rng: &mut R) {
+    let mut a = F::one(e);
+    for _ in 0..10000 {
+        let mut b = a.inverse(e).unwrap();
+        b.mul_assign(e, &a);
+        assert_eq!(b, F::one(e));
+        a.add_assign(e, &F::one(e));
+    }
+    a = F::one(e);
+    a.negate(e);
+    for _ in 0..10000 {
+        let mut b = a.inverse(e).unwrap();
+        b.mul_assign(e, &a);
+        assert_eq!(b, F::one(e));
+        a.sub_assign(e, &F::one(e));
+    }
+    a = F::zero();
+    assert!(a.inverse(e).is_none());
+    for _ in 0..10000 {
+        let r = F::random(e, rng);
+        assert!(!r.is_zero());
+        let mut rinv = r.inverse(e).unwrap();
+        rinv.mul_assign(e, &r);
+        assert_eq!(rinv, F::one(e));
+    }
+}
+
+fn expansion_tests<E: Engine, F: Field<E>, R: Rng>(e: &E, rng: &mut R) {
+    for _ in 0..100 {
+        let a = F::random(e, rng);
+        let b = F::random(e, rng);
+        let c = F::random(e, rng);
+        let d = F::random(e, rng);
+
+        let lhs;
+        {
+            let mut t0 = a;
+            t0.add_assign(e, &b);
+            let mut t1 = c;
+            t1.add_assign(e, &d);
+            t0.mul_assign(e, &t1);
+            lhs = t0;
+        }
+
+        let rhs;
+        {
+            let mut t0 = a;
+            t0.mul_assign(e, &c);
+            let mut t1 = b;
+            t1.mul_assign(e, &c);
+            let mut t2 = a;
+            t2.mul_assign(e, &d);
+            let mut t3 = b;
+            t3.mul_assign(e, &d);
+            t0.add_assign(e, &t1);
+            t0.add_assign(e, &t2);
+            t0.add_assign(e, &t3);
+            rhs = t0;
+        }
+
+        assert_eq!(lhs, rhs);
+    }
+}
+
+fn squaring_tests<E: Engine, F: Field<E>, R: Rng>(e: &E, rng: &mut R) {
+    for _ in 0..100 {
+        let mut a = F::random(e, rng);
+        let mut b = a;
+        b.mul_assign(e, &a);
+        a.square(e);
+
+        assert_eq!(a, b);
+    }
+
+    let mut cur = F::zero();
+    for _ in 0..100 {
+        let mut a = cur;
+        a.square(e);
+        let mut b = cur;
+        b.mul_assign(e, &cur);
+
+        assert_eq!(a, b);
+
+        cur.add_assign(e, &F::one(e));
+    }
+}
+
+fn operation_tests<E: Engine, F: Field<E>, R: Rng>(e: &E, rng: &mut R) {
+    {
+        let mut acc = F::zero();
+        for _ in 0..1000 {
+            let mut a = acc;
+            a.negate(e);
+            a.add_assign(e, &acc);
+
+            assert_eq!(a, F::zero());
+            acc.add_assign(e, &F::one(e));
+        }
+    }
+    {
+        for _ in 0..1000 {
+            let mut a = F::random(e, rng);
+            let mut at = a;
+            let mut b = F::random(e, rng);
+
+            a.sub_assign(e, &b);
+            b.negate(e);
+            at.add_assign(e, &b);
+
+            assert_eq!(a, at);
+        }
+    }
+}
+
+pub fn test_field<E: Engine, F: Field<E>>(e: &E) {
+    let rng = &mut rand::thread_rng();
+
+    inversion_tests::<E, F, _>(e, rng);
+    expansion_tests::<E, F, _>(e, rng);
+    squaring_tests::<E, F, _>(e, rng);
+    operation_tests::<E, F, _>(e, rng);
+}
+
+pub fn test_sqrt_field<E: Engine, F: SqrtField<E>>(e: &E) {
+    const SAMPLES: isize = 10000;
+
+    {
+        let mut acc = F::one(e);
+
+        for _ in 0..SAMPLES {
+            let mut b = acc;
+            b.square(e);
+            let mut c = b.sqrt(e).unwrap();
+            if c != acc {
+                c.negate(e);
+            }
+
+            assert_eq!(acc, c);
+
+            acc.add_assign(e, &F::one(e));
+        }
+    }
+
+    {
+        let mut acc = F::one(e);
+
+        for _ in 0..SAMPLES {
+            match acc.sqrt(e) {
+                Some(mut a) => {
+                    a.square(e);
+
+                    assert_eq!(a, acc);
+                },
+                None => {}
+            }
+
+            acc.add_assign(e, &F::one(e));
+        }
+    }
+
+    {
+        let rng = &mut rand::thread_rng();
+
+        for _ in 0..SAMPLES {
+            let a = F::random(e, rng);
+            let mut b = a;
+            b.square(e);
+            let mut c = b.sqrt(e).unwrap();
+            if c != a {
+                c.negate(e);
+            }
+
+            assert_eq!(a, c);
+        }
+    }
+
+    {
+        let rng = &mut rand::thread_rng();
+
+        let mut qr: isize = 0;
+        let mut nqr: isize = 0;
+        for _ in 0..SAMPLES {
+            let a = F::random(e, rng);
+            match a.sqrt(e) {
+                Some(mut b) => {
+                    qr += 1;
+                    b.square(e);
+                    assert_eq!(a, b);
+                },
+                None => {
+                    nqr += 1;
+                }
+            }
+        }
+
+        assert!((qr - nqr < (SAMPLES / 20)) || (qr - nqr > -(SAMPLES / 20)));
+    }
+}
+
+pub fn test_prime_field<E: Engine, F: PrimeField<E>>(e: &E) {
+    let rng = &mut rand::thread_rng();
+
+    for _ in 0..100 {
+        let a = F::random(e, rng);
+        let b = F::random(e, rng);
+        let mut c = a;
+        c.mul_assign(e, &b);
+        let a = a.into_repr(e);
+        let b = b.into_repr(e);
+        let expected_a = F::from_repr(e, a).unwrap();
+        let expected_b = F::from_repr(e, b).unwrap();
+        let mut expected_c = expected_a;
+        expected_c.mul_assign(e, &expected_b);
+        assert_eq!(c, expected_c);
+    }
+}

src/curves/tests/groups.rs

@@ -0,0 +1,277 @@
+use rand;
+use super::super::{Engine, Field, PrimeField, Group, GroupAffine};
+
+fn random_test_mixed_addition<E: Engine, G: Group<E>>(e: &E)
+{
+    let rng = &mut rand::thread_rng();
+
+    // affine is zero
+    {
+        let a = G::zero(e).to_affine(e);
+        let mut b = G::random(e, rng);
+        let bcpy = b;
+
+        b.add_assign_mixed(e, &a);
+
+        assert!(bcpy.is_equal(e, &b));
+        assert_eq!(bcpy.to_affine(e), b.to_affine(e));
+    }
+
+    // self is zero
+    {
+        let a = G::random(e, rng).to_affine(e);
+        let mut b = G::zero(e);
+        let acpy = a.to_jacobian(e);
+
+        b.add_assign_mixed(e, &a);
+
+        assert!(acpy.is_equal(e, &b));
+        assert_eq!(acpy.to_affine(e), b.to_affine(e));
+    }
+
+    // both are zero
+    {
+        let a = G::zero(e).to_affine(e);
+        let mut b = G::zero(e);
+        let acpy = a.to_jacobian(e);
+
+        b.add_assign_mixed(e, &a);
+
+        assert!(acpy.is_equal(e, &b));
+        assert_eq!(acpy.to_affine(e), b.to_affine(e));
+    }
+
+    // one is negative of the other
+    {
+        let a = G::random(e, rng);
+        let mut b = a;
+        b.negate(e);
+        let a = a.to_affine(e);
+
+        b.add_assign_mixed(e, &a);
+        assert!(b.is_zero());
+        assert_eq!(b.to_affine(e), G::zero(e).to_affine(e));
+    }
+
+    // doubling case
+    {
+        let a = G::random(e, rng);
+        let b = a.to_affine(e);
+        let mut acpy = a;
+        acpy.add_assign_mixed(e, &b);
+
+        let mut t = a;
+        t.double(e);
+        assert!(acpy.is_equal(e, &t));
+    }
+
+    for _ in 0..100 {
+        let mut x = G::random(e, rng);
+        let mut y = x;
+        let b = G::random(e, rng);
+        let baffine = b.to_affine(e);
+
+        x.add_assign(e, &b);
+        y.add_assign_mixed(e, &baffine);
+
+        assert!(x.is_equal(e, &y));
+    }
+}
+
+fn random_test_addition<E: Engine, G: Group<E>>(e: &E) {
+    let rng = &mut rand::thread_rng();
+
+    for _ in 0..50 {
+        let r1 = G::random(e, rng);
+        let r2 = G::random(e, rng);
+        let r3 = G::random(e, rng);
+
+        {
+            let mut tmp1 = r1;
+            tmp1.add_assign(e, &r2);
+            tmp1.add_assign(e, &r3);
+
+            let mut tmp2 = r2;
+            tmp2.add_assign(e, &r3);
+            tmp2.add_assign(e, &r1);
+
+            assert!(tmp1.is_equal(e, &tmp2));
+        }
+
+        {
+            let mut tmp = r1;
+            tmp.add_assign(e, &r2);
+            tmp.add_assign(e, &r3);
+            tmp.sub_assign(e, &r1);
+            tmp.sub_assign(e, &r2);
+            tmp.sub_assign(e, &r3);
+
+            assert!(tmp.is_zero());
+        }
+    }
+}
+
+fn random_test_doubling<E: Engine, G: Group<E>>(e: &E) {
+    let rng = &mut rand::thread_rng();
+
+    for _ in 0..50 {
+        let r1 = G::random(e, rng);
+        let r2 = G::random(e, rng);
+        let ti = E::Fr::from_str(e, "2").unwrap().inverse(e).unwrap();
+
+        {
+            let mut tmp_1 = r1;
+            tmp_1.add_assign(e, &r2);
+            tmp_1.add_assign(e, &r1);
+
+            let mut tmp_2 = r1;
+            tmp_2.double(e);
+            tmp_2.add_assign(e, &r2);
+
+            assert!(tmp_1.is_equal(e, &tmp_2));
+        }
+
+        {
+            let mut tmp = r1;
+            tmp.double(e);
+            tmp.mul_assign(e, &ti);
+
+            assert!(tmp.is_equal(e, &r1));
+        }
+    }
+}
+
+fn random_test_dh<E: Engine, G: Group<E>>(e: &E) {
+    let rng = &mut rand::thread_rng();
+
+    for _ in 0..50 {
+        let alice_sk = E::Fr::random(e, rng);
+        let bob_sk = E::Fr::random(e, rng);
+
+        let mut alice_pk = G::one(e);
+        alice_pk.mul_assign(e, &alice_sk);
+        let mut bob_pk = G::one(e);
+        bob_pk.mul_assign(e, &bob_sk);
+
+        let mut alice_shared = bob_pk;
+        alice_shared.mul_assign(e, &alice_sk);
+        let mut bob_shared = alice_pk;
+        bob_shared.mul_assign(e, &bob_sk);
+
+        assert!(alice_shared.is_equal(e, &bob_shared));
+    }
+}
+
+fn random_mixed_addition<E: Engine, G: Group<E>>(e: &E) {
+    let rng = &mut rand::thread_rng();
+
+    for _ in 0..50 {
+        let a = G::random(e, rng);
+
+        let mut res = a;
+        res.double(e);
+
+        let affine = a.to_affine(e);
+        let mut jacobian = affine.to_jacobian(e);
+        jacobian.double(e);
+
+        assert!(jacobian.is_equal(e, &res));
+    }
+}
+
+fn random_test_equality<E: Engine, G: Group<E>>(e: &E) {
+    let rng = &mut rand::thread_rng();
+
+    for _ in 0..50 {
+        let begin = G::random(e, rng);
+
+        let mut acc = begin;
+
+        let a = E::Fr::random(e, rng);
+        let b = G::random(e, rng);
+        let c = E::Fr::random(e, rng);
+        let d = G::random(e, rng);
+
+        for _ in 0..10 {
+            acc.mul_assign(e, &a);
+            acc.negate(e);
+            acc.add_assign(e, &b);
+            acc.mul_assign(e, &c);
+            acc.negate(e);
+            acc.sub_assign(e, &d);
+            acc.double(e);
+        }
+
+        assert!(!acc.is_equal(e, &begin));
+
+        let ai = a.inverse(e).unwrap();
+        let ci = c.inverse(e).unwrap();
+        let ti = E::Fr::from_str(e, "2").unwrap().inverse(e).unwrap();
+
+        for _ in 0..10 {
+            acc.mul_assign(e, &ti);
+            acc.add_assign(e, &d);
+            acc.negate(e);
+            acc.mul_assign(e, &ci);
+            acc.sub_assign(e, &b);
+            acc.negate(e);
+            acc.mul_assign(e, &ai);
+        }
+
+        assert!(acc.is_equal(e, &begin));
+    }
+}
+
+pub fn test_group<E: Engine, G: Group<E>>(e: &E) {
+    {
+        let rng = &mut rand::thread_rng();
+        let mut g = G::random(e, rng);
+        let order = <E::Fr as PrimeField<E>>::char(e);
+        g.mul_assign(e, &order);
+
+        assert!(g.is_zero());
+    }
+    {
+        let rng = &mut rand::thread_rng();
+        let mut neg1 = E::Fr::one(e);
+        neg1.negate(e);
+        for _ in 0..1000 {
+            let orig = G::random(e, rng);
+            let mut a = orig;
+            a.mul_assign(e, &neg1);
+            assert!(!a.is_zero());
+            a.add_assign(e, &orig);
+            assert!(a.is_zero());
+        }
+    }
+    {
+        let mut o = G::one(e);
+        o.sub_assign(e, &G::one(e));
+        assert!(o.is_zero());
+    }
+    {
+        let mut o = G::one(e);
+        o.add_assign(e, &G::one(e));
+        let mut r = G::one(e);
+        r.mul_assign(e, &E::Fr::from_str(e, "2").unwrap());
+        assert!(o.is_equal(e, &r));
+    }
+    {
+        let mut z = G::zero(e);
+        assert!(z.is_zero());
+        z.double(e);
+        assert!(z.is_zero());
+
+        let zaffine = z.to_affine(e);
+        let zjacobian = zaffine.to_jacobian(e);
+
+        assert!(zjacobian.is_zero());
+    }
+
+    random_test_equality::<E, G>(e);
+    random_test_dh::<E, G>(e);
+    random_test_doubling::<E, G>(e);
+    random_test_addition::<E, G>(e);
+    random_mixed_addition::<E, G>(e);
+    random_test_mixed_addition::<E, G>(e);
+}

src/curves/tests/mod.rs

@@ -0,0 +1,136 @@
+use super::{Engine, Group, GroupAffine, Field, PrimeField};
+use rand;
+
+mod fields;
+mod groups;
+
+fn test_multiexp<E: Engine, G: Group<E>>(e: &E) {
+    fn naiveexp<E: Engine, G: Group<E>>(e: &E, g: &[G::Affine], s: &[E::Fr]) -> G
+    {
+        assert!(g.len() == s.len());
+
+        let mut expected = G::zero(e);
+        for (g, s) in g.iter().zip(s.iter()) {
+            expected.add_assign(e, &g.mul(e, s));
+        }
+
+        expected
+    }
+
+    {
+        let rng = &mut rand::thread_rng();
+
+        let g: Vec<G::Affine> = (0..1000).map(|_| G::random(e, rng).to_affine(e)).collect();
+        let s: Vec<E::Fr> = (0..1000).map(|_| E::Fr::random(e, rng)).collect();
+
+        let naive = naiveexp::<E, G>(e, &g, &s);
+        let multi = e.multiexp::<G>(&g, &s);
+
+        assert!(naive.is_equal(e, &multi));
+        assert!(multi.is_equal(e, &naive));
+    }
+
+    {
+        let rng = &mut rand::thread_rng();
+        let g: Vec<G::Affine> = (0..2).map(|_| G::random(e, rng).to_affine(e)).collect();
+        let s = vec![E::Fr::from_str(e, "3435973836800000000000000000000000").unwrap(), E::Fr::from_str(e, "3435973836700000000000000000000000").unwrap()];
+
+        let naive = naiveexp::<E, G>(e, &g, &s);
+        let multi = e.multiexp::<G>(&g, &s);
+
+        assert!(naive.is_equal(e, &multi));
+        assert!(multi.is_equal(e, &naive));
+    }
+}
+
+fn test_bilinearity<E: Engine>(e: &E) {
+    let rng = &mut rand::thread_rng();
+
+    let a = E::G1::random(e, rng);
+    let b = E::G2::random(e, rng);
+    let s = E::Fr::random(e, rng);
+
+    let mut a_s = a;
+    a_s.mul_assign(e, &s);
+
+    let mut b_s = b;
+    b_s.mul_assign(e, &s);
+
+    let test1 = e.pairing(&a_s, &b);
+    assert!(test1 != E::Fqk::one(e));
+    let test2 = e.pairing(&a, &b_s);
+    assert_eq!(test1, test2);
+    
+    let mut test4 = e.pairing(&a, &b);
+    assert!(test4 != test1);
+    test4 = test4.pow(e, &s);
+    assert_eq!(test1, test4);
+}
+
+fn test_multimiller<E: Engine>(e: &E) {
+    let rng = &mut rand::thread_rng();
+
+    let a1 = E::G1::random(e, rng);
+    let a2 = E::G2::random(e, rng);
+
+    let b1 = E::G1::random(e, rng);
+    let b2 = E::G2::random(e, rng);
+
+    let mut p1 = e.pairing(&a1, &a2);
+    let p2 = e.pairing(&b1, &b2);
+    p1.mul_assign(e, &p2);
+
+    let mm = e.final_exponentiation(&e.miller_loop(
+        [
+            (&a1.prepare(e), &a2.prepare(e)),
+            (&b1.prepare(e), &b2.prepare(e))
+        ].into_iter()
+    ));
+
+    assert_eq!(p1, mm);
+}
+
+pub fn test_engine<E: Engine>() {
+    let engine = E::new();
+
+    fields::test_prime_field::<E, E::Fq>(&engine);
+    fields::test_prime_field::<E, E::Fr>(&engine);
+    fields::test_sqrt_field::<E, E::Fq>(&engine);
+    fields::test_sqrt_field::<E, E::Fr>(&engine);
+    fields::test_sqrt_field::<E, E::Fqe>(&engine);
+
+    fields::test_field::<E, E::Fq>(&engine);
+    fields::test_field::<E, E::Fr>(&engine);
+    fields::test_field::<E, E::Fqe>(&engine);
+    fields::test_field::<E, E::Fqk>(&engine);
+
+    groups::test_group::<E, E::G1>(&engine);
+    groups::test_group::<E, E::G2>(&engine);
+
+    test_bilinearity(&engine);
+    test_multimiller(&engine);
+    test_frobenius(&engine);
+    test_multiexp::<E, E::G1>(&engine);
+    test_multiexp::<E, E::G2>(&engine);
+}
+
+fn test_frobenius<E: Engine>(e: &E) {
+    let rng = &mut rand::thread_rng();
+    let modulus = E::Fq::char(e);
+
+    let a = E::Fqk::random(e, rng);
+    let mut acpy = a;
+    acpy.frobenius_map(e, 0);
+    assert_eq!(acpy, a);
+
+    let mut a_q = a.pow(e, &modulus);
+
+    for p in 1..12 {
+        acpy = a;
+        acpy.frobenius_map(e, p);
+
+        assert_eq!(acpy, a_q);
+
+        a_q = a_q.pow(e, &modulus);
+    }
+}

src/lib.rs

@@ -0,0 +1,8 @@
+#![feature(i128_type)]
+
+extern crate rand;
+extern crate rayon;
+extern crate byteorder;
+extern crate serde;
+
+pub mod curves;