tommath/bn_mp_jacobi.c

113 lines
2.4 KiB
C
Raw Normal View History

#include <tommath_private.h>
2004-10-29 18:07:18 -04:00
#ifdef BN_MP_JACOBI_C
2003-02-28 11:08:34 -05:00
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
2003-08-04 21:24:44 -04:00
* LibTomMath is a library that provides multiple-precision
2003-02-28 11:08:34 -05:00
* integer arithmetic as well as number theoretic functionality.
*
2003-08-04 21:24:44 -04:00
* The library was designed directly after the MPI library by
2003-02-28 11:08:34 -05:00
* Michael Fromberger but has been written from scratch with
* additional optimizations in place.
*
* The library is free for all purposes without any express
* guarantee it works.
*
* Tom St Denis, tstdenis82@gmail.com, http://libtom.org
2003-02-28 11:08:34 -05:00
*/
2003-05-17 08:33:54 -04:00
/* computes the jacobi c = (a | n) (or Legendre if n is prime)
2003-02-28 11:08:34 -05:00
* HAC pp. 73 Algorithm 2.149
* HAC is wrong here, as the special case of (0 | 1) is not
* handled correctly.
2003-02-28 11:08:34 -05:00
*/
2015-11-25 15:24:13 -05:00
int mp_jacobi (mp_int * a, mp_int * n, int *c)
2003-02-28 11:08:34 -05:00
{
2003-07-02 11:39:39 -04:00
mp_int a1, p1;
int k, s, r, res;
2003-02-28 11:09:08 -05:00
mp_digit residue;
2003-02-28 11:08:34 -05:00
2015-11-25 15:24:13 -05:00
/* if n <= 0 return MP_VAL */
if (mp_cmp_d(n, 0) != MP_GT) {
2003-08-04 21:24:44 -04:00
return MP_VAL;
}
/* step 1. handle case of a == 0 */
2015-10-25 11:34:43 -04:00
if (mp_iszero (a) == MP_YES) {
/* special case of a == 0 and n == 1 */
if (mp_cmp_d (n, 1) == MP_EQ) {
*c = 1;
} else {
*c = 0;
}
return MP_OKAY;
2003-02-28 11:08:34 -05:00
}
/* step 2. if a == 1, return 1 */
if (mp_cmp_d (a, 1) == MP_EQ) {
*c = 1;
return MP_OKAY;
}
/* default */
2003-08-04 21:24:44 -04:00
s = 0;
2003-02-28 11:08:34 -05:00
2003-07-02 11:39:39 -04:00
/* step 3. write a = a1 * 2**k */
2003-02-28 11:08:34 -05:00
if ((res = mp_init_copy (&a1, a)) != MP_OKAY) {
return res;
}
2003-07-02 11:39:39 -04:00
if ((res = mp_init (&p1)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_A1;
2003-02-28 11:08:34 -05:00
}
2003-08-04 21:24:44 -04:00
/* divide out larger power of two */
k = mp_cnt_lsb(&a1);
if ((res = mp_div_2d(&a1, k, &a1, NULL)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_P1;
2003-02-28 11:08:34 -05:00
}
/* step 4. if e is even set s=1 */
2003-07-02 11:39:39 -04:00
if ((k & 1) == 0) {
2003-02-28 11:08:34 -05:00
s = 1;
} else {
2003-07-02 11:39:39 -04:00
/* else set s=1 if p = 1/7 (mod 8) or s=-1 if p = 3/5 (mod 8) */
2015-11-25 15:24:13 -05:00
residue = n->dp[0] & 7;
2003-02-28 11:08:34 -05:00
if (residue == 1 || residue == 7) {
s = 1;
} else if (residue == 3 || residue == 5) {
s = -1;
}
}
2003-07-02 11:39:39 -04:00
/* step 5. if p == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */
2015-11-25 15:24:13 -05:00
if ( ((n->dp[0] & 3) == 3) && ((a1.dp[0] & 3) == 3)) {
2003-07-02 11:39:39 -04:00
s = -s;
2003-02-28 11:08:34 -05:00
}
/* if a1 == 1 we're done */
if (mp_cmp_d (&a1, 1) == MP_EQ) {
*c = s;
} else {
/* n1 = n mod a1 */
2015-11-25 15:24:13 -05:00
if ((res = mp_mod (n, &a1, &p1)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_P1;
2003-02-28 11:08:34 -05:00
}
2003-07-02 11:39:39 -04:00
if ((res = mp_jacobi (&p1, &a1, &r)) != MP_OKAY) {
2004-12-22 21:40:37 -05:00
goto LBL_P1;
2003-02-28 11:08:34 -05:00
}
*c = s * r;
}
/* done */
res = MP_OKAY;
2004-12-22 21:40:37 -05:00
LBL_P1:mp_clear (&p1);
LBL_A1:mp_clear (&a1);
2003-02-28 11:08:34 -05:00
return res;
}
2004-10-29 18:07:18 -04:00
#endif
2005-08-01 12:37:28 -04:00
/* $Source$ */
/* $Revision$ */
/* $Date$ */