ECDH from compressed points with "free" sqrt

.gitignore

@@ -1,5 +1,6 @@
 bench_inv
 bench_ecdh
+bench_ecdh_opt
 bench_sign
 bench_verify
 bench_schnorr_verify

include/secp256k1_ecdh.h

@@ -24,6 +24,14 @@ SECP256K1_API SECP256K1_WARN_UNUSED_RESULT int secp256k1_ecdh(
   const unsigned char *privkey
 ) SECP256K1_ARG_NONNULL(1) SECP256K1_ARG_NONNULL(2) SECP256K1_ARG_NONNULL(3) SECP256K1_ARG_NONNULL(4);
 
+SECP256K1_API SECP256K1_WARN_UNUSED_RESULT int secp256k1_ecdh_opt(
+  const secp256k1_context* ctx,
+  unsigned char *result,
+  const unsigned char *pub,
+  size_t publen,
+  const unsigned char *privkey
+) SECP256K1_ARG_NONNULL(1) SECP256K1_ARG_NONNULL(2) SECP256K1_ARG_NONNULL(3) SECP256K1_ARG_NONNULL(5);
+
 # ifdef __cplusplus
 }
 # endif

src/bench_ecdh.c

@@ -13,13 +13,13 @@
 
 typedef struct {
     secp256k1_context *ctx;
-    secp256k1_pubkey point;
+    secp256k1_pubkey pubkey;
     unsigned char scalar[32];
-} bench_ecdh_t;
+} bench_ecdh_data;
 
 static void bench_ecdh_setup(void* arg) {
     int i;
-    bench_ecdh_t *data = (bench_ecdh_t*)arg;
+    bench_ecdh_data *data = (bench_ecdh_data*)arg;
     const unsigned char point[] = {
         0x03,
         0x54, 0x94, 0xc1, 0x5d, 0x32, 0x09, 0x97, 0x06,
@@ -28,25 +28,25 @@ static void bench_ecdh_setup(void* arg) {
         0xa2, 0xba, 0xd1, 0x84, 0xf8, 0x83, 0xc6, 0x9f
     };
 
-    data->ctx = secp256k1_context_create(0);
+    data->ctx = secp256k1_context_create(SECP256K1_CONTEXT_NONE);
     for (i = 0; i < 32; i++) {
         data->scalar[i] = i + 1;
     }
-    CHECK(secp256k1_ec_pubkey_parse(data->ctx, &data->point, point, sizeof(point)) == 1);
+    CHECK(secp256k1_ec_pubkey_parse(data->ctx, &data->pubkey, point, sizeof(point)) == 1);
 }
 
 static void bench_ecdh(void* arg) {
     int i;
     unsigned char res[32];
-    bench_ecdh_t *data = (bench_ecdh_t*)arg;
+    bench_ecdh_data *data = (bench_ecdh_data*)arg;
 
     for (i = 0; i < 20000; i++) {
-        CHECK(secp256k1_ecdh(data->ctx, res, &data->point, data->scalar) == 1);
+        CHECK(secp256k1_ecdh(data->ctx, res, &data->pubkey, data->scalar) == 1);
     }
 }
 
 int main(void) {
-    bench_ecdh_t data;
+    bench_ecdh_data data;
 
     run_benchmark("ecdh", bench_ecdh, bench_ecdh_setup, NULL, &data, 10, 20000);
     return 0;

src/bench_ecdh_opt.c

@@ -0,0 +1,54 @@
+/**********************************************************************
+ * Copyright (c) 2015 Pieter Wuille, Andrew Poelstra                  *
+ * Distributed under the MIT software license, see the accompanying   *
+ * file COPYING or http://www.opensource.org/licenses/mit-license.php.*
+ **********************************************************************/
+
+#include <string.h>
+
+#include "include/secp256k1.h"
+#include "include/secp256k1_ecdh.h"
+#include "util.h"
+#include "bench.h"
+
+typedef struct {
+    secp256k1_context *ctx;
+    unsigned char point[33];
+    unsigned char scalar[32];
+} bench_ecdh_data;
+
+static void bench_ecdh_setup(void* arg) {
+    int i;
+    bench_ecdh_data *data = (bench_ecdh_data*)arg;
+    const unsigned char point[] = {
+        0x03,
+        0x54, 0x94, 0xc1, 0x5d, 0x32, 0x09, 0x97, 0x06,
+        0xc2, 0x39, 0x5f, 0x94, 0x34, 0x87, 0x45, 0xfd,
+        0x75, 0x7c, 0xe3, 0x0e, 0x4e, 0x8c, 0x90, 0xfb,
+        0xa2, 0xba, 0xd1, 0x84, 0xf8, 0x83, 0xc6, 0x9f
+    };
+
+    data->ctx = secp256k1_context_create(SECP256K1_CONTEXT_VERIFY);
+    for (i = 0; i < 32; i++) {
+        data->scalar[i] = i + 1;
+    }
+    CHECK(sizeof(point) == sizeof(data->point));
+    memcpy(data->point, point, sizeof(point));
+}
+
+static void bench_ecdh(void* arg) {
+    int i;
+    unsigned char res[32];
+    bench_ecdh_data *data = (bench_ecdh_data*)arg;
+
+    for (i = 0; i < 20000; i++) {
+        CHECK(secp256k1_ecdh_opt(data->ctx, res, data->point, sizeof(data->point), data->scalar) == 1);
+    }
+}
+
+int main(void) {
+    bench_ecdh_data data;
+
+    run_benchmark("ecdh", bench_ecdh, bench_ecdh_setup, NULL, &data, 10, 20000);
+    return 0;
+}

src/bench_internal.c

@@ -191,6 +191,26 @@ void bench_field_sqrt_var(void* arg) {
     }
 }
 
+void bench_field_rsqrt_var(void* arg) {
+    int i;
+    bench_inv_t *data = (bench_inv_t*)arg;
+
+    for (i = 0; i < 20000; i++) {
+        secp256k1_fe_rsqrt_var(&data->fe_x, &data->fe_y, &data->fe_x);
+        secp256k1_fe_add(&data->fe_x, &data->fe_y);
+    }
+}
+
+void bench_field_par_rsqrt_inv_var(void* arg) {
+    int i;
+    bench_inv_t *data = (bench_inv_t*)arg;
+
+    for (i = 0; i < 20000; i++) {
+        secp256k1_fe_par_rsqrt_inv_var(&data->fe_x, &data->fe_y, &data->fe_x, &data->fe_y);
+        secp256k1_fe_add(&data->fe_x, &data->fe_y);
+    }
+}
+
 void bench_group_double_var(void* arg) {
     int i;
     bench_inv_t *data = (bench_inv_t*)arg;
@@ -334,6 +354,8 @@ int main(int argc, char **argv) {
     if (have_flag(argc, argv, "field") || have_flag(argc, argv, "inverse")) run_benchmark("field_inverse", bench_field_inverse, bench_setup, NULL, &data, 10, 20000);
     if (have_flag(argc, argv, "field") || have_flag(argc, argv, "inverse")) run_benchmark("field_inverse_var", bench_field_inverse_var, bench_setup, NULL, &data, 10, 20000);
     if (have_flag(argc, argv, "field") || have_flag(argc, argv, "sqrt")) run_benchmark("field_sqrt_var", bench_field_sqrt_var, bench_setup, NULL, &data, 10, 20000);
+    if (have_flag(argc, argv, "field") || have_flag(argc, argv, "sqrt")) run_benchmark("field_rsqrt_var", bench_field_rsqrt_var, bench_setup, NULL, &data, 10, 20000);
+    if (have_flag(argc, argv, "field") || have_flag(argc, argv, "sqrt") || have_flag(argc, argv, "inverse")) run_benchmark("field_par_rsqrt_inv_var", bench_field_par_rsqrt_inv_var, bench_setup, NULL, &data, 10, 20000);
 
     if (have_flag(argc, argv, "group") || have_flag(argc, argv, "double")) run_benchmark("group_double_var", bench_group_double_var, bench_setup, NULL, &data, 10, 200000);
     if (have_flag(argc, argv, "group") || have_flag(argc, argv, "add")) run_benchmark("group_add_var", bench_group_add_var, bench_setup, NULL, &data, 10, 200000);

src/field.h

@@ -39,13 +39,11 @@ static void secp256k1_fe_normalize_weak(secp256k1_fe *r);
 /** Normalize a field element, without constant-time guarantee. */
 static void secp256k1_fe_normalize_var(secp256k1_fe *r);
 
-/** Verify whether a field element represents zero i.e. would normalize to a zero value. The field
- *  implementation may optionally normalize the input, but this should not be relied upon. */
-static int secp256k1_fe_normalizes_to_zero(secp256k1_fe *r);
+/** Verify whether a field element represents zero i.e. would normalize to a zero value. */
+static int secp256k1_fe_normalizes_to_zero(const secp256k1_fe *r);
 
-/** Verify whether a field element represents zero i.e. would normalize to a zero value. The field
- *  implementation may optionally normalize the input, but this should not be relied upon. */
-static int secp256k1_fe_normalizes_to_zero_var(secp256k1_fe *r);
+/** Verify whether a field element represents zero i.e. would normalize to a zero value. */
+static int secp256k1_fe_normalizes_to_zero_var(const secp256k1_fe *r);
 
 /** Set a field element equal to a small integer. Resulting field element is normalized. */
 static void secp256k1_fe_set_int(secp256k1_fe *r, int a);
@@ -88,12 +86,28 @@ static void secp256k1_fe_mul(secp256k1_fe *r, const secp256k1_fe *a, const secp2
 static void secp256k1_fe_sqr(secp256k1_fe *r, const secp256k1_fe *a);
 
 /** If a has a square root, it is computed in r and 1 is returned. If a does not
- *  have a square root, the root of its negation is computed and 0 is returned.
+ *  have a square root, the root of -a is computed and 0 is returned.
  *  The input's magnitude can be at most 8. The output magnitude is 1 (but not
  *  guaranteed to be normalized). The result in r will always be a square
  *  itself. */
 static int secp256k1_fe_sqrt_var(secp256k1_fe *r, const secp256k1_fe *a);
 
+/** If a has a square root, the square root is computed in rs, its reciprocal square root is
+ *  calculated in rr, and 1 is returned. If a does not have a square root, the root (and recip. root)
+ *  of -a are computed and 0 is returned. The input's magnitude can be at most 8. The
+ *  outputs' magnitudes are 1 (but not guaranteed to be normalized). The result in rs will always be
+ *  a square itself. The result in rr will be a square if, and only if, a is a square.
+ */
+static int secp256k1_fe_rsqrt_var(secp256k1_fe *rs, secp256k1_fe *rr, const secp256k1_fe *a);
+
+/** Parallel reciprocal square root and inverse. Sets ri to be the (modular) inverse of b. If a has a
+ *  square root, the reciprocal of its square root is computed in rr and 1 is returned. If a does not
+ *  have a square root, the reciprocal root of -a is computed and 0 is returned. The inputs'
+ *  magnitudes can be at most 8. The outputs' magnitudes are 1 (but not guaranteed to be normalized).
+ *  The result in rr will be a square if, and only if, a is a square.
+ */
+static int secp256k1_fe_par_rsqrt_inv_var(secp256k1_fe *rr,  secp256k1_fe *ri, const secp256k1_fe *a, const secp256k1_fe *b);
+
 /** Sets a field element to be the (modular) inverse of another. Requires the input's magnitude to be
  *  at most 8. The output magnitude is 1 (but not guaranteed to be normalized). */
 static void secp256k1_fe_inv(secp256k1_fe *r, const secp256k1_fe *a);

src/field_10x26_impl.h

@@ -188,7 +188,7 @@ static void secp256k1_fe_normalize_var(secp256k1_fe *r) {
 #endif
 }
 
-static int secp256k1_fe_normalizes_to_zero(secp256k1_fe *r) {
+static int secp256k1_fe_normalizes_to_zero(const secp256k1_fe *r) {
     uint32_t t0 = r->n[0], t1 = r->n[1], t2 = r->n[2], t3 = r->n[3], t4 = r->n[4],
              t5 = r->n[5], t6 = r->n[6], t7 = r->n[7], t8 = r->n[8], t9 = r->n[9];
 
@@ -217,7 +217,7 @@ static int secp256k1_fe_normalizes_to_zero(secp256k1_fe *r) {
     return (z0 == 0) | (z1 == 0x3FFFFFFUL);
 }
 
-static int secp256k1_fe_normalizes_to_zero_var(secp256k1_fe *r) {
+static int secp256k1_fe_normalizes_to_zero_var(const secp256k1_fe *r) {
     uint32_t t0, t1, t2, t3, t4, t5, t6, t7, t8, t9;
     uint32_t z0, z1;
     uint32_t x;

src/field_5x52_impl.h

@@ -167,7 +167,7 @@ static void secp256k1_fe_normalize_var(secp256k1_fe *r) {
 #endif
 }
 
-static int secp256k1_fe_normalizes_to_zero(secp256k1_fe *r) {
+static int secp256k1_fe_normalizes_to_zero(const secp256k1_fe *r) {
     uint64_t t0 = r->n[0], t1 = r->n[1], t2 = r->n[2], t3 = r->n[3], t4 = r->n[4];
 
     /* z0 tracks a possible raw value of 0, z1 tracks a possible raw value of P */
@@ -190,7 +190,7 @@ static int secp256k1_fe_normalizes_to_zero(secp256k1_fe *r) {
     return (z0 == 0) | (z1 == 0xFFFFFFFFFFFFFULL);
 }
 
-static int secp256k1_fe_normalizes_to_zero_var(secp256k1_fe *r) {
+static int secp256k1_fe_normalizes_to_zero_var(const secp256k1_fe *r) {
     uint64_t t0, t1, t2, t3, t4;
     uint64_t z0, z1;
     uint64_t x;