1 /* $NetBSD: bn_mp_jacobi.c,v 1.2 2017/01/28 21:31:47 christos Exp $ */
2
3 #include <tommath.h>
4 #ifdef BN_MP_JACOBI_C
5 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6 *
7 * LibTomMath is a library that provides multiple-precision
8 * integer arithmetic as well as number theoretic functionality.
9 *
10 * The library was designed directly after the MPI library by
11 * Michael Fromberger but has been written from scratch with
12 * additional optimizations in place.
13 *
14 * The library is free for all purposes without any express
15 * guarantee it works.
16 *
17 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
18 */
19
20 /* computes the jacobi c = (a | n) (or Legendre if n is prime)
21 * HAC pp. 73 Algorithm 2.149
22 */
mp_jacobi(mp_int * a,mp_int * p,int * c)23 int mp_jacobi (mp_int * a, mp_int * p, int *c)
24 {
25 mp_int a1, p1;
26 int k, s, r, res;
27 mp_digit residue;
28
29 /* if p <= 0 return MP_VAL */
30 if (mp_cmp_d(p, 0) != MP_GT) {
31 return MP_VAL;
32 }
33
34 /* step 1. if a == 0, return 0 */
35 if (mp_iszero (a) == 1) {
36 *c = 0;
37 return MP_OKAY;
38 }
39
40 /* step 2. if a == 1, return 1 */
41 if (mp_cmp_d (a, 1) == MP_EQ) {
42 *c = 1;
43 return MP_OKAY;
44 }
45
46 /* default */
47 s = 0;
48
49 /* step 3. write a = a1 * 2**k */
50 if ((res = mp_init_copy (&a1, a)) != MP_OKAY) {
51 return res;
52 }
53
54 if ((res = mp_init (&p1)) != MP_OKAY) {
55 goto LBL_A1;
56 }
57
58 /* divide out larger power of two */
59 k = mp_cnt_lsb(&a1);
60 if ((res = mp_div_2d(&a1, k, &a1, NULL)) != MP_OKAY) {
61 goto LBL_P1;
62 }
63
64 /* step 4. if e is even set s=1 */
65 if ((k & 1) == 0) {
66 s = 1;
67 } else {
68 /* else set s=1 if p = 1/7 (mod 8) or s=-1 if p = 3/5 (mod 8) */
69 residue = p->dp[0] & 7;
70
71 if (residue == 1 || residue == 7) {
72 s = 1;
73 } else if (residue == 3 || residue == 5) {
74 s = -1;
75 }
76 }
77
78 /* step 5. if p == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */
79 if ( ((p->dp[0] & 3) == 3) && ((a1.dp[0] & 3) == 3)) {
80 s = -s;
81 }
82
83 /* if a1 == 1 we're done */
84 if (mp_cmp_d (&a1, 1) == MP_EQ) {
85 *c = s;
86 } else {
87 /* n1 = n mod a1 */
88 if ((res = mp_mod (p, &a1, &p1)) != MP_OKAY) {
89 goto LBL_P1;
90 }
91 if ((res = mp_jacobi (&p1, &a1, &r)) != MP_OKAY) {
92 goto LBL_P1;
93 }
94 *c = s * r;
95 }
96
97 /* done */
98 res = MP_OKAY;
99 LBL_P1:mp_clear (&p1);
100 LBL_A1:mp_clear (&a1);
101 return res;
102 }
103 #endif
104
105 /* Source: /cvs/libtom/libtommath/bn_mp_jacobi.c,v */
106 /* Revision: 1.4 */
107 /* Date: 2006/12/28 01:25:13 */
108