Adding shamir's trick and windowing method

src/secp256r1/ec_mulmuladd.cairo

@@ -0,0 +1,205 @@
+//*************************************************************************************/
+///* Copyright (C) 2022 - Renaud Dubois - This file is part of Cairo_musig2 project   */
+///* License: This software is licensed under a dual BSD and GPL v2 license.          */
+///* See LICENSE file at the root folder of the project.                              */
+///* FILE: multipoint.cairo                                                           */
+///*                                                                                  */
+///*                                                                                  */
+///* DESCRIPTION: optimization of dual base multiplication                            */
+///* the algorithm combines the so called Shamir's trick with Windowing method        */
+//*************************************************************************************/
+from starkware.cairo.common.cairo_secp.bigint import BigInt3
+
+from src.secp256r1.ec import EcPoint, ec_add, ec_mul, ec_double
+
+
+//Structure storing all aP+b.Q for (a,b) in [0..3]x[0..3]
+struct Window {
+    G: EcPoint,
+    Q: EcPoint,
+    W3: EcPoint,
+    W4: EcPoint,
+    W5: EcPoint,
+    W6: EcPoint,
+    W7: EcPoint,
+    W8: EcPoint,
+    W9: EcPoint,
+    W10: EcPoint,
+    W11: EcPoint,
+    W12: EcPoint,
+    W13: EcPoint,
+    W14: EcPoint,
+    W15: EcPoint,
+}
+
+
+//https://crypto.stackexchange.com/questions/99975/strauss-shamir-trick-on-ec-multiplication-by-scalar,
+//* Internal call for recursion of point multiplication via Shamir's trick */
+func ec_mulmuladd_inner{range_check_ptr}(
+    R: EcPoint, G: EcPoint, Q: EcPoint, H: EcPoint,
+    scalar_u: felt, scalar_v: felt, m: felt
+) -> (res: EcPoint){
+    alloc_locals;
+
+      // this means if m=-1, beware if felt definition changes
+    if (m == -1) {
+        return(res=R);
+    }
+
+    let (double_point) = ec_double(R);
+
+    let mm1 = m-1;
+    local dibit;
+    //extract MSB values of both exponents
+     %{ ids.dibit = ((ids.scalar_u >>ids.m)&1) +2*((ids.scalar_v >>ids.m)&1) %}
+
+    //set R:=R+R
+    if (dibit==0){
+      let (res)= ec_mulmuladd_inner(double_point,G,Q,H,scalar_u,scalar_v,mm1);
+      return(res=res);
+    }
+    //if ui=1 and vi=0, set R:=R+G
+    if (dibit==1){
+      let (res10)=ec_add(double_point,G);
+      let (res)= ec_mulmuladd_inner(res10,G,Q,H,scalar_u,scalar_v,mm1);
+      return(res=res);
+    }
+    //(else) if ui=0 and vi=1, set R:=R+Q
+    if (dibit==2){
+      let (res01)=ec_add(double_point,Q);
+      let (res)= ec_mulmuladd_inner(res01,G,Q,H,scalar_u,scalar_v,mm1);
+      return(res=res);
+        }
+    //(else) if ui=1 and vi=1, set R:=R+Q
+    if (dibit==3){
+     let (res11)=ec_add(double_point,H);
+      let (res)= ec_mulmuladd_inner(res11,G,Q,H,scalar_u,scalar_v,mm1);
+      return(res=res);
+    }
+
+     //you shall never end up here
+     return(res=R);
+}
+
+
+
+//https://crypto.stackexchange.com/questions/99975/strauss-shamir-trick-on-ec-multiplication-by-scalar,
+//* Internal call for recursion of point multiplication via Shamir's trick+Windowed method */
+func ec_mulmuladd_W_inner{range_check_ptr}(
+    R: EcPoint, Prec:Window,
+    scalar_u: felt, scalar_v: felt, m: felt
+) -> (res: EcPoint){
+    alloc_locals;
+    let mm2 = m-2;
+
+    //(8*v1 4*u1+ 2*v0 + u0), where (u1,u0) represents two bit at index m of scalar u, (resp for v)
+    local quad_bit;
+
+    if (m == -1) {
+        return (res=R);
+    }
+
+    let (double_point) = ec_double(R);
+
+     //still have to make the last addition over 1 bit (initial length was odd)
+    if(m == 0){
+     let (res)=ec_mulmuladd_inner(R, Prec.G, Prec.Q, Prec.W3, scalar_u, scalar_v, m);
+     return (res=res);
+    }
+
+    let (quadruple_point) = ec_double(double_point);
+
+    //compute quadruple (8*v1 4*u1+ 2*v0 + u0)
+    %{
+        ids.quad_bit = (
+            8 * ((ids.scalar_v >> ids.m) & 1)
+            + 4 * ((ids.scalar_u >> ids.m) & 1)
+            + 2 * ((ids.scalar_v >> (ids.m - 1)) & 1)
+            + ((ids.scalar_u >> (ids.m - 1)) & 1)
+        )
+    %}
+
+    if (quad_bit == 0) {
+        let (res) = ec_mulmuladd_W_inner(quadruple_point, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 1) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.G);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 2) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.Q);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+
+    if (quad_bit == 3) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W3);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 4) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W4);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 5) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W5);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 6) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W6);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 7) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W7);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 8) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W8);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 9) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W9);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 10) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W10);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 11) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W11);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 12) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W12);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 13) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W13);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 14) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W14);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+    if (quad_bit == 15) {
+        let (ecTemp) = ec_add(quadruple_point,Prec.W15);
+        let (res) = ec_mulmuladd_W_inner(ecTemp, Prec, scalar_u, scalar_v, mm2);
+        return (res=res);
+    }
+
+    //shall not be reach
+    return (res=R);
+}

src/secp256r1/ec_mulmuladd_secp256r1.cairo

@@ -0,0 +1,74 @@
+//*************************************************************************************/
+///* Copyright (C) 2022 - Renaud Dubois - This file is part of Cairo_musig2 project   */
+///* License: This software is licensed under a dual BSD and GPL v2 license.          */
+///* See LICENSE file at the root folder of the project.                              */
+///* FILE: multipoint.cairo                                                           */
+///*                                                                                  */
+///*                                                                                  */
+///* DESCRIPTION: optimization of dual base multiplication                            */
+///* the algorithm combines the so called Shamir's trick with Windowing method        */
+//*************************************************************************************/
+
+// Shamir's trick:https://crypto.stackexchange.com/questions/99975/strauss-shamir-trick-on-ec-multiplication-by-scalar,
+// Windowing method : https://en.wikipedia.org/wiki/Exponentiation_by_squaring, section 'sliding window'
+// The implementation use a 2 bits window with trick, leading to a 16 points elliptic point precomputation
+
+from starkware.cairo.common.cairo_secp.bigint import BigInt3
+
+from src.secp256r1.ec import (
+    ec_add,
+    ec_double,
+    ec_mul,
+    EcPoint,
+)
+from src.secp256r1.ec_mulmuladd import (
+    Window,
+    ec_mulmuladd_W_inner,
+)
+
+
+func ec_mulmuladdW_bg3{range_check_ptr}(
+    G: EcPoint, Q: EcPoint,
+    scalar_u: BigInt3, scalar_v: BigInt3
+) -> (res: EcPoint){
+    alloc_locals;
+    local len_hi; //hi 84 bits part of scalar
+    local len_med; //med 86 bits part
+    local len_low;//low bits part
+
+    // Precompute a 4-bit window , W0=infty, W1=P, W2=Q,
+    // the window is indexed by (8*v1 4*u1+ 2*v0 + u0), where (u1,u0) represents two bit of scalar u,
+    // (resp for v)
+
+    let (W3) = ec_add(G,Q); //3:G+Q
+    let (W4) = ec_double(G); //4:2G
+    let (W5) = ec_add(G,W4); //5:3G
+    let (W6) = ec_add(W4,Q); //6:2G+Q
+    let (W7) = ec_add(W5,Q); //7:3G+Q
+    let (W8) = ec_double(Q); //8:2Q
+
+    let (W9) = ec_add(W8,G); //9:2Q+G
+    let (W10) = ec_add(W8,Q); //10:3Q
+    let (W11) = ec_add(W10,G); //11:3Q+G
+    let (W12) = ec_add(W8,W4); //12:2Q+2G
+    let (W13) = ec_add(W8,W5); //13:2Q+3G
+    let (W14) = ec_add(W10,W4); //14:3Q+2G
+    let (W15) = ec_add(W10,W5); //15:3Q+3G
+
+    let PrecPoint = Window( G,Q,W3, W4, W5, W6, W7, W8, W9, W10, W11, W12, W13, W14, W15);
+
+    //initialize R with infinity point
+    let R = EcPoint(BigInt3(0, 0, 0), BigInt3(0, 0, 0));
+
+    %{ ids.len_hi=max(ids.scalar_u.d2.bit_length(), ids.scalar_v.d2.bit_length())-1 %}
+    assert [range_check_ptr] = len_hi;
+    assert [range_check_ptr + 1] = 86 - len_hi;
+    let range_check_ptr = range_check_ptr + 2;
+
+    let (hiR )= ec_mulmuladd_W_inner(R, PrecPoint, scalar_u.d2, scalar_v.d2, len_hi);
+    let (medR) = ec_mulmuladd_W_inner(hiR, PrecPoint, scalar_u.d1, scalar_v.d1, 85);
+    let (lowR) = ec_mulmuladd_W_inner(medR, PrecPoint, scalar_u.d0, scalar_v.d0, 85);
+
+    return (res=lowR);
+
+}

src/secp256r1/signature.cairo

@@ -20,6 +20,7 @@ from src.secp256r1.constants import (
     GX0, GX1, GX2, GY0, GY1, GY2
 )
 from src.secp256r1.ec import ec_add, ec_mul
+from src.secp256r1.ec_mulmuladd_secp256r1 import ec_mulmuladdW_bg3
 from src.secp256r1.field import (
     unreduced_mul,
     unreduced_sqr,
@@ -133,10 +134,7 @@ func verify_secp256r1_signature{range_check_ptr}(
         let (u1: BigInt3) = div_mod_n(msg_hash, s);
         let (u2: BigInt3) = div_mod_n(r, s);
 
-        let (point1) = ec_mul(generator_point, u1);
-        let (point2) = ec_mul(public_key, u2);
-
-        let (point3) = ec_add(point1, point2);
+        let (point3) = ec_mulmuladdW_bg3(generator_point, public_key, u1, u2);
 
         let (x_mod_N) = div_mod_n(point3.x, BigInt3(d0=1, d1=0 ,d2=0));
         // We already validated r in [1, N) so no need to mod N it