add jacobian algo

src/lib.rs

@@ -22,6 +22,15 @@ struct Point {
     curve: CurveFp,
 }
 
+/// 椭圆曲线上的点（雅可比坐标）
+#[derive(Clone)]
+struct JacobianPoint {
+    x: BigUint,     // X 坐标
+    y: BigUint,     // Y 坐标
+    z: BigUint,     // Z 坐标
+    curve: CurveFp, // 椭圆曲线的参数
+}
+
 // Initialize the SM2 Curve
 fn sm2p256v1() -> CurveFp {
     CurveFp {
@@ -139,10 +148,8 @@ fn mod_inverse(a: &BigUint, m: &BigUint) -> BigUint {
 
     t
 }
-fn point_addition(p1: &Point, p2: &Point) -> Point {
-    let curve = &p1.curve;
-    let p = &curve.p;
 
+fn point_addition(p1: &Point, p2: &Point) -> Point {
     if p1.x.is_zero() && p1.y.is_zero() {
         return p2.clone();
     }
@@ -150,46 +157,162 @@ fn point_addition(p1: &Point, p2: &Point) -> Point {
         return p1.clone();
     }
 
-    let lambda = if p1.x == p2.x && p1.y == p2.y {
+    from_jacobian(jacobian_add(to_jacobian(p1), to_jacobian(p2)))
-        let num = (BigUint::from(3u32) * &p1.x * &p1.x + &curve.a) % p;
+}
-        let denom = (BigUint::from(2u32) * &p1.y) % p;
-        (num * mod_inverse(&denom, p)) % p
-    } else {
-        let num = ((&p2.y + p) - &p1.y) % p;
-        let denom = ((&p2.x + p) - &p1.x) % p;
 
-        (num * mod_inverse(&denom, p)) % p
+/// 将仿射坐标转换为雅可比坐标 (X, Y, Z)
+fn to_jacobian(p: &Point) -> JacobianPoint {
+    JacobianPoint {
+        x: p.x.clone(),
+        y: p.y.clone(),
+        z: BigUint::one(), // Z = 1 表示仿射坐标
+        curve: sm2p256v1(),
+    }
+}
+
+/// 将雅可比坐标转换为仿射坐标
+fn from_jacobian(p: JacobianPoint) -> Point {
+    if p.z.is_zero() {
+        return Point {
+            x: BigUint::zero(),
+            y: BigUint::zero(),
+            curve: sm2p256v1(),
         };
+    }
+    let p_mod = &p.curve.p;
 
-    let x3 = (lambda.clone() * &lambda - &p1.x - &p2.x) % p;
+    // 计算 Z 的模反
-    let y3 = (lambda * (&p1.x + p - &x3) - &p1.y) % p;
+    let z_inv = mod_inverse(&p.z, p_mod);
+    let z_inv2 = (&z_inv * &z_inv) % p_mod; // Z_inv^2
+    let z_inv3 = (&z_inv2 * &z_inv) % p_mod; // Z_inv^3
+
+    // 计算 x = X * Z_inv^2, y = Y * Z_inv^3
+    let x_affine = (&p.x * &z_inv2) % p_mod;
+    let y_affine = (&p.y * &z_inv3) % p_mod;
 
     Point {
+        x: x_affine,
+        y: y_affine,
+        curve: sm2p256v1(),
+    }
+}
+
+/// 雅可比坐标下的点加法
+fn jacobian_add(p1: JacobianPoint, p2: JacobianPoint) -> JacobianPoint {
+    if p1.z.is_zero() {
+        return p2.clone();
+    }
+    if p2.z.is_zero() {
+        return p1.clone();
+    }
+
+    let p_mod = &p1.curve.p;
+
+    // U1 = X1 * Z2^2, U2 = X2 * Z1^2
+    let z1z1 = (&p1.z * &p1.z) % p_mod;
+    let z2z2 = (&p2.z * &p2.z) % p_mod;
+    let u1 = (&p1.x * &z2z2) % p_mod;
+    let u2 = (&p2.x * &z1z1) % p_mod;
+
+    // S1 = Y1 * Z2^3, S2 = Y2 * Z1^3
+    let z1z1z1 = (&z1z1 * &p1.z) % p_mod;
+    let z2z2z2 = (&z2z2 * &p2.z) % p_mod;
+    let s1 = (&p1.y * &z2z2z2) % p_mod;
+    let s2 = (&p2.y * &z1z1z1) % p_mod;
+
+    if u1 == u2 && s1 == s2 {
+        // 点倍运算 (p1 == p2)
+        return jacobian_double(p1);
+    }
+
+    // H = U2 - U1, R = S2 - S1
+    let h = (u2 + p_mod - &u1) % p_mod;
+    let r = (s2 + p_mod - &s1) % p_mod;
+
+    // X3 = R^2 - H^3 - 2 * U1 * H^2
+    let h2 = (&h * &h) % p_mod;
+    let h3 = (&h2 * &h) % p_mod;
+    let x3 = (&r * &r + p_mod - &h3 - (&BigUint::from(2u32) * &u1 * &h2) % p_mod) % p_mod;
+
+    // Y3 = R * (U1 * H^2 - X3) - S1 * H^3
+    let y3 =
+        (&r * ((&u1 * &h2 + p_mod - (&x3 % p_mod)) % p_mod) + p_mod - (&s1 * &h3) % p_mod) % p_mod;
+
+    // Z3 = Z1 * Z2 * H
+    let z3 = (&p1.z * &p2.z * &h) % p_mod;
+
+    JacobianPoint {
         x: x3,
         y: y3,
+        z: z3,
+        curve: sm2p256v1(),
+    }
+}
+
+/// 雅可比坐标下的点倍运算
+fn jacobian_double(p: JacobianPoint) -> JacobianPoint {
+    let curve = &p.curve;
+    let p_mod = &curve.p;
+
+    if p.y.is_zero() {
+        return JacobianPoint {
+            x: BigUint::zero(),
+            y: BigUint::zero(),
+            z: BigUint::zero(),
+            curve: curve.clone(),
+        };
+    }
+
+    // S = 4 * X * Y^2
+    let y2 = (&p.y * &p.y) % p_mod;
+    let s = (&p.x * &y2 * BigUint::from(4u32)) % p_mod;
+
+    // M = 3 * X^2 + a * Z^4
+    let z2 = (&p.z * &p.z) % p_mod;
+    let z4 = (&z2 * &z2) % p_mod;
+    let m = ((&p.x * &p.x * BigUint::from(3u32)) + &curve.a * &z4) % p_mod;
+
+    // X3 = M^2 - 2 * S
+    let x3 = (&m * &m + p_mod - &s * BigUint::from(2u32)) % p_mod;
+
+    // Y3 = M * (S - X3) - 8 * Y^4
+    let y4 = (&y2 * &y2) % p_mod;
+    let y3 = (&m * (&s + p_mod - &x3) + p_mod - BigUint::from(8u32) * &y4) % p_mod;
+
+    // Z3 = 2 * Y * Z
+    let z3 = (&p.y * &p.z * BigUint::from(2u32)) % p_mod;
+
+    JacobianPoint {
+        x: x3,
+        y: y3,
+        z: z3,
         curve: curve.clone(),
     }
 }
 
 fn point_multiplication(p: &Point, n: &BigUint) -> Point {
-    let mut result = Point {
+    let mut result = JacobianPoint {
         x: BigUint::zero(),
-        y: BigUint::zero(),
+        y: BigUint::one(), // 无穷远点的雅可比坐标表示
+        z: BigUint::zero(),
         curve: p.curve.clone(),
     };
 
-    let mut addend = p.clone();
+    // 将输入点从仿射坐标转换为雅可比坐标
+    let mut addend = to_jacobian(p);
     let mut k = n.clone();
 
+    // 使用二进制展开法进行点乘运算
     while !k.is_zero() {
         if &k % 2u32 == BigUint::one() {
-            result = point_addition(&result, &addend);
+            result = jacobian_add(result, addend.clone());
         }
-        addend = point_addition(&addend, &addend);
+        addend = jacobian_double(addend); // 倍点运算
         k >>= 1;
     }
 
-    result
+    // 将结果从雅可比坐标转换为仿射坐标
+    from_jacobian(result)
 }
 
 /// SM2 addition
@@ -213,13 +336,13 @@ fn add(p1: (BigUint, BigUint), p2: (BigUint, BigUint)) -> (BigUint, BigUint) {
     let point1 = Point {
         x: x1,
         y: y1,
-        curve: curve.clone(),
+        curve: sm2p256v1(),
     };
 
     let point2 = Point {
         x: x2,
         y: y2,
-        curve: curve.clone(),
+        curve: sm2p256v1(),
     };
 
     let result = point_addition(&point1, &point2);
@@ -232,9 +355,9 @@ fn multiply(point: (BigUint, BigUint), n: BigUint) -> (BigUint, BigUint) {
     let curve = sm2p256v1();
     // Construct the point with BigUint values
     let point = Point {
-        x: point.0.clone(),
+        x: point.0,
-        y: point.1.clone(),
+        y: point.1,
-        curve: curve.clone(),
+        curve,
     };
     // Perform point multiplication
     let result = point_multiplication(&point, &n);