Perrin Primality Tests
by Dana Jacobsen
8 August 2016
The Perrin sequence (OEIS A001608) is:
It can be generalized to other third-order recurrences, typically using the two parameters r and s. This page considers just the original Perrin sequence, where r = 0 and s = -1. Hence Δ = -23 and f(x) = x3 - x - 1.
The use for primality tests (technically compositeness tests) comes from the observation that n | P(n) for all primes and relatively few composites. Adams and Shanks AS82 is the primary reference for primality testing with the Perrin sequence. They also showed that if the sequence is reversed, then n | P(-n)+1. These may be equivalently stated as P(n) = P(1) mod n and P(-n) = P(-1) mod n.
Composites which pass these tests are called Perrin pseudoprimes. Some authors Arno91;Grantham00 restrict the composite inputs to odd numbers, while others AS82;KSW86 do not. There are at least four variations that have been published, which are described below.
The signature mod m of an integer n is the 6-tuple ( P(-n-1),P(-n),P(-n+1), P(n-1),P(n),P(n+1) ) mod m. Typically m = n, and this is just called the signature of n. Each acceptable signature is of the form (__ P(-1) __ __ P(1) __). There are efficient ways to calculate this in O(log n) time, as well as simple divisibility tests using the periods modulo small primes that can quickly discard many composites without having to calculate the sequence.
The congruences for Q- and I-type signatures come from Adams and Shanks AS82. The alternate congurences are useful for the more general case of non-Perrin cubic recurrences, and are the ones used by Arno Arno91 and Grantham Grantham00.
Barring the quadratic form for I-types used by Adams and Shanks, which is not used by other authors, each variation above produces pseudoprimes which are a subset of the previous.
Unrestricted: The "unrestricted" Perrin pseudoprimes use the simple P(n) = 0 mod n relation, which is sometimes expressed as Trace(Mn) = Trace(M) mod n, with M=[0 1 0 ; 0 0 1 ; 1 1 0]. The first pseudoprime is 271441. There are 1700 pseudoprimes less than 1014. These comprise the sequence A013998.
Minimal restricted: Adding the P(-n) = -1 mod n condition is called the "minimal test" in KSW86 and the "strong Perrin test" in Wagon99 . The first pseudoprime is 27664033. There are 942 pseudoprimes less than 1014. These comprise the sequence A018187.
Adams/Shanks: Adams and Shanks (1982) note the signature must be one of only three distinct forms, each with specific allowable Jacobi symbols. They further give a quadratic form equivalence for I-type signatures, which is intentionally removed by later authors. There are 936 pseudoprimes less than 1014, with the six discarded from the previous set being:
These comprise the sequence A275612.
Grantham: Arno (1991) and Grantham (2000) describe a concise variation of the AS82 test. They explicitly limit pseudoprimes to odd composites and disallow (-23|n) = 0. Note that these changes require a primality test to ensure the primes 2 and 23 are returned correctly. There are 931 pseudoprimes less than 1014, with the five discarded from the previous set being:
These comprise the sequence A275613.
Pari/GP code for each of the four types. This is example code, and while not inefficient, is certainly not optimal. Perl/ntheory has a custom function, is_perrin_pseudoprime, which for 64-bit inputs is 30-80x faster.
The matrix method used below takes at least 27×log2(n) mulmods for a half-signature, double that for a full signature. Using the doubling rule from Adams and Shanks reduces this to 6×log2(n) mulmods for a full signature.. For comparison, well-implemented Lucas tests such as the ones in Pari/GP and Perl/ntheory take 2×log2(n) mulmods.
perrin0(n) = { lift(trace(Mod([0,1,0; 0,0,1; 1,1,0],n)^n)) == 0; }
perrin1(n) = { lift(trace(Mod([0,1,0; 0,0,1; 1,1,0],n)^n)) == 0 && lift(trace(Mod([0,1,0; 0,0,1; 1,0,-1],n)^n)) == n-1; }
perrin2(n) = {
my(M,L,S,j,A,B,C,D);
M=Mod( [0,1,0; 0,0,1; 1,1,0], n)^n;
L=Mod( [0,1,0; 0,0,1; 1,0,-1], n)^n;
S=[ 3*L[3,2]-L[3,3], 3*L[2,2]-L[2,3], 3*L[1,2]-L[1,3], \
3*M[3,1]+2*M[3,3], 3*M[1,1]+2*M[1,3], 3*M[2,1]+2*M[2,3] ];
if (S[5] != 0 || S[2] != n-1,return(0));
j = kronecker(-23,n);
if (j == -1, B=S[3];A=1+3*B-B^2;C=3*B^2-2; if(S[1]==A && S[3]==B && S[4]==B && S[6] == C && B != 3 && B^3-B==1, return(1), return(0)));
if (S[1] == 1 && S[3] == 3 && S[4] == 3 && S[6] == 2, return(1));
if (j == 1 && S[1] == 0 && S[6] == n-1 && S[3] != S[4] && S[3]+S[4] == n-3 && (S[3]-S[4])^2 == Mod(-23,n), return(1));
return(0);
}
perrin3(n) = {
my(M,L,S,j,A,B,C,D);
if(n==2||n==23,return(1));
if(n%2==0,return(0));
M=Mod( [0,1,0; 0,0,1; 1,1,0], n)^n;
L=Mod( [0,1,0; 0,0,1; 1,0,-1], n)^n;
S=[ 3*L[3,2]-L[3,3], 3*L[2,2]-L[2,3], 3*L[1,2]-L[1,3], \
3*M[3,1]+2*M[3,3], 3*M[1,1]+2*M[1,3], 3*M[2,1]+2*M[2,3] ];
if (S[5] != 0 || S[2] != n-1,return(0));
j = kronecker(-23,n);
if (j == 0,return(0));
if (j == -1, B=S[3];A=1+3*B-B^2;C=3*B^2-2; if(S[1]==A && S[3]==B && S[4]==B && S[6] == C && B != 3 && B^3-B==1, return(1), return(0)));
if (S[1] == 1 && S[3] == 3 && S[4] == 3 && S[6] == 2, return(1));
if (S[1] == 0 && S[6] == n-1 && S[3] != S[4] && S[3]+S[4] == n-3 && (S[3]-S[4])^2 == Mod(-23,n), return(1));
return(0);
}
forcomposite(n=1,1e6,if(perrin0(n),print(n)))
In practice the Perrin sequence is not used for primality testing. They are more computationally expensive to calculate than Lucas sequences and overlap significantly with strong pseudoprimes. This makes them undesirable in both performance and strength compared to the BPSW test which is the most commonly used test, or the Frobenius test with quadratic polynomials..
In practical primality tests, some amount of trial division is done to quickly remove many composites. For the purposes of finding Perrin pseudoprimes, this is not done. However, it is possible to use the periods of the Perrin sequence modulo small primes, as described in Kurtz et al. KSW86, to quickly remove many composites that can not satisfy P(n) = 0 mod n. This filters out 70-90% of composites without having the run the test.
The naive and extremely inefficient way to compute P(n) mod n is to find the n-th Perrin number by iterating the rule P(n) = P(n-2) + P(n-3), then checking divisibility. The first observation is that we can perform all math mod n, hence there is no need to deal with big numbers if the input n is of reasonable size. This still leaves O(n) operations, which is unacceptable.
It can be shown that modular exponentiation can be used to generate P(n), P(n-1), and P(n+1) using the square matrix [0,1,0; 0,0,1; 1,1,0]. This method can also generate the negative values as well using the matrix [0,1,0; 0,0,1; 1,0,-1]. Modular exponentiation can be done in O(log(n)) steps (with 1 or 2 modular 3x3 matrix multiplications in each step). This is not an uncommon operation to have available in some languages/packages (e.g. Pari/GP makes this trivial).
Using the doubling method outlined in Adams and Shanks (1982), the full signature can be constructed with only 6 mulmods per bit of the input. Matrix exponentiation takes 23-27 × (1 + number of 1 bits) mulmods per bit, and two are required to get the full signature (albeit almost all composites are rejected after the first).
Colin Wright has an excellent blog post about these, with iterated examples.
Programming Praxis wrote an exercise about them: Perrin Pseudoprimes.