libs/libtommath/bn_mp_jacobi.c

   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  */
  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.4 $ */
 105 /* $Date: 2006/12/28 01:25:13 $ */