Use s2n-bignum P-521 scalar multiplication and modular inverse

This replaces the general (fresh, not precomputed, point) scalar
multiplication with the corresponding function p521_jscalarmul or
p521_jscalarmul_alt from s2n-bignum, and also replaces the Fermat
inverse in p521.c with the markedly faster and formally verified
divstep-based code from s2n-bignum, bignum_inv_p521.

crypto/fipsmodule/CMakeLists.txt

@@ -226,6 +226,11 @@ if((((ARCH STREQUAL "x86_64") AND NOT MY_ASSEMBLER_IS_TOO_OLD_FOR_512AVX) OR
     p521/bignum_sqr_p521_alt.S
     p521/bignum_tolebytes_p521.S
     p521/bignum_fromlebytes_p521.S
+    p521/bignum_inv_p521.S
+    p521/p521_jdouble.S
+    p521/p521_jdouble_alt.S
+    p521/p521_jscalarmul.S
+    p521/p521_jscalarmul_alt.S
 
     curve25519/bignum_mod_n25519.S
     curve25519/bignum_neg_p25519.S

crypto/fipsmodule/ec/p521.c

@@ -196,6 +196,14 @@ static void p521_to_generic(EC_FELEM *out, const p521_felem in) {
 #endif
 }
 
+#if defined(EC_NISTP_USE_S2N_BIGNUM)
+
+static void p521_felem_inv(p521_felem output, const p521_felem t1) {
+  bignum_inv_p521(output,t1);
+}
+
+#else
+
 // Finite field inversion using Fermat Little Theorem.
 // The code is autogenerated by the ECCKiila project:
 //   https://arxiv.org/abs/2007.11481
@@ -259,6 +267,32 @@ static void p521_felem_inv(p521_felem output, const p521_felem t1) {
     p521_felem_mul(output, acc, t1);
 }
 
+#endif
+
+#if defined(EC_NISTP_USE_S2N_BIGNUM)
+
+static void p521_point_double(p521_felem x_out,
+                              p521_felem y_out,
+                              p521_felem z_out,
+                              const p521_felem x_in,
+                              const p521_felem y_in,
+                              const p521_felem z_in) {
+
+  uint64_t s2n_point[27], s2n_result[27];
+
+  OPENSSL_memcpy(s2n_point,x_in,72);
+  OPENSSL_memcpy(s2n_point+9,y_in,72);
+  OPENSSL_memcpy(s2n_point+18,z_in,72);
+
+  p521_jdouble_selector(s2n_result,s2n_point);
+
+  OPENSSL_memcpy(x_out,s2n_result,72);
+  OPENSSL_memcpy(y_out,s2n_result+9,72);
+  OPENSSL_memcpy(z_out,s2n_result+18,72);
+}
+
+#else
+
 static void p521_point_double(p521_felem x_out,
                               p521_felem y_out,
                               p521_felem z_out,
@@ -268,6 +302,8 @@ static void p521_point_double(p521_felem x_out,
   ec_nistp_point_double(p521_methods(), x_out, y_out, z_out, x_in, y_in, z_in);
 }
 
+#endif
+
 // p521_point_add calculates (x1, y1, z1) + (x2, y2, z2)
 //
 // The method is taken from:
@@ -446,6 +482,28 @@ static void p521_select_point_affine(p521_felem out[2],
 }
 
 // Multiplication of an arbitrary point by a scalar, r = [scalar]P.
+
+#if defined(EC_NISTP_USE_S2N_BIGNUM)
+
+static void ec_GFp_nistp521_point_mul(const EC_GROUP *group, EC_JACOBIAN *r,
+                                      const EC_JACOBIAN *p,
+                                      const EC_SCALAR *scalar) {
+
+  uint64_t s2n_point[27], s2n_result[27];
+
+  OPENSSL_memcpy(s2n_point,p->X.words,72);
+  OPENSSL_memcpy(s2n_point+9,p->Y.words,72);
+  OPENSSL_memcpy(s2n_point+18,p->Z.words,72);
+
+  p521_jscalarmul_selector(s2n_result,(uint64_t*)scalar,s2n_point);
+
+  OPENSSL_memcpy(r->X.words,s2n_result,72);
+  OPENSSL_memcpy(r->Y.words,s2n_result+9,72);
+  OPENSSL_memcpy(r->Z.words,s2n_result+18,72);
+}
+
+#else
+
 static void ec_GFp_nistp521_point_mul(const EC_GROUP *group, EC_JACOBIAN *r,
                                       const EC_JACOBIAN *p,
                                       const EC_SCALAR *scalar) {
@@ -463,6 +521,8 @@ static void ec_GFp_nistp521_point_mul(const EC_GROUP *group, EC_JACOBIAN *r,
   p521_to_generic(&r->Z, res[2]);
 }
 
+#endif
+
 // Include the precomputed table for the based point scalar multiplication.
 #include "p521_table.h"