New arithmetic and G1/G2 implementations.

benches/fp.rs

@@ -7,6 +7,42 @@ extern crate bn;
 
 use bn::*;
 
+#[bench]
+fn fr_addition(b: &mut test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| Fr::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| Fr::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] + v2[ctr % SAMPLES]
+    });
+}
+
+#[bench]
+fn fr_subtraction(b: &mut test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| Fr::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| Fr::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] - v2[ctr % SAMPLES]
+    });
+}
+
 #[bench]
 fn fr_multiplication(b: &mut test::Bencher) {
     const SAMPLES: usize = 1000;
@@ -42,6 +78,42 @@ fn fr_inverses(b: &mut test::Bencher) {
     });
 }
 
+#[bench]
+fn fq_addition(b: &mut test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| Fq::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| Fq::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] + v2[ctr % SAMPLES]
+    });
+}
+
+#[bench]
+fn fq_subtraction(b: &mut test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| Fq::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| Fq::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] - v2[ctr % SAMPLES]
+    });
+}
+
 #[bench]
 fn fq_multiplication(b: &mut test::Bencher) {
     const SAMPLES: usize = 1000;
@@ -95,6 +167,42 @@ fn fq2_multiplication(b: &mut test::Bencher) {
     });
 }
 
+#[bench]
+fn fq2_addition(b: &mut test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| Fq2::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| Fq2::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] + v2[ctr % SAMPLES]
+    });
+}
+
+#[bench]
+fn fq2_subtraction(b: &mut test::Bencher) {
+    const SAMPLES: usize = 1000;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| Fq2::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| Fq2::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] - v2[ctr % SAMPLES]
+    });
+}
+
 #[bench]
 fn fq2_inverses(b: &mut test::Bencher) {
     const SAMPLES: usize = 1000;
@@ -111,3 +219,111 @@ fn fq2_inverses(b: &mut test::Bencher) {
         v1[ctr % SAMPLES].inverse()
     });
 }
+
+#[bench]
+fn g1_addition(b: &mut test::Bencher) {
+    const SAMPLES: usize = 100;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| G1::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| G1::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] + v2[ctr % SAMPLES]
+    });
+}
+
+#[bench]
+fn g1_subtraction(b: &mut test::Bencher) {
+    const SAMPLES: usize = 100;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| G1::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| G1::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] - v2[ctr % SAMPLES]
+    });
+}
+
+#[bench]
+fn g1_scalar_multiplication(b: &mut test::Bencher) {
+    const SAMPLES: usize = 100;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| G1::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| Fr::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] * v2[ctr % SAMPLES]
+    });
+}
+
+#[bench]
+fn g2_addition(b: &mut test::Bencher) {
+    const SAMPLES: usize = 100;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| G2::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| G2::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] + v2[ctr % SAMPLES]
+    });
+}
+
+#[bench]
+fn g2_subtraction(b: &mut test::Bencher) {
+    const SAMPLES: usize = 100;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| G2::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| G2::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] - v2[ctr % SAMPLES]
+    });
+}
+
+#[bench]
+fn g2_scalar_multiplication(b: &mut test::Bencher) {
+    const SAMPLES: usize = 100;
+
+    let rng = &mut rand::thread_rng();
+
+    let v1: Vec<_> = (0..SAMPLES).map(|_| G2::random(rng)).collect();
+    let v2: Vec<_> = (0..SAMPLES).map(|_| Fr::random(rng)).collect();
+
+    let mut ctr = 0;
+
+    b.iter(|| {
+        ctr += 1;
+
+        v1[ctr % SAMPLES] * v2[ctr % SAMPLES]
+    });
+}

src/arith.rs

@@ -1,13 +1,10 @@
 use std::cmp::Ordering;
 use rand::Rng;
 
-mod primitives;
-use self::primitives::*;
-
 /// 256-bit, stack allocated biginteger for use in prime field
 /// arithmetic.
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]
-pub struct U256([Digit; LIMBS]);
+pub struct U256([u64; 4]);
 
 impl Ord for U256 {
     #[inline]
@@ -25,21 +22,15 @@ impl Ord for U256 {
 }
 
 impl PartialOrd for U256 {
+    #[inline]
     fn partial_cmp(&self, other: &U256) -> Option<Ordering> {
         Some(self.cmp(other))
     }
 }
 
-impl From<u32> for U256 {
+impl From<[u64; 4]> for U256 {
     #[inline]
-    fn from(f: u32) -> U256 {
+    fn from(a: [u64; 4]) -> U256 {
-        U256([f, 0, 0, 0, 0, 0, 0, 0])
-    }
-}
-
-impl From<[u32; 8]> for U256 {
-    #[inline]
-    fn from(a: [u32; 8]) -> U256 {
         U256(a)
     }
 }
@@ -47,19 +38,16 @@ impl From<[u32; 8]> for U256 {
 impl U256 {
     #[inline]
     pub fn zero() -> U256 {
-        0.into()
+        [0, 0, 0, 0].into()
     }
 
     #[inline]
     pub fn one() -> U256 {
-        1.into()
+        [1, 0, 0, 0].into()
     }
 
     /// Produce a random number (mod `modulo`)
-    pub fn rand<R: Rng>(
+    pub fn random<R: Rng>(rng: &mut R, modulo: &U256) -> U256
-        rng: &mut R,
-        modulo: &U256
-    ) -> U256
     {
         let mut res;
 
@@ -94,8 +82,8 @@ impl U256 {
     {
         assert!(n < 256);
 
-        let part = n / 32;
+        let part = n / 64;
-        let bit = n - (32 * part);
+        let bit = n - (64 * part);
 
         if to {
             self.0[part] |= 1 << bit;
@@ -124,7 +112,7 @@ impl U256 {
 
     /// Multiply `self` by `other` (mod `modulo`) via the Montgomery
     /// multiplication method.
-    pub fn mul(&mut self, other: &U256, modulo: &U256, inv: u32) {
+    pub fn mul(&mut self, other: &U256, modulo: &U256, inv: u64) {
         mul_reduce(&mut self.0, &other.0, &modulo.0, inv);
 
         if *self >= *modulo {
@@ -221,20 +209,156 @@ impl<'a> Iterator for BitIterator<'a> {
         else {
             self.n -= 1;
 
-            let part = self.n / 32;
+            let part = self.n / 64;
-            let bit = self.n - (32 * part);
+            let bit = self.n - (64 * part);
 
             Some(self.int.0[part] & (1 << bit) > 0)
         }
     }
 }
 
+/// Divide by two
+#[inline]
+fn div2(a: &mut [u64; 4]) {
+    let mut b = 0;
+    for a in a.iter_mut().rev() {
+        let t = *a << 63;
+        *a = *a >> 1;
+        *a = *a | b;
+        b = t;
+    }
+}
+
+#[inline]
+fn split_u64(i: u64) -> (u64, u64) {
+    (i >> 32, i & 0xFFFFFFFF)
+}
+
+#[inline]
+fn combine_u64(hi: u64, lo: u64) -> u64 {
+    (hi << 32) | lo
+}
+
+#[inline]
+fn adc(a: u64, b: u64, carry: &mut u64) -> u64 {
+    let (a1, a0) = split_u64(a);
+    let (b1, b0) = split_u64(b);
+    let (c, r0) = split_u64(a0 + b0 + *carry);
+    let (c, r1) = split_u64(a1 + b1 + c);
+    *carry = c;
+
+    combine_u64(r1, r0)
+}
+
+#[inline]
+fn add_nocarry(a: &mut [u64; 4], b: &[u64; 4]) {
+    let mut carry = 0;
+
+    for (a, b) in a.into_iter().zip(b.iter()) {
+        *a = adc(*a, *b, &mut carry);
+    }
+
+    debug_assert!(0 == carry);
+}
+
+#[inline]
+fn sub_noborrow(a: &mut [u64; 4], b: &[u64; 4]) {
+    #[inline]
+    fn sbb(a: u64, b: u64, borrow: &mut u64) -> u64 {
+        let (a1, a0) = split_u64(a);
+        let (b1, b0) = split_u64(b);
+        let (b, r0) = split_u64((1 << 32) + a0 - b0 - *borrow);
+        let (b, r1) = split_u64((1 << 32) + a1 - b1 - ((b == 0) as u64));
+
+        *borrow = (b == 0) as u64;
+
+        combine_u64(r1, r0)
+    }
+
+    let mut borrow = 0;
+
+    for (a, b) in a.into_iter().zip(b.iter()) {
+        *a = sbb(*a, *b, &mut borrow);
+    }
+
+    debug_assert!(0 == borrow);
+}
+
+#[inline]
+fn mul_reduce(
+    this: &mut [u64; 4],
+    by: &[u64; 4],
+    modulus: &[u64; 4],
+    inv: u64
+)
+{
+    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 (b_hi, b_lo) = split_u64(b);
+            let (c_hi, c_lo) = split_u64(c);
+
+            let (x_hi, x_lo) = split_u64(b_lo * c_lo);
+            let (y_hi, y_lo) = split_u64(b_lo * c_hi);
+            let (z_hi, z_lo) = split_u64(b_hi * c_lo);
+            let (r_hi, r_lo) = split_u64(x_hi + y_lo + z_lo);
+            let mut c0 = (r_lo << 32) | x_lo;
+            let mut c1 = (b_hi * c_hi) + r_hi + y_hi + z_hi;
+
+            let mut c = 0;
+            c0 = adc(c0, *carry, &mut c);
+            c1 = adc(c1, 0, &mut c);
+            c0 = adc(c0, a, &mut c);
+            c1 = adc(c1, 0, &mut c);
+
+            *carry = c1;
+            c0
+        }
+
+        if c == 0 {
+            return;
+        }
+
+        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 if carry != 0 {
+                *ai = mac_with_carry(*ai, 0, c, &mut carry);
+            } else {
+                break;
+            }
+        }
+
+        debug_assert!(carry == 0);
+    }
+
+    // The Montgomery reduction here is based on Algorithm 14.32 in
+    // Handbook of Applied Cryptography
+    // <http://cacr.uwaterloo.ca/hac/about/chap14.pdf>.
+
+    let mut res = [0; 2*4];
+    for (i, xi) in this.iter().enumerate() {
+        mac_digit(&mut res[i..], by, *xi);
+    }
+
+    for i in 0..4 {
+        let k = inv.wrapping_mul(res[i]);
+        mac_digit(&mut res[i..], modulus, k);
+    }
+
+    this.copy_from_slice(&res[4..]);
+}
+
 #[test]
 fn setting_bits() {
     let rng = &mut ::rand::thread_rng();
-    let modulo = U256([0xffffffff; LIMBS]);
+    let modulo = U256([0xffffffffffffffff; 4]);
 
-    let a = U256::rand(rng, &modulo);
+    let a = U256::random(rng, &modulo);
     let mut e = U256::zero();
     for (i, b) in a.bits().enumerate() {
         e.set_bit(255 - i, b);

src/arith/primitives.rs

@@ -1,130 +0,0 @@
-pub type Digit = u32;
-pub const LIMBS: usize = 8;
-type DoubleDigit = u64;
-const BITS: usize = 32;
-const BASE: DoubleDigit = 1 << BITS;
-const LO_MASK: DoubleDigit = (-1i32 as DoubleDigit) >> BITS;
-
-#[inline]
-fn get_hi(n: DoubleDigit) -> Digit {
-    (n >> BITS) as Digit
-}
-
-#[inline]
-fn get_lo(n: DoubleDigit) -> Digit {
-    (n & LO_MASK) as Digit
-}
-
-#[inline]
-fn split(n: DoubleDigit) -> (Digit, Digit) {
-    (get_hi(n), get_lo(n))
-}
-
-/// Divide by two
-#[inline]
-pub fn div2(a: &mut [Digit; LIMBS]) {
-    let mut b = 0;
-    for a in a.iter_mut().rev() {
-        let t = *a << 31;
-        *a = *a >> 1;
-        *a = *a | b;
-        b = t;
-    }
-}
-
-#[inline]
-pub fn add_nocarry(a: &mut [Digit; LIMBS], b: &[Digit; LIMBS]) {
-    #[inline]
-    fn adc(a: Digit, b: Digit, carry: &mut Digit) -> Digit {
-        let (hi, lo) = split((a as DoubleDigit) + (b as DoubleDigit) +
-                             (*carry as DoubleDigit));
-
-        *carry = hi;
-        lo
-    }
-
-    let mut carry = 0;
-
-    for (a, b) in a.into_iter().zip(b.iter()) {
-        *a = adc(*a, *b, &mut carry);
-    }
-
-    assert!(0 == carry);
-}
-
-#[inline]
-pub fn sub_noborrow(a: &mut [Digit; LIMBS], b: &[Digit; LIMBS]) {
-    #[inline]
-    fn sbb(a: Digit, b: Digit, borrow: &mut Digit) -> Digit {
-        let (hi, lo) = split(BASE + (a as DoubleDigit) -
-                             (b as DoubleDigit) -
-                             (*borrow as DoubleDigit));
-        *borrow = (hi == 0) as Digit;
-        lo
-    }
-
-    let mut borrow = 0;
-
-    for (a, b) in a.into_iter().zip(b.iter()) {
-        *a = sbb(*a, *b, &mut borrow);
-    }
-
-    assert!(0 == borrow);
-}
-
-
-pub fn mul_reduce(
-    this: &mut [Digit; LIMBS],
-    by: &[Digit; LIMBS],
-    modulus: &[Digit; LIMBS],
-    inv: Digit
-)
-{
-    #[inline]
-    fn mac_digit(acc: &mut [Digit], b: &[Digit], c: Digit)
-    {
-        #[inline]
-        fn mac_with_carry(a: Digit, b: Digit, c: Digit, carry: &mut Digit) -> Digit {
-            let (hi, lo) = split((a as DoubleDigit) +
-                                 (b as DoubleDigit) * (c as DoubleDigit) +
-                                 (*carry as DoubleDigit));
-            *carry = hi;
-            lo
-        }
-
-        if c == 0 {
-            return;
-        }
-
-        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 if carry != 0 {
-                *ai = mac_with_carry(*ai, 0, c, &mut carry);
-            } else {
-                break;
-            }
-        }
-
-        assert!(carry == 0);
-    }
-
-    // The Montgomery reduction here is based on Algorithm 14.32 in
-    // Handbook of Applied Cryptography
-    // <http://cacr.uwaterloo.ca/hac/about/chap14.pdf>.
-
-    let mut res = [0; 2*LIMBS];
-    for (i, xi) in this.iter().enumerate() {
-        mac_digit(&mut res[i..], by, *xi);
-    }
-
-    for i in 0..LIMBS {
-        let k = inv.wrapping_mul(res[i]);
-        mac_digit(&mut res[i..], modulus, k);
-    }
-
-    this.copy_from_slice(&res[LIMBS..]);
-}

src/fields/fp.rs

@@ -9,7 +9,7 @@ use arith::U256;
 pub trait FpParams {
     fn name() -> &'static str;
     fn modulus() -> U256;
-    fn inv() -> u32;
+    fn inv() -> u64;
     fn rsquared() -> U256;
     fn rcubed() -> U256;
     fn one() -> U256;
@@ -69,11 +69,15 @@ impl<P: FpParams> Fp<P> {
 }
 
 impl<P: FpParams> Fp<P> {
-    /// Assumes input is mod p, not exposed publicly
+    /// Converts a U256 to an Fp so long as it's below the modulus.
-    fn new_checked(mut a: U256) -> Self {
+    pub fn new(mut a: U256) -> Option<Self> {
+        if a < P::modulus() {
             a.mul(&P::rsquared(), &P::modulus(), P::inv());
 
-        Fp(a, PhantomData)
+            Some(Fp(a, PhantomData))
+        } else {
+            None
+        }
     }
 }
 
@@ -87,20 +91,22 @@ impl<P: FpParams> FieldElement for Fp<P> {
     }
     
     fn random<R: Rng>(rng: &mut R) -> Self {
-        Fp::new_checked(U256::rand(rng, &P::modulus()))
+        const_fp(U256::random(rng, &P::modulus()))
     }
 
     fn is_zero(&self) -> bool {
         self.0.is_zero()
     }
 
-    fn inverse(mut self) -> Self {
+    fn inverse(mut self) -> Option<Self> {
-        assert!(!self.is_zero());
+        if self.is_zero() {
-
+            None
+        } else {
             self.0.invert(&P::modulus());
             self.0.mul(&P::rcubed(), &P::modulus(), P::inv());
 
-        self
+            Some(self)
+        }
     }
 }
 
@@ -152,29 +158,29 @@ impl FpParams for FrParams {
     #[inline]
     fn modulus() -> U256 {
         // 21888242871839275222246405745257275088548364400416034343698204186575808495617
-        [0xf0000001, 0x43e1f593, 0x79b97091, 0x2833e848, 0x8181585d, 0xb85045b6, 0xe131a029, 0x30644e72].into()
+        [0x43e1f593f0000001, 0x2833e84879b97091, 0xb85045b68181585d, 0x30644e72e131a029].into()
     }
 
     #[inline]
-    fn inv() -> u32 {
+    fn inv() -> u64 {
-        0xefffffff
+        0xc2e1f593efffffff
     }
 
     #[inline]
     fn rsquared() -> U256 {
         // 944936681149208446651664254269745548490766851729442924617792859073125903783
-        [0xae216da7, 0x1bb8e645, 0xe35c59e3, 0x53fe3ab1, 0x53bb8085, 0x8c49833d, 0x7f4e44a5, 0x0216d0b1].into()
+        [0x1bb8e645ae216da7, 0x53fe3ab1e35c59e3, 0x8c49833d53bb8085, 0x0216d0b17f4e44a5].into()
     }
 
     #[inline]
     fn rcubed() -> U256 {
         // 5866548545943845227489894872040244720403868105578784105281690076696998248512
-        [0xb4bf0040, 0x5e94d8e1, 0x1cfbb6b8, 0x2a489cbe, 0xa19fcfed, 0x893cc664, 0x7fcc657c, 0x0cf8594b].into()
+        [0x5e94d8e1b4bf0040, 0x2a489cbe1cfbb6b8, 0x893cc664a19fcfed, 0x0cf8594b7fcc657c].into()
     }
 
     #[inline]
     fn one() -> U256 {
-        [0x4ffffffb, 0xac96341c, 0x9f60cd29, 0x36fc7695, 0x7879462e, 0x666ea36f, 0x9a07df2f, 0x0e0a77c1].into()
+        [0xac96341c4ffffffb, 0x36fc76959f60cd29, 0x666ea36f7879462e, 0xe0a77c19a07df2f].into()
     }
 }
 
@@ -186,29 +192,29 @@ impl FpParams for FqParams {
     #[inline]
     fn modulus() -> U256 {
         // 21888242871839275222246405745257275088696311157297823662689037894645226208583
-        [0xd87cfd47, 0x3c208c16, 0x6871ca8d, 0x97816a91, 0x8181585d, 0xb85045b6, 0xe131a029, 0x30644e72].into()
+        [0x3c208c16d87cfd47, 0x97816a916871ca8d, 0xb85045b68181585d, 0x30644e72e131a029].into()
     }
 
     #[inline]
-    fn inv() -> u32 {
+    fn inv() -> u64 {
-        0xe4866389
+        0x87d20782e4866389
     }
 
     #[inline]
     fn rsquared() -> U256 {
         // 3096616502983703923843567936837374451735540968419076528771170197431451843209
-        [0x538afa89, 0xf32cfc5b, 0xd44501fb, 0xb5e71911, 0x0a417ff6, 0x47ab1eff, 0xcab8351f, 0x06d89f71].into()
+        [0xf32cfc5b538afa89, 0xb5e71911d44501fb, 0x47ab1eff0a417ff6, 0x06d89f71cab8351f].into()
     }
 
     #[inline]
     fn rcubed() -> U256 {
         // 14921786541159648185948152738563080959093619838510245177710943249661917737183
-        [0xda1530df, 0xb1cd6daf, 0xa7283db6, 0x62f210e6, 0x0ada0afb, 0xef7f0b0c, 0x2d592544, 0x20fd6e90].into()
+        [0xb1cd6dafda1530df, 0x62f210e6a7283db6, 0xef7f0b0c0ada0afb, 0x20fd6e902d592544].into()
     }
 
     #[inline]
     fn one() -> U256 {
-        [0xc58f0d9d, 0xd35d438d, 0xf5c70b3d, 0x0a78eb28, 0x7879462c, 0x666ea36f, 0x9a07df2f, 0x0e0a77c1].into()
+        [0xd35d438dc58f0d9d, 0xa78eb28f5c70b3d, 0x666ea36f7879462c, 0xe0a77c19a07df2f].into()
     }
 }
 
@@ -219,7 +225,7 @@ fn test_rsquared() {
     for _ in 0..1000 {
         let a = Fr::random(rng);
         let b: U256 = a.into();
-        let c = Fr::new_checked(b);
+        let c = Fr::new(b).unwrap();
 
         assert_eq!(a, c);
     }
@@ -227,7 +233,7 @@ fn test_rsquared() {
     for _ in 0..1000 {
         let a = Fq::random(rng);
         let b: U256 = a.into();
-        let c = Fq::new_checked(b);
+        let c = Fq::new(b).unwrap();
 
         assert_eq!(a, c);
     }

src/fields/fq2.rs

@@ -6,7 +6,7 @@ use rand::Rng;
 fn non_residue() -> Fq {
     // (q - 1) is a quadratic nonresidue in Fq
     // 21888242871839275222246405745257275088696311157297823662689037894645226208582
-    const_fp([317583274, 1757628553, 1923792719, 2366144360, 151523889, 1373741639, 1193918714, 576313009])
+    const_fp([0x68c3488912edefaa, 0x8d087f6872aabf4f, 0x51e1a24709081231, 0x2259d6b14729c0fa])
 }
 
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]
@@ -63,15 +63,16 @@ impl FieldElement for Fq2 {
         }
     }
 
-    fn inverse(self) -> Self {
+    fn inverse(self) -> Option<Self> {
         // "High-Speed Software Implementation of the Optimal Ate Pairing
         // over Barreto–Naehrig Curves"; Algorithm 8
 
-        let t = (self.c0.squared() - (self.c1.squared() * non_residue())).inverse();
+        match (self.c0.squared() - (self.c1.squared() * non_residue())).inverse() {
-
+            Some(t) => Some(Fq2 {
-        Fq2 {
                 c0: self.c0 * t,
                 c1: -(self.c1 * t)
+            }),
+            None => None
         }
     }
 }

src/fields/mod.rs

@@ -27,7 +27,7 @@ pub trait FieldElement: Sized
     fn squared(&self) -> Self {
         (*self) * (*self)
     }
-    fn inverse(self) -> Self;
+    fn inverse(self) -> Option<Self>;
     fn pow<I: Into<U256>>(&self, by: I) -> Self {
         let mut res = Self::one();
 
@@ -48,21 +48,11 @@ mod tests;
 #[test]
 fn test_fr() {
     tests::field_trials::<Fr>();
-
-    assert_eq!(
-        const_fp::<self::fp::FrParams, _>([0xf0000000, 0x43e1f593, 0x79b97091, 0x2833e848, 0x8181585d, 0xb85045b6, 0xe131a029, 0x30644e72]).inverse(),
-        const_fp::<self::fp::FrParams, _>([1105105498, 673779534, 2522683054, 3560287638, 767940567, 738640505, 1642290052, 776830401])
-    );
 }
 
 #[test]
 fn test_fq() {
     tests::field_trials::<Fq>();
-
-    assert_eq!(
-        const_fp::<self::fp::FqParams, _>([0xd87cfd46, 0x3c208c16, 0x6871ca8d, 0x97816a91, 0x8181585d, 0xb85045b6, 0xe131a029, 0x30644e72]).inverse(),
-        const_fp::<self::fp::FqParams, _>([2230452926, 1223921595, 2485962897, 3784987007, 2000672870, 1889871543, 377056010, 697020161])
-    );
 }
 
 #[test]
@@ -75,21 +65,3 @@ fn test_str() {
     assert_eq!(-Fr::one(), Fr::from_str("21888242871839275222246405745257275088548364400416034343698204186575808495616").unwrap());
     assert_eq!(-Fq::one(), Fq::from_str("21888242871839275222246405745257275088696311157297823662689037894645226208582").unwrap());
 }
-
-#[test]
-fn test_fr_squaring() {
-    // Test vector.
-    assert_eq!(
-        Fr::from_str("21888242871839275222246405745257275088548364400416034343698204186575808495610").unwrap().pow([0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff]),
-        const_fp([2572677806, 4136501712, 2816341579, 1079830324, 2386666211, 2796237710, 2115470140, 673230670])
-    );
-}
-
-#[test]
-fn test_fq_squaring() {
-    // Test vector
-    assert_eq!(
-        Fq::from_str("21888242871839275222246405745257275088696311157297823662689037894645226208540").unwrap().pow([0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff]),
-        const_fp([1128196772, 3636718463, 2555618142, 559034227, 1481851987, 3987127805, 2583602581, 88292310])
-    );
-}

src/fields/tests.rs

@@ -5,17 +5,33 @@ fn can_invert<F: FieldElement>() {
     let mut a = F::one();
 
     for _ in 0..10000 {
-        assert_eq!(a * a.inverse(), F::one());
+        assert_eq!(a * a.inverse().unwrap(), F::one());
 
         a = a + F::one();
     }
 
     a = -F::one();
     for _ in 0..10000 {
-        assert_eq!(a * a.inverse(), F::one());
+        assert_eq!(a * a.inverse().unwrap(), F::one());
 
         a = a - F::one();
     }
+
+    assert_eq!(F::zero().inverse(), None);
+}
+
+fn rand_element_eval<F: FieldElement, R: Rng>(rng: &mut R) {
+    for _ in 0..100 {
+        let a = F::random(rng);
+        let b = F::random(rng);
+        let c = F::random(rng);
+        let d = F::random(rng);
+
+        assert_eq!(
+            (a + b) * (c + d),
+            (a * c) + (b * c) + (a * d) + (b * d)
+        );
+    }
 }
 
 fn rand_element_squaring<F: FieldElement, R: Rng>(rng: &mut R) {
@@ -34,6 +50,12 @@ fn rand_element_squaring<F: FieldElement, R: Rng>(rng: &mut R) {
 }
 
 fn rand_element_addition_and_negation<F: FieldElement, R: Rng>(rng: &mut R) {
+    for _ in 0..100 {
+        let a = F::random(rng);
+
+        assert_eq!(a + (-a), F::zero());
+    }
+
     for _ in 0..100 {
         let mut a = F::random(rng);
         let r = F::random(rng);
@@ -69,9 +91,9 @@ fn rand_element_addition_and_negation<F: FieldElement, R: Rng>(rng: &mut R) {
 fn rand_element_inverse<F: FieldElement, R: Rng>(rng: &mut R) {
     for _ in 0..10000 {
         let a = F::random(rng);
-        assert!(a.inverse() * a == F::one());
+        assert!(a.inverse().unwrap() * a == F::one());
         let b = F::random(rng);
-        assert_eq!((a * b) * (a.inverse()), b);
+        assert_eq!((a * b) * (a.inverse().unwrap()), b);
     }
 }
 
@@ -102,4 +124,5 @@ pub fn field_trials<F: FieldElement>() {
     rand_element_addition_and_negation::<F, StdRng>(&mut rng);
     rand_element_multiplication::<F, StdRng>(&mut rng);
     rand_element_inverse::<F, StdRng>(&mut rng);
+    rand_element_eval::<F, StdRng>(&mut rng);
 }

src/groups/mod.rs

@@ -0,0 +1,258 @@
+use std::ops::{Add,Sub,Neg,Mul};
+use fields::{FieldElement, Fq, Fq2, Fr, const_fp};
+use arith::U256;
+use std::fmt;
+use rand::Rng;
+
+pub trait GroupElement: Sized +
+                    Copy +
+                    Clone +
+                    PartialEq +
+                    Eq +
+                    fmt::Debug +
+                    Add<Output=Self> +
+                    Sub<Output=Self> +
+                    Neg<Output=Self> +
+                    Mul<Fr, Output=Self>
+{
+    fn zero() -> Self;
+    fn one() -> Self;
+    fn random<R: Rng>(rng: &mut R) -> Self;
+    fn is_zero(&self) -> bool;
+    fn double(&self) -> Self;
+}
+
+pub trait GroupParams: Sized {
+    type Base: FieldElement;
+
+    fn name() -> &'static str;
+    fn one() -> Jacobian<Self>;
+}
+
+pub struct Jacobian<P: GroupParams> {
+    x: P::Base,
+    y: P::Base,
+    z: P::Base
+}
+
+impl<P: GroupParams> fmt::Debug for Jacobian<P> {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        write!(f, "{}({:?}, {:?}, {:?})", P::name(), self.x, self.y, self.z)
+    }
+}
+
+impl<P: GroupParams> Clone for Jacobian<P> {
+    fn clone(&self) -> Self {
+        Jacobian {
+            x: self.x,
+            y: self.y,
+            z: self.z
+        }
+    }
+}
+impl<P: GroupParams> Copy for Jacobian<P> {}
+
+impl<P: GroupParams> PartialEq for Jacobian<P> {
+    fn eq(&self, other: &Self) -> bool {
+        if self.is_zero() {
+            return other.is_zero()
+        }
+
+        if other.is_zero() {
+            return false;
+        }
+
+        let z1_squared = self.z.squared();
+        let z2_squared = other.z.squared();
+
+        if self.x * z2_squared != other.x * z1_squared {
+            return false;
+        }
+
+        let z1_cubed = self.z * z1_squared;
+        let z2_cubed = other.z * z2_squared;
+
+        if self.y * z2_cubed != other.y * z1_cubed {
+            return false;
+        }
+
+        return true;
+    }
+}
+impl<P: GroupParams> Eq for Jacobian<P> { }
+
+impl<P: GroupParams> GroupElement for Jacobian<P> {
+    fn zero() -> Self {
+        Jacobian {
+            x: P::Base::zero(),
+            y: P::Base::one(),
+            z: P::Base::zero()
+        }
+    }
+
+    fn one() -> Self {
+        P::one()
+    }
+
+    fn random<R: Rng>(rng: &mut R) -> Self {
+        P::one() * Fr::random(rng)
+    }
+
+    fn is_zero(&self) -> bool {
+        self.z.is_zero()
+    }
+
+    fn double(&self) -> Self {
+        let a = self.x.squared();
+        let b = self.y.squared();
+        let c = b.squared();
+        let mut d = (self.x + b).squared() - a - c;
+        d = d + d;
+        let e = a + a + a;
+        let f = e.squared();
+        let x3 = f - (d + d);
+        let mut eight_c = c + c;
+        eight_c = eight_c + eight_c;
+        eight_c = eight_c + eight_c;
+        let y1z1 = self.y * self.z;
+
+        Jacobian {
+            x: x3,
+            y: e * (d - x3) - eight_c,
+            z: y1z1 + y1z1
+        }
+    }
+}
+
+impl<P: GroupParams> Mul<Fr> for Jacobian<P> {
+    type Output = Jacobian<P>;
+
+    fn mul(self, other: Fr) -> Jacobian<P> {
+        let mut res = Jacobian::zero();
+        let mut found_one = false;
+
+        for i in U256::from(other).bits() {
+            if found_one {
+                res = res.double();
+            }
+
+            if i {
+                found_one = true;
+                res = res + self;
+            }
+        }
+
+        res
+    }
+}
+
+impl<P: GroupParams> Add<Jacobian<P>> for Jacobian<P> {
+    type Output = Jacobian<P>;
+
+    fn add(self, other: Jacobian<P>) -> Jacobian<P> {
+        if self.is_zero() {
+            return other;
+        }
+
+        if other.is_zero() {
+            return self;
+        }
+
+        let z1_squared = self.z.squared();
+        let z2_squared = other.z.squared();
+        let u1 = self.x * z2_squared;
+        let u2 = other.x * z1_squared;
+        let z1_cubed = self.z * z1_squared;
+        let z2_cubed = other.z * z2_squared;
+        let s1 = self.y * z2_cubed;
+        let s2 = other.y * z1_cubed;
+
+        if u1 == u2 && s1 == s2 {
+            self.double()
+        } else {
+            let h = u2 - u1;
+            let s2_minus_s1 = s2 - s1;
+            let i = (h + h).squared();
+            let j = h * i;
+            let r = s2_minus_s1 + s2_minus_s1;
+            let v = u1 * i;
+            let s1_j = s1 * j;
+            let x3 = r.squared() - j - (v + v);
+
+            Jacobian {
+                x: x3,
+                y: r * (v - x3) - (s1_j + s1_j),
+                z: ((self.z + other.z).squared() - z1_squared - z2_squared) * h
+            }
+        }
+    }
+}
+
+impl<P: GroupParams> Neg for Jacobian<P> {
+    type Output = Jacobian<P>;
+
+    fn neg(self) -> Jacobian<P> {
+        Jacobian {
+            x: self.x,
+            y: -self.y,
+            z: self.z
+        }
+    }
+}
+
+impl<P: GroupParams> Sub<Jacobian<P>> for Jacobian<P> {
+    type Output = Jacobian<P>;
+
+    fn sub(self, other: Jacobian<P>) -> Jacobian<P> {
+        self + (-other)
+    }
+}
+
+pub struct G1Params;
+
+impl GroupParams for G1Params {
+    type Base = Fq;
+
+    fn name() -> &'static str { "G1" }
+
+    fn one() -> Jacobian<Self> {
+        Jacobian {
+            x: Fq::one(),
+            y: const_fp([0xa6ba871b8b1e1b3a, 0x14f1d651eb8e167b, 0xccdd46def0f28c58, 0x1c14ef83340fbe5e]),
+            z: Fq::one()
+        }
+    }
+}
+
+pub type G1 = Jacobian<G1Params>;
+
+pub struct G2Params;
+
+impl GroupParams for G2Params {
+    type Base = Fq2;
+
+    fn name() -> &'static str { "G2" }
+
+    fn one() -> Jacobian<Self> {
+        Jacobian {
+            x: Fq2::new(const_fp([0x8e83b5d102bc2026, 0xdceb1935497b0172, 0xfbb8264797811adf, 0x19573841af96503b]), const_fp([0xafb4737da84c6140, 0x6043dd5a5802d8c4, 0x09e950fc52a02f86, 0x14fef0833aea7b6b])),
+            y: Fq2::new(const_fp([0x619dfa9d886be9f6, 0xfe7fd297f59e9b78, 0xff9e1a62231b7dfe, 0x28fd7eebae9e4206]), const_fp([0x64095b56c71856ee, 0xdc57f922327d3cbb, 0x55f935be33351076, 0x0da4a0e693fd6482])),
+            z: Fq2::one()
+        }
+    }
+}
+
+pub type G2 = Jacobian<G2Params>;
+
+#[cfg(test)]
+mod tests;
+
+#[test]
+fn test_g1() {
+    tests::group_trials::<G1>();
+}
+
+#[test]
+fn test_g2() {
+    tests::group_trials::<G2>();
+}

src/groups/tests.rs

@@ -0,0 +1,100 @@
+use super::GroupElement;
+use fields::{FieldElement, Fr};
+use rand::Rng;
+
+fn random_test_addition<G: GroupElement, R: Rng>(rng: &mut R) {
+    for _ in 0..50 {
+        let r1 = G::random(rng);
+        let r2 = G::random(rng);
+        let r3 = G::random(rng);
+
+        assert_eq!((r1 + r2) + r3, r1 + (r2 + r3));
+        assert!(((r1 + r2 + r3) - r2 - r3 - r1).is_zero());
+    }
+}
+
+fn random_test_doubling<G: GroupElement, R: Rng>(rng: &mut R) {
+    for _ in 0..50 {
+        let r1 = G::random(rng);
+        let r2 = G::random(rng);
+        let ti = Fr::from_str("2").unwrap().inverse().unwrap();
+
+        assert_eq!((r1 + r2) + r1, r1.double() + r2);
+        assert_eq!(r1, r1.double() * ti);
+    }
+}
+
+fn random_test_dh<G: GroupElement, R: Rng>(rng: &mut R) {
+    for _ in 0..50 {
+        let alice_sk = Fr::random(rng);
+        let bob_sk = Fr::random(rng);
+
+        let alice_pk = G::one() * alice_sk;
+        let bob_pk = G::one() * bob_sk;
+
+        let alice_shared = bob_pk * alice_sk;
+        let bob_shared = alice_pk * bob_sk;
+
+        assert_eq!(alice_shared, bob_shared);
+    }
+}
+
+fn random_test_equality<G: GroupElement, R: Rng>(rng: &mut R) {
+    for _ in 0..50 {
+        let begin = G::random(rng);
+
+        let mut acc = begin;
+
+        // Do a bunch of random things.
+
+        let a = Fr::random(rng);
+        let b = G::random(rng);
+        let c = Fr::random(rng);
+        let d = G::random(rng);
+
+        for _ in 0..10 {
+            acc = acc * a;
+            acc = -acc;
+            acc = acc + b;
+            acc = acc * c;
+            acc = -acc;
+            acc = acc - d;
+            acc = acc.double();
+        }
+
+        // Then reverse the operations
+
+        let ai = a.inverse().unwrap();
+        let ci = c.inverse().unwrap();
+        let ti = Fr::from_str("2").unwrap().inverse().unwrap();
+
+        for _ in 0..10 {
+            acc = acc * ti;
+            acc = acc + d;
+            acc = -acc;
+            acc = acc * ci;
+            acc = acc - b;
+            acc = -acc;
+            acc = acc * ai;
+        }
+
+        assert_eq!(acc, begin);
+    }
+}
+
+pub fn group_trials<G: GroupElement>() {
+    assert!(G::zero().is_zero());
+    assert!((G::one() - G::one()).is_zero());
+    assert_eq!(G::one() + G::one(), G::one() * Fr::from_str("2").unwrap());
+    assert!(G::zero().double().is_zero());
+
+    use rand::{SeedableRng,StdRng};
+    let seed: [usize; 4] = [103245, 191922, 1293, 192103];
+    let mut rng = StdRng::from_seed(&seed);
+
+    random_test_addition::<G, _>(&mut rng);
+    random_test_doubling::<G, _>(&mut rng);
+    random_test_dh::<G, _>(&mut rng);
+    random_test_equality::<G, _>(&mut rng);
+}
+

src/lib.rs

@@ -2,6 +2,8 @@ extern crate rand;
 
 mod arith;
 mod fields;
+mod groups;
 
 pub use arith::U256;
 pub use fields::{Fq, Fr, Fq2, FieldElement};
+pub use groups::{G1, G2, GroupElement};