Home | History | Annotate | Line # | Download | only in libtommath 
      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 (at) 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  */
     23 int 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;
     99 LBL_P1:mp_clear (&p1);
    100 LBL_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