Add Jacobi symbol test via GMP

Also add native Jacobi symbol test (Andrew)

Rebased-by: Andrew Poelstra

src/bench_internal.c

@@ -299,6 +299,19 @@ void bench_context_sign(void* arg) {
     }
 }
 
+void bench_num_jacobi(void* arg) {
+    int i;
+    bench_inv_t *data = (bench_inv_t*)arg;
+    secp256k1_num nx, norder;
+
+    secp256k1_scalar_get_num(&nx, &data->scalar_x);
+    secp256k1_scalar_order_get_num(&norder);
+    secp256k1_scalar_get_num(&norder, &data->scalar_y);
+
+    for (i = 0; i < 200000; i++) {
+        secp256k1_num_jacobi(&nx, &norder);
+    }
+}
 
 int have_flag(int argc, char** argv, char *flag) {
     char** argm = argv + argc;
@@ -350,5 +363,6 @@ int main(int argc, char **argv) {
     if (have_flag(argc, argv, "context") || have_flag(argc, argv, "verify")) run_benchmark("context_verify", bench_context_verify, bench_setup, NULL, &data, 10, 20);
     if (have_flag(argc, argv, "context") || have_flag(argc, argv, "sign")) run_benchmark("context_sign", bench_context_sign, bench_setup, NULL, &data, 10, 200);
 
+    if (have_flag(argc, argv, "num") || have_flag(argc, argv, "jacobi")) run_benchmark("num_jacobi", bench_num_jacobi, bench_setup, NULL, &data, 10, 200000);
     return 0;
 }

src/num.h

@@ -34,6 +34,9 @@ static void secp256k1_num_set_bin(secp256k1_num *r, const unsigned char *a, unsi
 /** Compute a modular inverse. The input must be less than the modulus. */
 static void secp256k1_num_mod_inverse(secp256k1_num *r, const secp256k1_num *a, const secp256k1_num *m);
 
+/** Compute the jacobi symbol (a|b). b must be positive and odd. */
+static int secp256k1_num_jacobi(const secp256k1_num *a, const secp256k1_num *b);
+
 /** Compare the absolute value of two numbers. */
 static int secp256k1_num_cmp(const secp256k1_num *a, const secp256k1_num *b);
 

src/num_5x64.h

@@ -12,6 +12,7 @@
 #define NUM_N_WORDS	5
 #define NUM_WORD_WIDTH	64
 #define NUM_WORD_CTLZ __builtin_clzl
+#define NUM_WORD_CTZ __builtin_ctzl
 typedef uint64_t secp256k1_num_word;
 typedef int64_t secp256k1_num_sword;
 typedef uint128_t secp256k1_num_dword;

src/num_9x32.h

@@ -12,6 +12,7 @@
 #define NUM_N_WORDS	9
 #define NUM_WORD_WIDTH	32
 #define NUM_WORD_CTLZ __builtin_clz
+#define NUM_WORD_CTZ __builtin_ctz
 typedef uint32_t secp256k1_num_word;
 typedef int32_t secp256k1_num_sword;
 typedef uint64_t secp256k1_num_dword;

src/num_gmp_impl.h

@@ -144,6 +144,28 @@ static void secp256k1_num_mod_inverse(secp256k1_num *r, const secp256k1_num *a,
     memset(v, 0, sizeof(v));
 }
 
+static int secp256k1_num_jacobi(const secp256k1_num *a, const secp256k1_num *b) {
+    int ret;
+    mpz_t ga, gb;
+    secp256k1_num_sanity(a);
+    secp256k1_num_sanity(b);
+    VERIFY_CHECK(!b->neg && (b->limbs > 0) && (b->data[0] & 1));
+
+    mpz_inits(ga, gb, NULL);
+
+    mpz_import(gb, b->limbs, -1, sizeof(mp_limb_t), 0, 0, b->data);
+    mpz_import(ga, a->limbs, -1, sizeof(mp_limb_t), 0, 0, a->data);
+    if (a->neg) {
+        mpz_neg(ga, ga);
+    }
+
+    ret = mpz_jacobi(ga, gb);
+
+    mpz_clears(ga, gb, NULL);
+
+    return ret;
+}
+
 static int secp256k1_num_is_zero(const secp256k1_num *a) {
     return (a->limbs == 1 && a->data[0] == 0);
 }

src/num_native_impl.h

@@ -96,6 +96,16 @@ static void secp256k1_num_shift(secp256k1_num *r, int bits) {
         r->data[i] = 0;
 }
 
+SECP256K1_INLINE static int secp256k1_num_is_one(const secp256k1_num *a) {
+    int i;
+    if (a->data[0] != 1)
+        return 0;
+    for (i = 1; i < NUM_N_WORDS - 1; ++i)
+        if (a->data[i] != 0)
+            return 0;
+    return 1;
+}
+
 SECP256K1_INLINE static int secp256k1_num_is_zero(const secp256k1_num *a) {
     int i;
     for (i = 0; i < NUM_N_WORDS - 1; ++i)
@@ -591,4 +601,80 @@ static void secp256k1_num_mod_inverse(secp256k1_num *rr, const secp256k1_num *a,
 }
 /* end mod inverse */
 
+/* start jacobi symbol */
+/* Compute a number modulo some power of 2 */
+SECP256K1_INLINE static int secp256k1_num_mod_2(const secp256k1_num *a, int m) {
+    VERIFY_CHECK(m > 0);
+    VERIFY_CHECK((m & (m - 1)) == 0);  /* check that m is a power of 2 */
+    /* Since our words are powers of 2 we only need to mod the lowest digit */
+    return a->data[0] % m;
+}
+
+static int secp256k1_num_jacobi_1(secp256k1_num_word a, secp256k1_num_word b) {
+    int ret = 1;
+    secp256k1_num_word t;
+    /* Iterate, left-multiplying it by [[0 1] [1 -w]] as many times as we can. */
+    while (1) {
+        a %= b;
+        if (a == 0)
+            return 0;
+        if (a % 2 == 0) {
+            int shift = NUM_WORD_CTZ(a);
+            a >>= shift;
+            if ((b % 8 == 3 || b % 8 == 5) && shift % 2 == 1)
+                ret *= -1;
+        }
+        if (a == 1)
+            break;
+        if (b % 4 == 3 && a % 4 == 3)
+            ret *= -1;
+        t = a; a = b; b = t;
+    }
+    return ret;
+}
+
+/* Compute the Jacobian symbol (a|b) assuming b is an odd prime */
+static int secp256k1_num_jacobi(const secp256k1_num *a, const secp256k1_num *b) {
+    secp256k1_num top = *a, bot = *b, scratch;
+    secp256k1_num_word x, y;
+    int index[2];
+    int ret = 1;
+
+    while (1) {
+        int mod8 = secp256k1_num_mod_2(&bot, 8);
+        secp256k1_num_leading_digit(&x, &index[0], &top);
+        secp256k1_num_leading_digit(&y, &index[1], &bot);
+
+        if (index[0] == 0 && index[1] == 0)
+            return ret * secp256k1_num_jacobi_1(x, y);
+
+        /* Algorithm from https://en.wikipedia.org/wiki/Jacobi_symbol#Calculating_the_Jacobi_symbol */
+        secp256k1_num_div_mod(&scratch, &top, &top, &bot); /* top <- top mod bottom */
+
+        /* are we done? */
+        if (secp256k1_num_is_zero(&top))
+            return 0;
+
+        /* cast out powers of two from the "numerator" */
+        while (secp256k1_num_mod_2(&top, 2) == 0) {
+            int shift = NUM_WORD_CTZ(top.data[0]);
+            secp256k1_num_shift(&top, shift);
+            if ((mod8 == 3 || mod8 == 5) && shift % 2 == 1)
+                ret *= -1;
+        }
+
+        /* are we done? */
+        if (secp256k1_num_is_one(&top))
+            return ret;
+        /* if not, iterate */
+        if (mod8 % 4 == 3 && secp256k1_num_mod_2(&top, 4) == 3)
+            ret *= -1;
+
+        scratch = top;
+        top = bot;
+        bot = scratch;
+    }
+}
+/* end jacobi symbol */
+
 #endif

src/tests.c

@@ -498,11 +498,32 @@ void test_num_add_sub(void) {
     CHECK(secp256k1_num_eq(&n2p1, &n1));
 }
 
+void test_num_jacobi(void) {
+    secp256k1_scalar sqr;
+    secp256k1_scalar five;  /* five is not a quadratic residue */
+    secp256k1_num order, n;
+
+    /* setup values */
+    random_scalar_order_test(&sqr);
+    secp256k1_scalar_sqr(&sqr, &sqr);
+    secp256k1_scalar_set_int(&five, 5);
+    secp256k1_scalar_order_get_num(&order);
+
+    /* test residue */
+    secp256k1_scalar_get_num(&n, &sqr);
+    CHECK(secp256k1_num_jacobi(&n, &order) == 1);
+    /* test nonresidue */
+    secp256k1_scalar_mul(&sqr, &sqr, &five);
+    secp256k1_scalar_get_num(&n, &sqr);
+    CHECK(secp256k1_num_jacobi(&n, &order) == -1);
+}
+
 void run_num_smalltests(void) {
     int i;
     for (i = 0; i < 100*count; i++) {
         test_num_negate();
         test_num_add_sub();
+        test_num_jacobi();
     }
 }