Bw6 miller loop optimization

ec/src/models/bw6/g2.rs

@@ -46,6 +46,16 @@ impl<P: BW6Config> Default for G2Prepared<P> {
     }
 }
 
+// impl into G2Affine from G2HomProjective
+impl<P: BW6Config> From<G2HomProjective<P>> for G2Affine<P> {
+    fn from(q: G2HomProjective<P>) -> Self {
+        let z_inv = q.z.inverse().unwrap();
+        let x = q.x * &z_inv;
+        let y = q.y * &z_inv;
+        G2Affine::<P>::new_unchecked(x, y)
+    }
+}
+
 impl<P: BW6Config> From<G2Affine<P>> for G2Prepared<P> {
     fn from(q: G2Affine<P>) -> Self {
         if q.infinity {
@@ -56,7 +66,7 @@ impl<P: BW6Config> From<G2Affine<P>> for G2Prepared<P> {
             };
         }
 
-        // f_{u+1,Q}(P)
+        // f_{u,Q}(P)
         let mut ell_coeffs_1 = vec![];
         let mut r = G2HomProjective::<P> {
             x: q.x,
@@ -71,23 +81,34 @@ impl<P: BW6Config> From<G2Affine<P>> for G2Prepared<P> {
                 ell_coeffs_1.push(r.add_in_place(&q));
             }
         }
+        // TODO: this is probably the slowest part
+        // While G2 preparation is overall faster due to shortened 2nd loop,
+        // The inversion could probably be avoided by using Hom(P) + Hom(Q) addition,
+        // instead of mixed addition as is currently done.
+        let r_affine: G2Affine<P> = r.into();
+        // Swap the signs of `qu`, `r` & `neg_qu` if the loop count is negative.
+        let (qu, neg_qu) = if P::ATE_LOOP_COUNT_1_IS_NEGATIVE {
+            (-r_affine, r_affine)
+        } else {
+            (r_affine, -r_affine)
+        };
 
-        // f_{u^3-u^2-u,Q}(P)
-        let mut ell_coeffs_2 = vec![];
-        let mut r = G2HomProjective::<P> {
-            x: q.x,
-            y: q.y,
+        r = G2HomProjective::<P> {
+            x: qu.x,
+            y: qu.y,
             z: P::Fp::one(),
         };
+        ell_coeffs_1.push(r.clone().add_in_place(&q));
 
-        let negq = -q;
+        let mut ell_coeffs_2 = vec![];
 
+        // f_{u^2-u-1,[u]Q}(P)
         for bit in P::ATE_LOOP_COUNT_2.iter().rev().skip(1) {
             ell_coeffs_2.push(r.double_in_place());
 
             match bit {
-                1 => ell_coeffs_2.push(r.add_in_place(&q)),
-                -1 => ell_coeffs_2.push(r.add_in_place(&negq)),
+                1 => ell_coeffs_2.push(r.add_in_place(&qu)),
+                -1 => ell_coeffs_2.push(r.add_in_place(&neg_qu)),
                 _ => continue,
             }
         }

ec/src/models/bw6/mod.rs

@@ -30,6 +30,7 @@ pub trait BW6Config: 'static + Eq + Sized {
     const X_MINUS_1_DIV_3: <Self::Fp as PrimeField>::BigInt;
     const ATE_LOOP_COUNT_1: &'static [u64];
     const ATE_LOOP_COUNT_1_IS_NEGATIVE: bool;
+    // X^2 - X - 1
     const ATE_LOOP_COUNT_2: &'static [i8];
     const ATE_LOOP_COUNT_2_IS_NEGATIVE: bool;
     const TWIST_TYPE: TwistType;
@@ -96,7 +97,8 @@ pub trait BW6Config: 'static + Eq + Sized {
             })
             .unzip::<_, _, Vec<_>, Vec<_>>();
 
-        let mut f_1 = cfg_chunks_mut!(pairs_1, 4)
+        // compute f_u which we can later re-use for the 2nd loop
+        let mut f_u = cfg_chunks_mut!(pairs_1, 4)
             .map(|pairs| {
                 let mut f = <BW6<Self> as Pairing>::TargetField::one();
                 for i in BitIteratorBE::without_leading_zeros(Self::ATE_LOOP_COUNT_1).skip(1) {
@@ -114,26 +116,45 @@ pub trait BW6Config: 'static + Eq + Sized {
             })
             .product::<<BW6<Self> as Pairing>::TargetField>();
 
+        let f_u_inv;
+
         if Self::ATE_LOOP_COUNT_1_IS_NEGATIVE {
-            f_1.cyclotomic_inverse_in_place();
+            f_u_inv = f_u;
+            f_u.cyclotomic_inverse_in_place();
+        } else {
+            f_u_inv = f_u.cyclotomic_inverse().unwrap();
         }
+
+        // f_1(P) = f_(u+1)(P) = f_u(P) * l([u]q, q)(P)
+        let mut f_1 = cfg_chunks_mut!(pairs_1, 4)
+            .map(|pairs| {
+                pairs.iter_mut().fold(f_u, |mut f, (p, coeffs)| {
+                    BW6::<Self>::ell(&mut f, &coeffs.next().unwrap(), &p.0);
+                    f
+                })
+            })
+            .product::<<BW6<Self> as Pairing>::TargetField>();
+
         let mut f_2 = cfg_chunks_mut!(pairs_2, 4)
             .map(|pairs| {
-                let mut f = <<BW6<Self> as Pairing>::TargetField>::one();
+                let mut f = f_u;
                 for i in (1..Self::ATE_LOOP_COUNT_2.len()).rev() {
-                    if i != Self::ATE_LOOP_COUNT_2.len() - 1 {
-                        f.square_in_place();
-                    }
+                    f.square_in_place();
 
                     for (p, ref mut coeffs) in pairs.iter_mut() {
                         BW6::<Self>::ell(&mut f, &coeffs.next().unwrap(), &p.0);
                     }
 
                     let bit = Self::ATE_LOOP_COUNT_2[i - 1];
-                    if bit == 1 || bit == -1 {
-                        for &mut (p, ref mut coeffs) in pairs.iter_mut() {
-                            BW6::<Self>::ell(&mut f, &coeffs.next().unwrap(), &p.0);
-                        }
+                    if bit == 1 {
+                        f *= &f_u;
+                    } else if bit == -1 {
+                        f *= &f_u_inv;
+                    } else {
+                        continue;
+                    }
+                    for &mut (p, ref mut coeffs) in pairs.iter_mut() {
+                        BW6::<Self>::ell(&mut f, &coeffs.next().unwrap(), &p.0);
                     }
                 }
                 f