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