Let m = u×2k, where u is odd. For notational convenience, let e (the even part) = 2k. Let y = xu and b = cu, giving ye = ±1 and be = -1.
If r prime divides n, b2e = 1 mod r, and r-1 is divisible by 2e. Thus all divisors of n are 1 mod 2e.
Assume d1 and d2 are divisors of n, and [y\d1] and [y\d2] agree with y(d1-1)/2 and y(d2-1)/2 mod n. Consider d3 = d1×d2. Jacobi symbols are multiplied, so [y\d3] = [y\d1]×[y\d2]. Write d1 = 2ev1+1 and d2 = 2ev2+1. The new exponent becomes (d1d2-1)/2, or 2e2v1v2+ev1+ev2. The last two terms are (d1-1)/2 and (d2-1)/2. Bring these together and you are multiplying x(d1-1)/2 by x(d2-1)/2. The first term in the exponent looks like d1-1 times ev2. Raise x to this power and get 1, which doesn't change a thing. If d1 and d2 match, then so does d3.
For each prime r, [y\r] = y(r-1)/2 mod r. This expression is also a power of ye, and is ±1 mod n, so [y\r] = y(r-1)/2 mod n. This holds for all primes r dividing n. Use the previous result to combine divisors, until [y\n] = y(n-1)/2 mod n.
Great, but we're really interested in x. Replace y with xu and get this.
[xu\n] = xu×(n-1)/2 mod n
The left side is [x\n]u, and since u is odd, this is the same as [x\n]. The right side is also an expression that equals ±1 raised to the u, so that doesn't change anything either. At the end of the day, [x\n] = x(n-1)/2. The jacobi symbol agrees with the exponential. All we need is c and m, as described above.
If x(n-1)/2 ≠ ±1, the pseudoprime test fails and n is composite.
let u be the odd portion of n-1. If xu ≠ 1, and its successive squares eventually produce 1, without passing through -1, 1 has a square root different from ±1, and n is composite.
Raise x to the u. If this is 1 it will remain 1 all the way up to xn-1. The exponential matches the jacobi symbol by condition 4. If xu is -1, invoke condition 3. If xu is anything else, keep squaring until you reach -1. When you find it, invoke condition 3. If you don't find it, n is composite.
bool StrongPrimeTest(x, n)
{
/* remove powers of 2 from n-1 */
for(k=0, u=n-1; even(u); u/=2) ++k;
x = (x^u) % n; /* x to the u mod n */
if(x == 1 || x == -1) return true; /* could be prime */
--k; /* don't need to go all the way up to n-1 */
while(k > 0) { /* at most log2(n) iterations */
x = (x^2) % n;
if(x == 1) return false; /* composite */
if(x == -1) return true;
--k;
}
/* We didn't even get +-1 for x^((n-1)/2). */
return false;
}
Let u be the largest odd divisor of n-1. Let k be the maximum for which x2k = -1. Note that k is at least 0, since -1 to the first power is indeed -1.
Let m = u×2k. Suppose x is not an mth root of ±1. Raise x to the u and we don't get ±1, so we keep going. Square, and square, and square, and we never get ±1, else x would be an mth root of ±1. The test fails. We only need look at the mth roots of ±1.
The aforementioned roots form a group, which I will call G. Since (-1)u = -1, and some x to the 2k = -1, G includes an element x that is a primitive mth root of -1.
Let H be all the (2m)th roots of 1. Clearly G is a subgroup of H. The homomorphism xm takes H into the square roots of 1, and it takes G onto ±1.
Split n into a product of s prime powers pibi. Everything in G becomes an mth root of ±1 mod pibi. A typical cycle of units, of order pb-1×(p-1) has an mth root of -1, inherited from G. This implies a primitive root of order 2×2k. Since p is odd, 2×2k divides p-1. In other words, p = 1 mod 2×2k. Multiply primes together and n = 1 mod 2×2k. this places an upper bound on k. It follows that 2m is at most n-1.
For each prime power, independently select an mth root of 1 or -1. When raised to the m power, we have all combinations of 1 and -1. This gives 2s square roots of 1. In other words, xm maps H onto the square roots of 1. Only two of these come out ±1 mod n, and this is the image of G. Running the homomorphism backwards, the index of G in H is 2s-1. when outside of G, the test fails, so if s > 2, the test fails 3/4 of the time, and we're done.
In the worst case, 2m = n-1, but the (n-1)th roots of 1 form a subgroup inside the units mod n. Call this group W, and let V be the group of units mod n. The efficiency of the algorithm is multiplied by the index of W in V.
Suppose s = 2, whence n cannot be a carmichael number. There is an x relatively prime to n that is not in W. The index of W in V is at least 2, and when combined with the index of G in H, we get 4, and we're done.
Finally, let s = 1, so that n = pb, and V is cyclic of order pb-1×(p-1). since p and n-1 are coprime, W is contained in the latter cycle. At worst, the entire cycle is W, whence the index of W in V is pb-1. this is less than 4 only when n = 32. (Not that anybody's going to work hard to see if 9 is prime.) Even in this case, we can still save the day. Bring in the elements that are not units, namely 3 and 6, and 6 elements out of 8 fail the pseudoprime test.
For any composite n, at least 3/4 of the values of x between 1 and n-1 prove that n is composite. For most values of n, the ratio is much better. However, when n is a carmichael number based on 3 primes, G has index 4 in H, and W has index 1 in V, giving a combined index of 4. It is conjectured that there are infinitely many numbers of this form, and if true, the test cannot be more than 3/4 accurate, so it looks like we're done.
Let b = xn-1, already computed. If b = 0 then x contains all the primes of n, and gcd(x,n) provides a factor of n. If b = 1 then x yields a nontrivial square root of 1. add 1 to this square root, giving a number that has some, but not all, of the primes of n. Use the gcd algorithm to find the common factor. Otherwise, run gcd(b-1,n). When n = pb, raising to the n-1 begins by raising to the p-1, which produces a value equal to 1 mod p. The result is not 1 mod n, so the gcd will pull out some, but not all, of the factors of p.