Implement point doubling for weierstrass curves

src/weierstrass/point.rs

@@ -139,12 +139,61 @@ impl Point {
             z: z3,
         })
     }
+    /// core_doubling implements exception free point doubling formulas for prime order groups.
+    // Reference: "Complete addition formulas for prime order elliptic curves" by
+    // Costello-Renes-Batina. [Alg.3] (eprint.iacr.org/2015/1060).
+    fn core_doubling(&self) -> <Curve as EllipticCurve>::Point {
+        let a = &self.e.a;
+        let b3 = &self.e.b + &self.e.b + &self.e.b;
+        let x = &self.c.x;
+        let y = &self.c.y;
+        let z = &self.c.z;
+        let (mut x3, mut y3, mut z3);
+        let (mut t0, t1, mut t2, mut t3);
+        t0 = x * x; //      1. t0 =  X *  X
+        t1 = y * y; //      2. t1 =  Y *  Y
+        t2 = z * z; //      3. t2 =  Z *  Z
+        t3 = x * y; //      4. t3 =  X *  Y
+        t3 = &t3 + &t3; //  5. t3 = t3 + t3
+        z3 = x * z; //      6. Z3 =  X *  Z
+        z3 = &z3 + &z3; //  7. Z3 = Z3 + Z3
+        x3 = a * &z3; //    8. X3 =  a * Z3
+        y3 = &b3 * &t2; //  9. Y3 = b3 * t2
+        y3 = &x3 + &y3; // 10. Y3 = X3 + Y3
+        x3 = &t1 - &y3; // 11. X3 = t1 - Y3
+        y3 = &t1 + &y3; // 12. Y3 = t1 + Y3
+        y3 = &x3 * &y3; // 13. Y3 = X3 * Y3
+        x3 = &t3 * &x3; // 14. X3 = t3 * X3
+        z3 = &b3 * &z3; // 15. Z3 = b3 * Z3
+        t2 = a * &t2; //   16. t2 =  a * t2
+        t3 = &t0 - &t2; // 17. t3 = t0 - t2
+        t3 = a * &t3; //   18. t3 =  a * t3
+        t3 = &t3 + &z3; // 19. t3 = t3 + Z3
+        z3 = &t0 + &t0; // 20. Z3 = t0 + t0
+        t0 = &z3 + &t0; // 21. t0 = Z3 + t0
+        t0 = &t0 + &t2; // 22. t0 = t0 + t2
+        t0 = &t0 * &t3; // 23. t0 = t0 * t3
+        y3 = &y3 + &t0; // 24. Y3 = Y3 + t0
+        t2 = y * z; //     25. t2 =  Y *  Z
+        t2 = &t2 + &t2; // 26. t2 = t2 + t2
+        t0 = &t2 * &t3; // 27. t2 = t2 * t3
+        x3 = &x3 - &t0; // 28. X3 = X3 - t0
+        z3 = &t2 * &t1; // 29. Z3 = t2 * t1
+        z3 = &z3 + &z3; // 30. Z3 = Z3 + Z3
+        z3 = &z3 + &z3; // 31. Z3 = Z3 + Z3
+        self.e.new_proy_point(ProyCoordinates {
+            x: x3,
+            y: y3,
+            z: z3,
+        })
+    }
+
     /// core_mul implements the double&add Scalar multiplication method.
     /// This function run in non-constant time.
     fn core_mul(&self, k: &Scalar) -> <Curve as EllipticCurve>::Point {
         let mut q = self.e.identity();
         for ki in k.iter_lr() {
-            q = &q + &q;
+            q = q.core_doubling();
             if ki {
                 q = q + self;
             }