| Alexandre Lision | 7fd5d3d | 2013-12-04 13:06:40 -0500 | [diff] [blame] | 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 | } |