Add Jacobi symbol test via GMP

Also add native Jacobi symbol test (Andrew)

Rebased-by: Andrew Poelstra

Author: peterdettman

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 < 2000; 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);
 
@@ -56,6 +59,9 @@ static void secp256k1_num_shift(secp256k1_num *r, int bits);
 /** Check whether a number is zero. */
 static int secp256k1_num_is_zero(const secp256k1_num *a);
 
+/** Check whether a number is one. */
+static int secp256k1_num_is_one(const secp256k1_num *a);
+
 /** Check whether a number is strictly negative. */
 static int secp256k1_num_is_neg(const secp256k1_num *a);
 

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

@@ -102,6 +102,16 @@ SECP256K1_INLINE static int secp256k1_num_is_zero(const secp256k1_num *a) {
     return 1;
 }
 
+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_neg(const secp256k1_num *a) {
     return a->data[NUM_N_WORDS - 1] >> (NUM_WORD_WIDTH - 1);
 }
@@ -599,4 +609,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 Jacobi 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

@@ -472,6 +472,8 @@ void test_num_negate(void) {
 }
 
 void test_num_add_sub(void) {
+    int i;
+    secp256k1_scalar s;
     secp256k1_num n1;
     secp256k1_num n2;
     secp256k1_num n1p2, n2p1, n1m2, n2m1;
@@ -498,6 +500,18 @@ void test_num_add_sub(void) {
     CHECK(!secp256k1_num_eq(&n2p1, &n1));
     secp256k1_num_sub(&n2p1, &n2p1, &n2); /* n2p1 = R2 + R1 - R2 = R1 */
     CHECK(secp256k1_num_eq(&n2p1, &n1));
+
+    /* check is_one */
+    secp256k1_scalar_set_int(&s, 1);
+    secp256k1_scalar_get_num(&n1, &s);
+    CHECK(secp256k1_num_is_one(&n1));
+    /* check that 2^n + 1 is never 1 */
+    secp256k1_scalar_get_num(&n2, &s);
+    for (i = 0; i < 250; ++i) {
+        secp256k1_num_add(&n1, &n1, &n1);    /* n1 *= 2 */
+        secp256k1_num_add(&n1p2, &n1, &n2);  /* n1p2 = n1 + 1 */
+        CHECK(!secp256k1_num_is_one(&n1p2));
+    }
 }
 
 void test_num_mod(void) {
@@ -531,12 +545,74 @@ void test_num_mod(void) {
     CHECK(secp256k1_num_is_zero(&n));
 }
 
+void test_num_jacobi(void) {
+    secp256k1_scalar sqr;
+    secp256k1_scalar small;
+    secp256k1_scalar five;  /* five is not a quadratic residue */
+    secp256k1_num order, n;
+    int i;
+    /* squares mod 5 are 1, 4 */
+    const int jacobi5[10] = { 0, 1, -1, -1, 1, 0, 1, -1, -1, 1 };
+
+    /* check some small values with 5 as the order */
+    secp256k1_scalar_set_int(&five, 5);
+    secp256k1_scalar_get_num(&order, &five);
+    for (i = 0; i < 10; ++i) {
+        secp256k1_scalar_set_int(&small, i);
+        secp256k1_scalar_get_num(&n, &small);
+        CHECK(secp256k1_num_jacobi(&n, &order) == jacobi5[i]);
+    }
+
+    /** test large values with 5 as group order */
+    secp256k1_scalar_get_num(&order, &five);
+    /* we first need a scalar which is not a multiple of 5 */
+    do {
+        secp256k1_num fiven;
+        random_scalar_order_test(&sqr); 
+        secp256k1_scalar_get_num(&fiven, &five);
+        secp256k1_scalar_get_num(&n, &sqr);
+        secp256k1_num_mod(&n, &fiven);
+    } while (secp256k1_num_is_zero(&n));
+    /* next force it to be a residue. 2 is a nonresidue mod 5 so we can
+     * just multiply by two, i.e. add the number to itself */
+    if (secp256k1_num_jacobi(&n, &order) == -1) {
+        secp256k1_num_add(&n, &n, &n);
+    }
+
+    /* test residue */
+    CHECK(secp256k1_num_jacobi(&n, &order) == 1);
+    /* test nonresidue */
+    secp256k1_num_add(&n, &n, &n);
+    CHECK(secp256k1_num_jacobi(&n, &order) == -1);
+
+    /** test with secp group order as order */
+    secp256k1_scalar_order_get_num(&order);
+    random_scalar_order_test(&sqr); 
+    secp256k1_scalar_sqr(&sqr, &sqr);
+    /* 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);
+    /* test multiple of the order*/
+    CHECK(secp256k1_num_jacobi(&order, &order) == 0);
+
+    /* check one less than the order */
+    secp256k1_scalar_set_int(&small, 1);
+    secp256k1_scalar_get_num(&n, &small);
+    secp256k1_num_sub(&n, &order, &n);
+    CHECK(secp256k1_num_jacobi(&n, &order) == 1);  /* sage confirms this is 1 */
+}
+
 void run_num_smalltests(void) {
     int i;
     for (i = 0; i < 100*count; i++) {
         test_num_negate();
         test_num_add_sub();
         test_num_mod();
+        test_num_jacobi();
     }
 }
 
@@ -677,6 +753,8 @@ void scalar_test(void) {
             secp256k1_scalar_inverse(&inv, &inv);
             /* Inverting one must result in one. */
             CHECK(secp256k1_scalar_is_one(&inv));
+            secp256k1_scalar_get_num(&invnum, &inv);
+            CHECK(secp256k1_num_is_one(&invnum));
         }
     }