Merge branch 'fix/jacobi' into develop

This closes #31

bn_mp_jacobi.c

@@ -17,22 +17,29 @@
 
 /* computes the jacobi c = (a | n) (or Legendre if n is prime)
  * HAC pp. 73 Algorithm 2.149
+ * HAC is wrong here, as the special case of (0 | 1) is not
+ * handled correctly.
  */
-int mp_jacobi (mp_int * a, mp_int * p, int *c)
+int mp_jacobi (mp_int * a, mp_int * n, int *c)
 {
   mp_int  a1, p1;
   int     k, s, r, res;
   mp_digit residue;
 
-  /* if p <= 0 return MP_VAL */
+  /* if n <= 0 return MP_VAL */
-  if (mp_cmp_d(p, 0) != MP_GT) {
+  if (mp_cmp_d(n, 0) != MP_GT) {
      return MP_VAL;
   }
 
-  /* step 1.  if a == 0, return 0 */
+  /* step 1. handle case of a == 0 */
   if (mp_iszero (a) == MP_YES) {
-    *c = 0;
+     /* special case of a == 0 and n == 1 */
-    return MP_OKAY;
+     if (mp_cmp_d (n, 1) == MP_EQ) {
+       *c = 1;
+     } else {
+       *c = 0;
+     }
+     return MP_OKAY;
   }
 
   /* step 2.  if a == 1, return 1 */
@@ -64,7 +71,7 @@ int mp_jacobi (mp_int * a, mp_int * p, int *c)
     s = 1;
   } else {
     /* else set s=1 if p = 1/7 (mod 8) or s=-1 if p = 3/5 (mod 8) */
-    residue = p->dp[0] & 7;
+    residue = n->dp[0] & 7;
 
     if (residue == 1 || residue == 7) {
       s = 1;
@@ -74,7 +81,7 @@ int mp_jacobi (mp_int * a, mp_int * p, int *c)
   }
 
   /* step 5.  if p == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */
-  if ( ((p->dp[0] & 3) == 3) && ((a1.dp[0] & 3) == 3)) {
+  if ( ((n->dp[0] & 3) == 3) && ((a1.dp[0] & 3) == 3)) {
     s = -s;
   }
 
@@ -83,7 +90,7 @@ int mp_jacobi (mp_int * a, mp_int * p, int *c)
     *c = s;
   } else {
     /* n1 = n mod a1 */
-    if ((res = mp_mod (p, &a1, &p1)) != MP_OKAY) {
+    if ((res = mp_mod (n, &a1, &p1)) != MP_OKAY) {
       goto LBL_P1;
     }
     if ((res = mp_jacobi (&p1, &a1, &r)) != MP_OKAY) {

demo/demo.c

@@ -12,7 +12,7 @@
  * Configuration
  */
 #ifndef LTM_DEMO_TEST_VS_MTEST
-#define LTM_DEMO_TEST_VS_MTEST 1
+#define LTM_DEMO_TEST_VS_MTEST 0
 #endif
 
 #ifndef LTM_DEMO_TEST_REDUCE_2K_L
@@ -114,6 +114,16 @@ struct mp_sqrtmod_prime_st sqrtmod_prime[] = {
       { 7, 9, 4 },
       { 113, 2, 62 }
 };
+struct mp_jacobi_st {
+   unsigned long n;
+   int c[16];
+};
+struct mp_jacobi_st jacobi[] = {
+      { 3, {  1, -1,  0,  1, -1,  0,  1, -1,  0,  1, -1,  0,  1, -1,  0,  1 } },
+      { 5, {  0,  1, -1, -1,  1,  0,  1, -1, -1,  1,  0,  1, -1, -1,  1,  0 } },
+      { 7, {  1, -1,  1, -1, -1,  0,  1,  1, -1,  1, -1, -1,  0,  1,  1, -1 } },
+      { 9, { -1,  1,  0,  1,  1,  0,  1,  1,  0,  1,  1,  0,  1,  1,  0,  1 } },
+};
 
 char cmd[4096], buf[4096];
 int main(void)
@@ -186,6 +196,35 @@ int main(void)
    mp_add_d(&a, 1, &b);
    mp_add_d(&a, 6, &b);
 
+
+   mp_set_int(&a, 0);
+   mp_set_int(&b, 1);
+   if ((ix = mp_jacobi(&a, &b, &i)) != MP_OKAY) {
+      printf("Failed executing mp_jacobi(0 | 1) %s.\n", mp_error_to_string(ix));
+      return EXIT_FAILURE;
+   }
+   if (i != 1) {
+      printf("Failed trivial mp_jacobi(0 | 1) %d != 1\n", i);
+      return EXIT_FAILURE;
+   }
+   for (cnt = 0; cnt < (int)(sizeof(jacobi)/sizeof(jacobi[0])); ++cnt) {
+      mp_set_int(&b, jacobi[cnt].n);
+      /* only test positive values of a */
+//      for (n = -5; n <= 10; ++n) {
+      for (n = 0; n <= 10; ++n) {
+         mp_set_int(&a, abs(n));
+         if (n < 0) mp_neg(&a, &a);
+         if ((ix = mp_jacobi(&a, &b, &i)) != MP_OKAY) {
+            printf("Failed executing mp_jacobi(%d | %lu) %s.\n", n, jacobi[cnt].n, mp_error_to_string(ix));
+            return EXIT_FAILURE;
+         }
+         if (i != jacobi[cnt].c[n + 5]) {
+            printf("Failed trivial mp_jacobi(%d | %lu) %d != %d\n", n, jacobi[cnt].n, i, jacobi[cnt].c[n + 5]);
+            return EXIT_FAILURE;
+         }
+      }
+   }
+
    // test montgomery
    printf("Testing: montgomery...\n");
    for (i = 1; i <= 10; i++) {