e{12}.c: add GLV and GLS scalar multiplication methods.

supranational · Oct 25, 2020 · 03c3867 · 03c3867
1 parent e1b59c7
commit 03c3867
Show file tree

Hide file tree

Showing 4 changed files with 120 additions and 95 deletions.
diff --git a/src/e1.c b/src/e1.c
@@ -371,18 +371,63 @@ POINT_LADDER_POST_IMPL_A0(POINTonE1, 384, fp, onE1)
 POINT_MULT_SCALAR_LADDER_IMPL(POINTonE1)
 #endif
 
+static const vec384 beta = {            /* such that beta^3 - 1 = 0  */
+    /* -1/2 * (1 + sqrt(-3)) = ((P-2)^(P-2)) * (1 + (P-3)^((P+1)/4)) */
+    /* (0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4
+          897d29650fb85f9b409427eb4f49fffd8bfd00000000aaac << 384) % P */
+    TO_LIMB_T(0xcd03c9e48671f071), TO_LIMB_T(0x5dab22461fcda5d2),
+    TO_LIMB_T(0x587042afd3851b95), TO_LIMB_T(0x8eb60ebe01bacb9e),
+    TO_LIMB_T(0x03f97d6e83d050d2), TO_LIMB_T(0x18f0206554638741)
+};
+
+static void sigma(POINTonE1 *out, const POINTonE1 *in)
+{
+    vec_copy(out->X, in->X, 2*sizeof(out->X));
+    mul_fp(out->Z, in->Z, beta);
+}
+
+/* Gallant-Lambert-Vanstone, ~24-30% faster than POINTonE1_mult_w5 */
+static void POINTonE1_mult_glv(POINTonE1 *out, const POINTonE1 *in,
+                               const pow256 SK)
+{
+    union { vec256 l; pow256 s; } val;
+
+    /* SK/z^2 [in constant time] */
+
+    limbs_from_le_bytes(val.l, SK, 32);
+    div_by_zz(val.l);
+    le_bytes_from_limbs(val.s, val.l, 32);
+
+    {
+        const POINTonE1 *points[2];
+        const byte *scalars[2];
+        POINTonE1 in_sigma[1];
+
+        sigma(in_sigma, in);
+        POINTonE1_cneg(in_sigma, 1);
+        points[0] = in,         scalars[0] = val.s + 16;
+        points[1] = in_sigma,   scalars[1] = val.s;
+        POINTonE1s_mult_w4(out, points, 2, scalars, 128, NULL);
+        POINTonE1_cneg(out, 1);
+        mul_fp(out->Z, out->Z, beta);
+        mul_fp(out->Z, out->Z, beta);
+    }
+
+    vec_zero(val.l, sizeof(val));   /* scrub the copy of SK */
+}
+
 void blst_sk_to_pk_in_g1(POINTonE1 *out, const pow256 SK)
-{   POINTonE1_mult_w5(out, &BLS12_381_G1, SK, 255);   }
+{   POINTonE1_mult_glv(out, &BLS12_381_G1, SK);   }
 
 void blst_sign_pk_in_g2(POINTonE1 *out, const POINTonE1 *msg, const pow256 SK)
-{   POINTonE1_mult_w5(out, msg, SK, 255);   }
+{   POINTonE1_mult_glv(out, msg, SK);   }
 
 void blst_sk_to_pk2_in_g1(unsigned char out[96], POINTonE1_affine *PK,
                           const pow256 SK)
 {
     POINTonE1 P[1];
 
-    POINTonE1_mult_w5(P, &BLS12_381_G1, SK, 255);
+    POINTonE1_mult_glv(P, &BLS12_381_G1, SK);
     POINTonE1_from_Jacobian(P, P);
     if (PK != NULL)
         vec_copy(PK, P, sizeof(*PK));
@@ -398,7 +443,7 @@ void blst_sign_pk2_in_g2(unsigned char out[96], POINTonE1_affine *sig,
 {
     POINTonE1 P[1];
 
-    POINTonE1_mult_w5(P, hash, SK, 255);
+    POINTonE1_mult_glv(P, hash, SK);
     POINTonE1_from_Jacobian(P, P);
     if (sig != NULL)
         vec_copy(sig, P, sizeof(*sig));

diff --git a/src/e2.c b/src/e2.c
@@ -428,18 +428,83 @@ POINT_LADDER_POST_IMPL_A0(POINTonE2, 384x, fp2, onE2)
 POINT_MULT_SCALAR_LADDER_IMPL(POINTonE2)
 #endif
 
+static void psi(POINTonE2 *out, const POINTonE2 *in)
+{
+    static const vec384x frobenius_x = { /* 1/(1 + i)^((P-1)/3) */
+      { 0 },
+      { /* (0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4
+              897d29650fb85f9b409427eb4f49fffd8bfd00000000aaad << 384) % P */
+        TO_LIMB_T(0x890dc9e4867545c3), TO_LIMB_T(0x2af322533285a5d5),
+        TO_LIMB_T(0x50880866309b7e2c), TO_LIMB_T(0xa20d1b8c7e881024),
+        TO_LIMB_T(0x14e4f04fe2db9068), TO_LIMB_T(0x14e56d3f1564853a) }
+    };
+    static const vec384x frobenius_y = { /* 1/(1 + i)^((P-1)/2) */
+      { /* (0x135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60
+              ef396489f61eb45e304466cf3e67fa0af1ee7b04121bdea2 << 384) % P */
+        TO_LIMB_T(0x3e2f585da55c9ad1), TO_LIMB_T(0x4294213d86c18183),
+        TO_LIMB_T(0x382844c88b623732), TO_LIMB_T(0x92ad2afd19103e18),
+        TO_LIMB_T(0x1d794e4fac7cf0b9), TO_LIMB_T(0x0bd592fc7d825ec8) },
+      { /* (0x06af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e
+              77f76e17009241c5ee67992f72ec05f4c81084fbede3cc09 << 384) % P */
+        TO_LIMB_T(0x7bcfa7a25aa30fda), TO_LIMB_T(0xdc17dec12a927e7c),
+        TO_LIMB_T(0x2f088dd86b4ebef1), TO_LIMB_T(0xd1ca2087da74d4a7),
+        TO_LIMB_T(0x2da2596696cebc1d), TO_LIMB_T(0x0e2b7eedbbfd87d2) },
+    };
+
+    vec_copy(out, in, sizeof(*out));
+    cneg_fp(out->X[1], out->X[1], 1);   mul_fp2(out->X, out->X, frobenius_x);
+    cneg_fp(out->Y[1], out->Y[1], 1);   mul_fp2(out->Y, out->Y, frobenius_y);
+    cneg_fp(out->Z[1], out->Z[1], 1);
+}
+
+/* Galbraith-Lin-Scott, ~36-39% faster than POINTonE2_mul_w5 */
+static void POINTonE2_mult_gls(POINTonE2 *out, const POINTonE2 *in,
+                               const pow256 SK)
+{
+    union { vec256 l; pow256 s; } val;
+
+    /* break down SK to "digits" with |z| as radix [in constant time] */
+
+    limbs_from_le_bytes(val.l, SK, 32);
+    div_by_zz(val.l);
+    div_by_z(val.l);
+    div_by_z(val.l + NLIMBS(256)/2);
+    le_bytes_from_limbs(val.s, val.l, 32);
+
+    {
+        const POINTonE2 *points[4];
+        const byte *scalars[4];
+        POINTonE2 P[3];
+
+        psi(&P[0], in);
+        psi(&P[1], &P[0]);
+        psi(&P[2], &P[1]);
+
+        POINTonE2_cneg(&P[0], 1);   /* account for z being negative */
+        POINTonE2_cneg(&P[2], 1);
+
+        points[0] = in,     scalars[0] = val.s;
+        points[1] = &P[0],  scalars[1] = val.s + 8;
+        points[2] = &P[1],  scalars[2] = val.s + 16;
+        points[3] = &P[2],  scalars[3] = val.s + 24;
+        POINTonE2s_mult_w4(out, points, 4, scalars, 64, NULL);
+    }
+
+    vec_zero(val.l, sizeof(val));   /* scrub the copy of SK */
+}
+
 void blst_sk_to_pk_in_g2(POINTonE2 *out, const pow256 SK)
-{   POINTonE2_mult_w5(out, &BLS12_381_G2, SK, 255);   }
+{   POINTonE2_mult_gls(out, &BLS12_381_G2, SK);   }
 
 void blst_sign_pk_in_g1(POINTonE2 *out, const POINTonE2 *msg, const pow256 SK)
-{   POINTonE2_mult_w5(out, msg, SK, 255);   }
+{   POINTonE2_mult_gls(out, msg, SK);   }
 
 void blst_sk_to_pk2_in_g2(unsigned char out[192], POINTonE2_affine *PK,
                           const pow256 SK)
 {
     POINTonE2 P[1];
 
-    POINTonE2_mult_w5(P, &BLS12_381_G2, SK, 255);
+    POINTonE2_mult_gls(P, &BLS12_381_G2, SK);
     POINTonE2_from_Jacobian(P, P);
     if (PK != NULL)
         vec_copy(PK, P, sizeof(*PK));
@@ -455,7 +520,7 @@ void blst_sign_pk2_in_g1(unsigned char out[192], POINTonE2_affine *sig,
 {
     POINTonE2 P[1];
 
-    POINTonE2_mult_w5(P, hash, SK, 255);
+    POINTonE2_mult_gls(P, hash, SK);
     POINTonE2_from_Jacobian(P, P);
     if (sig != NULL)
         vec_copy(sig, P, sizeof(*sig));

diff --git a/src/map_to_g1.c b/src/map_to_g1.c
@@ -494,23 +494,9 @@ static void POINTonE1_times_zz_minus_1_div_by_3(POINTonE1 *out,
 }
 #endif
 
-static void sigma(POINTonE1 *out, const POINTonE1 *in)
-{
-    static const vec384 beta = {            /* such that beta^3 - 1 = 0  */
-        /* -1/2 * (1 + sqrt(-3)) = ((P-2)^(P-2)) * (1 + (P-3)^((P+1)/4)) */
-        /* (0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4
-              897d29650fb85f9b409427eb4f49fffd8bfd00000000aaac << 384) % P */
-        TO_LIMB_T(0xcd03c9e48671f071), TO_LIMB_T(0x5dab22461fcda5d2),
-        TO_LIMB_T(0x587042afd3851b95), TO_LIMB_T(0x8eb60ebe01bacb9e),
-        TO_LIMB_T(0x03f97d6e83d050d2), TO_LIMB_T(0x18f0206554638741)
-    };
-
-    mul_fp(out->X, in->X, beta);
-    vec_copy(out->Y, in->Y, 2*sizeof(out->Y));
-}
-
 static limb_t POINTonE1_in_G1(const POINTonE1_affine *p)
 {
+    void sigma(POINTonE1 *out, const POINTonE1 *in);
     POINTonE1 t0, t1, t2;
 
     vec_copy(t2.X, p->X, 2*sizeof(t0.X));

diff --git a/src/map_to_g2.c b/src/map_to_g2.c
@@ -307,78 +307,6 @@ static void clear_cofactor(POINTonE2 *out, const POINTonE2 *p)
  * As per suggestions in "7. Clearing the cofactor" at
  * https://tools.ietf.org/html/draft-irtf-cfrg-hash-to-curve-06
  */
-static const vec384x iwsc = {           /* 1/(1 + i) = 1/2 - 1/2*i      */
-  { TO_LIMB_T(0x1804000000015554), TO_LIMB_T(0x855000053ab00001),
-    TO_LIMB_T(0x633cb57c253c276f), TO_LIMB_T(0x6e22d1ec31ebb502),
-    TO_LIMB_T(0xd3916126f2d14ca2), TO_LIMB_T(0x17fbb8571a006596) },
-  { TO_LIMB_T(0xa1fafffffffe5557), TO_LIMB_T(0x995bfff976a3fffe),
-    TO_LIMB_T(0x03f41d24d174ceb4), TO_LIMB_T(0xf6547998c1995dbd),
-    TO_LIMB_T(0x778a468f507a6034), TO_LIMB_T(0x020559931f7f8103) }
-};
-
-static const vec384x k_cx = {           /* iwsc^((P-1)/3)               */
-  { 0 },
-  { /* (0x1a0111ea397fe699ec02408663d4de85aa0d857d89759ad4
-          897d29650fb85f9b409427eb4f49fffd8bfd00000000aaad << 384) % P */
-    TO_LIMB_T(0x890dc9e4867545c3), TO_LIMB_T(0x2af322533285a5d5),
-    TO_LIMB_T(0x50880866309b7e2c), TO_LIMB_T(0xa20d1b8c7e881024),
-    TO_LIMB_T(0x14e4f04fe2db9068), TO_LIMB_T(0x14e56d3f1564853a) }
-};
-
-static const vec384x k_cy = {           /* iwsc^((P-1)/2)               */
-  { /* (0x135203e60180a68ee2e9c448d77a2cd91c3dedd930b1cf60
-          ef396489f61eb45e304466cf3e67fa0af1ee7b04121bdea2 << 384) % P */
-    TO_LIMB_T(0x3e2f585da55c9ad1), TO_LIMB_T(0x4294213d86c18183),
-    TO_LIMB_T(0x382844c88b623732), TO_LIMB_T(0x92ad2afd19103e18),
-    TO_LIMB_T(0x1d794e4fac7cf0b9), TO_LIMB_T(0x0bd592fc7d825ec8) },
-  { /* (0x06af0e0437ff400b6831e36d6bd17ffe48395dabc2d3435e
-          77f76e17009241c5ee67992f72ec05f4c81084fbede3cc09 << 384) % P */
-    TO_LIMB_T(0x7bcfa7a25aa30fda), TO_LIMB_T(0xdc17dec12a927e7c),
-    TO_LIMB_T(0x2f088dd86b4ebef1), TO_LIMB_T(0xd1ca2087da74d4a7),
-    TO_LIMB_T(0x2da2596696cebc1d), TO_LIMB_T(0x0e2b7eedbbfd87d2) },
-};
-
-static void qi_x_iwsc(vec384x out, const vec384x in)
-{
-    mul_fp2(out, in, iwsc);
-    mul_fp(out[0], out[0], k_cx[1]);
-    mul_fp(out[1], out[1], k_cx[1]);
-    neg_fp(out[1], out[1]);
-}
-
-static void qi_y_iwsc(vec384x out, const vec384x in)
-{
-    vec384x t;
-
-    mul_fp2(t, in, iwsc);
-    add_fp(out[0], t[0], t[1]);
-    sub_fp(out[1], t[0], t[1]);
-    mul_fp(out[0], out[0], k_cy[1]);
-    mul_fp(out[1], out[1], k_cy[1]);
-}
-
-static void psi(POINTonE2 *out, const POINTonE2 *in)
-{
-    vec384x xn, xd, yn, yd;
-
-    sqr_fp2(xd, in->Z);                 /* xd = Z^2                     */
-    mul_fp2(yd, xd, in->Z);             /* yd = Z^3                     */
-
-    qi_x_iwsc(xn, in->X);   mul_fp2(xn, xn, k_cx);
-    qi_x_iwsc(xd, xd);
-
-    qi_y_iwsc(yn, in->Y);   mul_fp2(yn, yn, k_cy);
-    qi_y_iwsc(yd, yd);
-
-    /* convert (xn, xd, yn, yd) to Jacobian coordinates                 */
-    mul_fp2(out->Z, xd, yd);            /* Z = xd * yd                  */
-    mul_fp2(out->X, xn, yd);
-    mul_fp2(out->X, out->X, out->Z);    /* X = xn * xd * yd^2           */
-    sqr_fp2(out->Y, out->Z);
-    mul_fp2(out->Y, out->Y, xd);
-    mul_fp2(out->Y, out->Y, yn);        /* Y = yn * xd^3 * yd^2         */
-}
-
 static void POINTonE2_add_n_dbl(POINTonE2 *out, const POINTonE2 *p, size_t n)
 {
     POINTonE2_add(out, out, p);
@@ -398,6 +326,7 @@ static void POINTonE2_times_minus_z(POINTonE2 *out, const POINTonE2 *in)
 
 static void clear_cofactor(POINTonE2 *out, const POINTonE2 *p)
 {
+    void psi(POINTonE2 *out, const POINTonE2 *in);
     POINTonE2 t0, t1;
 
     /* A.Budroni, F.Pintore, "Efficient hash maps to G2 on BLS curves"  */