Mercurial > dropbear

comparison bn_mp_jacobi.c @ 282:91fbc376f010 libtommath-orig libtommath-0.35

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import of libtommath 0.35 From ltm-0.35.tar.bz2 SHA1 of 3f193dbae9351e92d02530994fa18236f7fde01c
author Matt Johnston <matt@ucc.asn.au>
date Wed, 08 Mar 2006 13:16:18 +0000
parents
children 97db060d0ef5
comparison
equal deleted inserted replaced
-1:000000000000 282:91fbc376f010
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, [email protected], http://math.libtomcrypt.org
16 */
17
18 /* computes the jacobi c = (a | n) (or Legendre if n is prime)
19 * HAC pp. 73 Algorithm 2.149
20 */
21 int 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;
97 LBL_P1:mp_clear (&p1);
98 LBL_A1:mp_clear (&a1);
99 return res;
100 }
101 #endif