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 */ 23int 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; 99LBL_P1:mp_clear (&p1); 100LBL_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