home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / crypl200.zip / BNLIB / JACOBI.C < prev    next >
C/C++ Source or Header  |  1996-05-16  |  2KB  |  68 lines
  1. /*
  2.  * Compute the Jacobi symbol (small prime case only).
  3.  *
  4.  * Copyright (c) 1995  Colin Plumb.  All rights reserved.
  5.  * For licensing and other legal details, see the file legal.c.
  6.  */
  7. #include "bn.h"
  8. #include "jacobi.h"
  9.  
  10. /*
  11.  * For a small (usually prime, but not necessarily) prime p,
  12.  * compute Jacobi(p,bn), which is -1, 0 or +1, using the following rules:
  13.  * Jacobi(x, y) = Jacobi(x mod y, y)
  14.  * Jacobi(0, y) = 0
  15.  * Jacobi(1, y) = 1
  16.  * Jacobi(2, y) = 0 if y is even, +1 if y is +/-1 mod 8, -1 if y = +/-3 mod 8
  17.  * Jacobi(x1*x2, y) = Jacobi(x1, y) * Jacobi(x2, y) (used with x1 = 2 & x1 = 4)
  18.  * If x and y are both odd, then
  19.  * Jacobi(x, y) = Jacobi(y, x) * (-1 if x = y = 3 mod 4, +1 otherwise)
  20.  */
  21. int
  22. bnJacobiQ(unsigned p, struct BigNum const *bn)
  23. {
  24.     int j = 1;
  25.     unsigned u = bnLSWord(bn);
  26.  
  27.     if (!(u & 1))
  28.         return 0;    /* Don't *do* that */
  29.  
  30.     /* First, get rid of factors of 2 in p */
  31.     while ((p & 3) == 0)
  32.         p >>= 2;
  33.     if ((p & 1) == 0) {
  34.         p >>= 1;
  35.         if ((u ^ u>>1) & 2)
  36.             j = -j;        /* 3 (011) or 5 (101) mod 8 */
  37.     }
  38.     if (p == 1)
  39.         return j;
  40.     /* Then, apply quadratic reciprocity */
  41.     if (p & u & 2)    /* p = u = 3 (mod 4? */
  42.         j = -j;
  43.     /* And reduce u mod p */
  44.     u = bnModQ(bn, p);
  45.  
  46.     /* Now compute Jacobi(u,p), u < p */
  47.     while (u) {
  48.         while ((u & 3) == 0)
  49.             u >>= 2;
  50.         if ((u & 1) == 0) {
  51.             u >>= 1;
  52.             if ((p ^ p>>1) & 2)
  53.                 j = -j;    /* 3 (011) or 5 (101) mod 8 */
  54.         }
  55.         if (u == 1)
  56.             return j;
  57.         /* Now both u and p are odd, so use quadratic reciprocity */
  58.         if (u < p) {
  59.             unsigned t = u; u = p; p = t;
  60.             if (u & p & 2)    /* u = p = 3 (mod 4? */
  61.                 j = -j;
  62.         }
  63.         /* Now u >= p, so it can be reduced */
  64.         u %= p;
  65.     }
  66.     return 0;
  67. }
  68.