Merge branch 'fix/jacobi' into develop This closes #31
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133
diff --git a/bn_mp_jacobi.c b/bn_mp_jacobi.c
index 160ac8f..6a1155d 100644
--- a/bn_mp_jacobi.c
+++ b/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 (mp_cmp_d(p, 0) != MP_GT) {
+ /* if n <= 0 return MP_VAL */
+ 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;
- return MP_OKAY;
+ /* special case of a == 0 and n == 1 */
+ 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) {
diff --git a/demo/demo.c b/demo/demo.c
index a5ac674..45be283 100644
--- a/demo/demo.c
+++ b/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++) {