Home | History | Annotate | Download | only in libtommath 
      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 (at) gmail.com, http://math.libtomcrypt.com
     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
    102 
    103 /* $Source: /cvs/libtom/libtommath/bn_mp_jacobi.c,v $ */
    104 /* $Revision: 1.3 $ */
    105 /* $Date: 2006/03/31 14:18:44 $ */
    106