crypto/external/bsd/heimdal/dist/lib/hcrypto/libtommath/bn_mp_jacobi.c

   1 /*      $NetBSD: bn_mp_jacobi.c,v 1.1.1.2 2014/04/24 12:45:31 pettai 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  */
  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  */