Jump to content

Perrin number: Difference between revisions

From Wikipedia, the free encyclopedia
Content deleted Content added
m fixed lint errors – missing end tag
An additional test and a C program
Line 151: Line 151:


===Pseudocode===
===Pseudocode===
The 1982 Adams and Shanks {{math|O(log n)}} Strong Perrin primality test.<ref>{{harvtxt|Adams|Shanks|1982|p=265, 269-270}}</ref>
The 1982 Adams and Shanks {{math|O(log n)}} Perrin primality test.<ref>{{harvtxt|Adams|Shanks|1982|p=265, 269-270}}</ref>


Two integer arrays u(3) and v(3) are initialized to the lowest terms of the Perrin sequence, with positive indices {{math|t {{=}} 0, 1, 2}} in u( ) and negative indices {{math|t {{=}} 0,−1,−2}} in v( ).
Two integer arrays u(3) and v(3) are initialized to the lowest terms of the Perrin sequence, with positive indices {{math|t {{=}} 0, 1, 2}} in u( ) and negative indices {{math|t {{=}} 0,−1,−2}} in v( ).
Line 211: Line 211:
Successively {{math|P(−n − 1), P(−n), P(−n + 1)}} and {{math|P(n − 1), P(n), P(n + 1) (mod n)}}.
Successively {{math|P(−n − 1), P(−n), P(−n + 1)}} and {{math|P(n − 1), P(n), P(n + 1) (mod n)}}.


If {{math|P(−n) {{=}} −1}} and {{math|P(n) {{=}} 0}} then ''n'' is a [[probable prime]], that is: actually prime or a strong Perrin pseudoprime.
If {{math|P(−n) {{=}} −1}} and {{math|P(n) {{=}} 0}} then ''n'' is a [[probable prime]], that is: actually prime or a restricted Perrin pseudoprime.


Shanks ''et al.'' observed that for all strong pseudoprimes they found, the final state of the above six registers (the "signature" of ''n'') equals the initial state 1,−1,3, 3,0,2.<ref>{{harvtxt|Adams|Shanks|1982|p=275}}, {{harvtxt|Kurtz|Shanks|Williams|1986|p=694}}. This was later confirmed for {{math|n < 10{{sup|14}}}} by Steven {{harvtxt|Arno|1991}}.</ref> The same occurs with {{math|≈ 1'''/'''6}} of all primes, so the two sets cannot be distinguished on the strength of this test alone.<ref>The signature does give discriminating information on the remaining two types of primes.
Shanks ''et al.'' observed that for all restricted pseudoprimes they found, the final state of the above six registers (the "signature" of ''n'') equals the initial state 1,−1,3, 3,0,2.<ref>{{harvtxt|Adams|Shanks|1982|p=275}}, {{harvtxt|Kurtz|Shanks|Williams|1986|p=694}}. This was later confirmed for {{math|n < 10{{sup|14}}}} by Steven {{harvtxt|Arno|1991}}.</ref> The same occurs with {{math|≈ 1'''/'''6}} of all primes, so the two sets cannot be distinguished on the strength of this test alone.<ref>The signature does give discriminating information on the remaining two types of primes.
For example, the smallest Q-type pseudoprime 50972694899204437633 computed by Holger {{harvtxt|Stephan|2019}} is exposed by signature conditions 14a and 14c in {{harvtxt|Adams|Shanks|1982|p=257}}.</ref> For those cases, they recommend to also use the [[Supergolden_ratio#Narayana sequence|Narayana-Lucas sister sequence]] with recurrence relation {{math|A(n) {{=}} A(n − 1) + A(n − 3)}} and initial values
For example, the smallest Q-type pseudoprime 50,972,694,899,204,437,633 computed by Holger {{harvtxt|Stephan|2019}} is exposed by signature conditions 14a and 14c in {{harvtxt|Adams|Shanks|1982|p=257}}.</ref> For those cases, they recommend to also use the [[Supergolden_ratio#Narayana sequence|Narayana-Lucas sister sequence]] with recurrence relation {{math|A(n) {{=}} A(n − 1) + A(n − 3)}} and initial values


u(0):= 3, u(1):= 1, u(2):= 1
u(0):= 3, u(1):= 1, u(2):= 1
Line 231: Line 231:
Here, ''n'' is a probable prime if {{math|A(−n) {{=}} 0}} and {{math|A(n) {{=}} 1}}.
Here, ''n'' is a probable prime if {{math|A(−n) {{=}} 0}} and {{math|A(n) {{=}} 1}}.


Kurtz ''et al.'' found no overlap between the odd pseudoprimes for the two sequences below {{math|50∙10{{sup|9}}}} and supposed that 2277740968903 = 1067179 ∙ 2134357 is the smallest composite number to pass both tests.<ref>{{harvtxt|Kurtz|Shanks|Williams|1986|p=697}}</ref>
Kurtz ''et al.'' found no overlap between the odd pseudoprimes for the two sequences below {{math|50∙10{{sup|9}}}} and supposed that 2,277,740,968,903 = 1067179 ∙ 2134357 is the smallest composite number to pass both tests.<ref>{{harvtxt|Kurtz|Shanks|Williams|1986|p=697}}</ref>

If the third-order Pell-Lucas recurrence {{math|A(n) {{=}} 2A(n − 1) + A(n − 3)}} is next tried,<ref>{{OEIS|A332647}}</ref> this bound will be pushed up to 4,057,052,731,496,380,171 = 1424263447 ∙ 2848526893.

Additionally, roots of the doubling rule-congruence <math>x^2 - 2x \equiv P(0) \equiv 3 \pmod{n}\,</math> other than {{math|−1}} or {{math|3}} expose composite numbers, like non-trivial square roots of {{math|1}} in the [[Miller–Rabin primality test|Miller–Rabin test]].<ref>{{harvtxt|Adams|Shanks|1982|p=280-283}}</ref> This reduces the number of pseudoprimes for each sequence by approximately one-third and is especially efficient in detecting [[Carmichael number]]s.

===C/C++ program===

If the cubic primality test is implemented using 64-bit [[primitive data type]]s, combining the above steps will leave no known pseudoprimes in its domain.

<syntaxhighlight lang="cpp">
/************************************************
Adams & Shanks 1982 Strong cubic primality test,
max. 63-bit argument. C99/C++
-----------------------------------------------*/
#include <inttypes.h>
#include <stdint.h>
#include <stdio.h>

//initial register values for
//x^3 = x + 1, x^3 = x^2 + 1 and x^3 = 2x^2 + 1.
const int64_t reg[3][6] = {
{1,-1,3, 3,0,2},
{-2,0,3, 3,1,1},
{-4,0,3, 3,2,4} };

//mask, leading zeroes, n to test
uint64_t msk, lz, n;


//a + b (mod n) with a + b < 2n
#define addm(r, a, b) r = a + b; if (r >= n) r -= n
//positive difference a - b (mod n)
#define subm(r, a, b) r = n - b + a; if (r >= n) r -= n

//double the index of u
uint64_t dblidx (uint64_t u, uint64_t v)
{
uint64_t r = 0;
//modular squaring without overflow
for (uint64_t a = u; a > 0; a >>= lz) {
r = (r + (a & msk) * u) % n;
u = (u << lz) % n;
}
//positive difference u^2 - 2v
addm(u, v, v);
subm(r, r, u);
return r;
}

//fill the gaps -/+ 1
void fill (uint64_t *a, int i, int j, int seq)
{
uint64_t u, v;
//copy to end of array
a[3] = a[2];

switch (seq) {
case 0:
subm(u, a[j], a[1]);
addm(a[j], a[i], u);
a[i] = u; break;
case 1:
addm(u, a[i], a[1]);
subm(a[i], a[j], u);
a[j] = u; break;
default:
addm(v, a[1], a[1]);
addm(u, a[i], v);
addm(v, u, u);
subm(a[i], a[j], v);
a[j] = u; break;
}
}

//centre and print signature triple (mod n)
#define triple(a) { \
for (int i = 0; i < 3; i++) { \
if (a > n - a) a -= n; \
printf(" % " PRId64, a); \
} printf("\n"); }

/* evaluate signature
return 0 for composite n
1 for probably prime n
3 the test is inconclusive,
repeat with the next sequence. */
int eval (int64_t u[], int64_t v[], int seq, int fl)
{
//result
triple(v[2 - i]);
triple(u[i]);

//exceptional cases 2, 3
if (n < 4) return 1;

//necessary for primes
int sw =(v[1]==reg[seq][1]) && (u[1]==reg[seq][4]);
if (!sw) return 0;

if (v[0]==u[0]) {
//repeat for S-type; exposed if fl==0
if (u[0]==3) sw = (fl) ? 3 : 0; }
else {
uint64_t addm(d, u[0], v[0] + 3);
sw &= d==0; }
//sufficient for Q and I-types below 2^63
return sw;
}

//cubic primality test
int Cubtest (int seq)
{
uint64_t t = 1, n1 = n - 1;
int64_t u[4], v[4];
int sp;

//trailing zeroes n - 1
while ((n1 & t)==0) t <<= 1;

//initialize registers
for (int i = 0; i < 3; i++) {
u[i] = reg[seq][i + 3];
v[i] = reg[seq][2 - i];
while (v[i] < 0) v[i] += n;
//exceptional cases
if (n < 5) { u[i] %= n; v[i] %= n; }
}

//double and add
for (uint64_t k = 1uLL << (62 - lz); \
k > 0; k >>= 1) {

//strong psp property (2)
if (k < t)
sp = sp || ((v[1]==n1) && (u[1]==n1));

//double
for (int i = 0; i < 3; i++) {
uint64_t d = dblidx(u[i], v[i]);
v[i] = dblidx(v[i], u[i]);
u[i] = d;
}

//rein in
fill((uint64_t *) u, 0, 2, seq);
fill((uint64_t *) v, 2, 0, seq);

if (n & k) {
//index + 1
for (int i = 0; i < 3; i++) {
u[i] = u[i + 1];
v[i] = v[i + 1];
}

//strong psp property (1)
if (k==t)
sp =(v[1]==3) && (u[1]==3);
}
}
return eval(u, v, seq, sp);
}

int main(void)
{
while (1) {
scanf("%" SCNu64, &n);
if (n==0) break;
printf("\nN = %" PRIu64 "\n", n);

//count leading zeroes
lz = 0;
for (msk = INT64_MIN; \
(n & msk)==0; msk >>= 1) lz++;

//max. 63-bit argument
if ((lz==0) || (lz==63)) continue;
msk = (1uLL << lz) - 1;

int sw =-1;
for (int i = 0; i < 3; i++) {
sw &= Cubtest(i);
if (sw < 2) break;
}

if (sw & 1)
printf("probably prime\n");
else
printf("composite\n");
}
return 0;
}</syntaxhighlight>


==Notes==
==Notes==

Revision as of 16:17, 19 April 2024

Spiral of equilateral triangles with side lengths equal to Perrin numbers.

In mathematics, the Perrin numbers are a doubly infinite constant-recursive integer sequence with characteristic equation x3 = x + 1. The Perrin numbers bear the same relationship to the Padovan sequence as the Lucas numbers do to the Fibonacci sequence.

Definition

The Perrin numbers are defined by the recurrence relation

and the reverse

The first few terms in both directions are

n 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
P(n) 3 0 2 3 2 5 5 7 10 12 17 22 29 39 51 68 90 119 ...[1]
P(−n) ... −1 1 2 −3 4 −2 −1 5 −7 6 −1 −6 12 −13 7 5 −18 ...[2]

Perrin numbers can be expressed as sums of the three initial terms

The first fourteen prime Perrin numbers are

n 2 3 4 5 6 7 10 12 20 21 24 34 38 75 ...[3]
P(n) 2 3 2 5 5 7 17 29 277 367 853 14197 43721 1442968193 ...[4]

History

In 1876 the sequence and its equation were initially mentioned by Édouard Lucas, who noted that the index n divides term P(n) if n is prime.[5] In 1899 Raoul Perrin asked if there were any counterexamples to this property.[6] The first P(n) divisible by composite index n was found only in 1982 by William Adams and Daniel Shanks.[7] They presented a detailed investigation of the sequence, with a sequel appearing four years later.[8]

Properties

A Perrin microcosm: the escape time algorithm is applied to the Newton map of the entire Perrin function F(z) around critical point z = 1.225432 with viewport width 0.05320. The basins of attraction emanating from the centre correspond to the infinite number of real (left) and complex roots (right) F(z) = 0.

The Perrin sequence also satisfies the recurrence relation Starting from this and the defining recurrence, one can create an infinite number of further relations, for example

The generating function of the Perrin sequence is

The sequence is related to sums of binomial coefficients by

[9]

Perrin numbers can be expressed in terms of partial sums

The Perrin numbers are obtained as integral powers n ≥ 0 of the matrix

and its inverse

The Perrin analogue of the Simson identity for Fibonacci numbers is given by the determinant

The number of different maximal independent sets in an n-vertex cycle graph is counted by the nth Perrin number for n ≥ 2.[10]

Binet formula

The Perrin function extends the sequence to real numbers.

The solution of the recurrence can be written in terms of the roots of characteristic equation . If the three solutions are real root (with approximate value 1.324718 and known as the plastic ratio) and complex conjugate roots and , the Perrin numbers can be computed with the Binet formula which also holds for negative n.

The polar form is with Since the formula reduces to either the first or the second term successively for large positive or negative n, and numbers with negative subscripts oscillate. Provided α is computed to sufficient precision, these formulas can be used to calculate Perrin numbers for large n.

Expanding the identity gives the important index-doubling rule by which the forward and reverse parts of the sequence are linked.

Prime index p divides P(p)

If the characteristic equation of the sequence is written as then the coefficients can be expressed in terms of roots with Vieta's formulas:

These integer valued functions are the elementary symmetric polynomials in

  • The fundamental theorem on symmetric polynomials states that every symmetric polynomial in the complex roots of monic can be represented as another polynomial function in the integer coefficients of
  • The analogue of Lucas's theorem for multinomial coefficients says that if then is divisible by prime

Given integers a, b, c and n > 0,

Rearrange into symmetric monomial summands, permuting exponents i, j, k:

Substitute prime for power and complex roots for integers and compute representations in terms of for all symmetric polynomial functions. For example, is and the power sum can be expressed in the coefficients with Newton's recursive scheme. It follows that the identity has integer terms and both sides divisible by prime

To show that prime divides in the reverse sequence, the characteristic equation has to be reflected. The roots are then the coefficients and the same reasoning applies.

Perrin primality test

Query 1484. The curious proposition of Chinese origin which is the subject of query 1401[11] would provide, if it is true, a more practical criterium than Wilson's theorem for verifying whether a given number m is prime or not; it would suffice to calculate the residues with respect to m of successive terms of the recurrence sequence
un = 3un−1 − 2un−2 with initial values u0 = −1, u1 = 0.[12]
I have found another recurrence sequence that seems to possess the same property; it is the one whose general term is
vn = vn−2 + vn−3 with initial values v0 = 3, v1 = 0, v2 = 2. It is easy to demonstrate that vn is divisible by n, if n is prime; I have verified, up to fairly high values of n, that in the opposite case it is not; but it would be interesting to know if this is really so, especially since the sequence vn gives much less rapidly increasing numbers than the sequence un (for n = 17, for example, one finds un = 131070, vn = 119), which leads to simpler calculations when n is a large number.
The same proof, applicable to one of the sequences, will undoubtedly bear upon the other, if the stated property is true for both: it is only a matter of discovering it.[13]

The Perrin sequence has the Fermat property: if p is prime, P(p) ≡ P(1) ≡ 0 (mod p). However, the converse is not true: some composite n may still divide P(n). A number with this property is called a Perrin pseudoprime.

The question of the existence of Perrin pseudoprimes was considered by Malo and Jarden,[14] but none were known until Adams and Shanks found the smallest one, 271441 = 5212 (the number P(271441) has 33150 decimal digits).[15] Jon Grantham later proved that there are infinitely many Perrin pseudoprimes.[16]

The seventeen Perrin pseudoprimes below 109 are

271441, 904631, 16532714, 24658561, 27422714, 27664033, 46672291, 102690901, 130944133, 196075949, 214038533, 517697641, 545670533, 801123451, 855073301, 903136901, 970355431.[17]

Adams and Shanks noted that primes also satisfy the congruence P(−p) ≡ P(−1) ≡ −1 (mod p). Composites for which both properties hold are called restricted Perrin pseudoprimes. There are only nine such numbers below 109.[18]

While Perrin pseudoprimes are rare, they overlap with Fermat pseudoprimes. Of the above seventeen numbers, four are base 2 Fermatians as well. In contrast, the Lucas pseudoprimes are anti-correlated.[19] Presumably, combining the Perrin and Lucas tests should make a primality test as strong as the reliable BPSW test which has no known pseudoprimes – though at higher computational cost.

Pseudocode

The 1982 Adams and Shanks O(log n) Perrin primality test.[20]

Two integer arrays u(3) and v(3) are initialized to the lowest terms of the Perrin sequence, with positive indices t = 0, 1, 2 in u( ) and negative indices t = 0,−1,−2 in v( ).

The main double-and-add loop, originally devised to run on an HP-41C pocket calculator, computes P(n) mod n and the reverse P(−n) mod n at the cost of six modular squarings for each bit of n.

The subscripts of the Perrin numbers are doubled using the identity P(2t) = P2(t) − 2P(−t). The resulting gaps between P(±2t) and P(±2t ± 2) are closed by applying the defining relation P(t) = P(t − 2) + P(t − 3).

Initial values
let int u(0):= 3, u(1):= 0, u(2):= 2
let int v(0):= 3, v(1):=−1, v(2):= 1
     
Test odd positive number n
input int n
     
set int h:= most significant bit of n
     
for k:= h − 1 downto 0
     
  Double the indices of
  the six Perrin numbers.
  for i = 0, 1, 2
    temp:= u(i)^2 − 2v(i) (mod n)
    v(i):= v(i)^2 − 2u(i) (mod n)
    u(i):= temp
  endfor
     
  Copy P(2t + 2) and P(−2t − 2)
  to the array ends and use
  in the if-statement below.
  u(3):= u(2)
  v(3):= v(2)
     
  Overwrite P(2t ± 2) with P(2t ± 1)
  temp:= u(2) − u(1)
  u(2):= u(0) + temp
  u(0):= temp
     
  Overwrite P(−2t ± 2) with P(−2t ± 1)
  temp:= v(0) − v(1)
  v(0):= v(2) + temp
  v(2):= temp
     
  if n has bit k set then
    Increase the indices of
    both Perrin triples by 1.
    for i = 0, 1, 2
      u(i):= u(i + 1)
      v(i):= v(i + 1)
    endfor
  endif
     
endfor
     
Result
print v(2), v(1), v(0)
print u(0), u(1), u(2)

Successively P(−n − 1), P(−n), P(−n + 1) and P(n − 1), P(n), P(n + 1) (mod n).

If P(−n) = −1 and P(n) = 0 then n is a probable prime, that is: actually prime or a restricted Perrin pseudoprime.

Shanks et al. observed that for all restricted pseudoprimes they found, the final state of the above six registers (the "signature" of n) equals the initial state 1,−1,3, 3,0,2.[21] The same occurs with ≈ 1/6 of all primes, so the two sets cannot be distinguished on the strength of this test alone.[22] For those cases, they recommend to also use the Narayana-Lucas sister sequence with recurrence relation A(n) = A(n − 1) + A(n − 3) and initial values

u(0):= 3, u(1):= 1, u(2):= 1
v(0):= 3, v(1):= 0, v(2):=−2

The same doubling rule applies and the formulas for filling the gaps are

temp:= u(0) + u(1)
u(0):= u(2) − temp
u(2):= temp
     
temp:= v(2) + v(1)
v(2):= v(0) − temp
v(0):= temp

Here, n is a probable prime if A(−n) = 0 and A(n) = 1.

Kurtz et al. found no overlap between the odd pseudoprimes for the two sequences below 50∙109 and supposed that 2,277,740,968,903 = 1067179 ∙ 2134357 is the smallest composite number to pass both tests.[23]

If the third-order Pell-Lucas recurrence A(n) = 2A(n − 1) + A(n − 3) is next tried,[24] this bound will be pushed up to 4,057,052,731,496,380,171 = 1424263447 ∙ 2848526893.

Additionally, roots of the doubling rule-congruence other than −1 or 3 expose composite numbers, like non-trivial square roots of 1 in the Miller–Rabin test.[25] This reduces the number of pseudoprimes for each sequence by approximately one-third and is especially efficient in detecting Carmichael numbers.

C/C++ program

If the cubic primality test is implemented using 64-bit primitive data types, combining the above steps will leave no known pseudoprimes in its domain.

/************************************************
Adams & Shanks 1982 Strong cubic primality test,
max. 63-bit argument. C99/C++
-----------------------------------------------*/
#include <inttypes.h>
#include <stdint.h>
#include <stdio.h>

//initial register values for
//x^3 = x + 1, x^3 = x^2 + 1 and x^3 = 2x^2 + 1.
const int64_t reg[3][6] = {
  {1,-1,3, 3,0,2},
  {-2,0,3, 3,1,1},
  {-4,0,3, 3,2,4} };

//mask, leading zeroes, n to test
uint64_t msk, lz, n;


//a + b (mod n) with a + b < 2n
#define addm(r, a, b) r = a + b; if (r >= n) r -= n
//positive difference a - b (mod n)
#define subm(r, a, b) r = n - b + a; if (r >= n) r -= n

//double the index of u
uint64_t dblidx (uint64_t u, uint64_t v)
{
uint64_t r = 0;
   //modular squaring without overflow
   for (uint64_t a = u; a > 0; a >>= lz) {
      r = (r + (a & msk) * u) % n;
      u = (u << lz) % n;
   }
   //positive difference u^2 - 2v
   addm(u, v, v);
   subm(r, r, u);
return r;
}

//fill the gaps -/+ 1
void fill (uint64_t *a, int i, int j, int seq)
{
uint64_t u, v;
   //copy to end of array
   a[3] = a[2];

   switch (seq) {
   case 0:
      subm(u, a[j], a[1]);
      addm(a[j], a[i], u);
      a[i] = u; break;
   case 1:
      addm(u, a[i], a[1]);
      subm(a[i], a[j], u);
      a[j] = u; break;
   default:
      addm(v, a[1], a[1]);
      addm(u, a[i], v);
      addm(v, u, u);
      subm(a[i], a[j], v);
      a[j] = u; break;
   }
}

//centre and print signature triple (mod n)
#define triple(a) {              \
   for (int i = 0; i < 3; i++) { \
      if (a > n - a) a -= n;     \
      printf(" % " PRId64, a);   \
   } printf("\n"); }

/* evaluate signature
return 0 for composite n
       1 for probably prime n
       3 the test is inconclusive,
         repeat with the next sequence. */
int eval (int64_t u[], int64_t v[], int seq, int fl)
{
   //result
   triple(v[2 - i]);
   triple(u[i]);

//exceptional cases 2, 3
if (n < 4) return 1;

//necessary for primes
int sw =(v[1]==reg[seq][1]) && (u[1]==reg[seq][4]);
if (!sw) return 0;

   if (v[0]==u[0]) {
      //repeat for S-type; exposed if fl==0
      if (u[0]==3) sw = (fl) ? 3 : 0; }
   else {
uint64_t addm(d, u[0], v[0] + 3);
      sw &= d==0; }
   //sufficient for Q and I-types below 2^63
return sw;
}

//cubic primality test
int Cubtest (int seq)
{
uint64_t t = 1, n1 = n - 1;
int64_t u[4], v[4];
int sp;

   //trailing zeroes n - 1
   while ((n1 & t)==0) t <<= 1;

   //initialize registers
   for (int i = 0; i < 3; i++) {
      u[i] = reg[seq][i + 3];
      v[i] = reg[seq][2 - i];
      while (v[i] < 0) v[i] += n;
      //exceptional cases
      if (n < 5) { u[i] %= n; v[i] %= n; }
   }

   //double and add
   for (uint64_t k = 1uLL << (62 - lz); \
      k > 0; k >>= 1) {

      //strong psp property (2)
      if (k < t)
         sp = sp || ((v[1]==n1) && (u[1]==n1));

      //double
      for (int i = 0; i < 3; i++) {
   uint64_t d = dblidx(u[i], v[i]);
         v[i] = dblidx(v[i], u[i]);
         u[i] = d;
      }

      //rein in
      fill((uint64_t *) u, 0, 2, seq);
      fill((uint64_t *) v, 2, 0, seq);

      if (n & k) {
         //index + 1
         for (int i = 0; i < 3; i++) {
            u[i] = u[i + 1];
            v[i] = v[i + 1];
         }

         //strong psp property (1)
         if (k==t)
            sp =(v[1]==3) && (u[1]==3);
      }
   }
return eval(u, v, seq, sp);
}

int main(void)
{
   while (1) {
      scanf("%" SCNu64, &n);
      if (n==0) break;
      printf("\nN = %" PRIu64 "\n", n);

      //count leading zeroes
      lz = 0;
      for (msk = INT64_MIN; \
        (n & msk)==0; msk >>= 1) lz++;

      //max. 63-bit argument
      if ((lz==0) || (lz==63)) continue;
      msk = (1uLL << lz) - 1;

      int sw =-1;
      for (int i = 0; i < 3; i++) {
         sw &= Cubtest(i);
         if (sw < 2) break;
      }

      if (sw & 1)
         printf("probably prime\n");
      else
         printf("composite\n");
   }
return 0;
}

Notes

  1. ^ Sloane, N. J. A. (ed.). "Sequence A001608". The On-Line Encyclopedia of Integer Sequences. OEIS Foundation.
  2. ^ (sequence A078712 in the OEIS)
  3. ^ (sequence A112881 in the OEIS)
  4. ^ (sequence A074788 in the OEIS)
  5. ^ Lucas (1878)
  6. ^ Perrin (1899)
  7. ^ Adams & Shanks (1982)
  8. ^ Kurtz, Shanks & Williams (1986)
  9. ^ (sequence A001608 in the OEIS)
  10. ^ Füredi (1987)
  11. ^ Tarry (1898)
  12. ^ equivalently un = 2n − 2. (sequence A000918 in the OEIS)
  13. ^ Perrin (1899) translated from the French
  14. ^ Malo (1900), Jarden (1966)
  15. ^ Adams & Shanks (1982, p. 255)
  16. ^ Grantham (2010), Stephan (2020)
  17. ^ (sequence A013998 in the OEIS)
  18. ^ (sequence A018187 in the OEIS), (sequence A275612 in the OEIS), (sequence A275613 in the OEIS)
  19. ^ None of the 2402549 Lucas-Selfridge pseudoprimes below 1015 listed by Dana Jacobsen (2020) is also a Perrin pseudoprime.
  20. ^ Adams & Shanks (1982, p. 265, 269-270)
  21. ^ Adams & Shanks (1982, p. 275), Kurtz, Shanks & Williams (1986, p. 694). This was later confirmed for n < 1014 by Steven Arno (1991).
  22. ^ The signature does give discriminating information on the remaining two types of primes. For example, the smallest Q-type pseudoprime 50,972,694,899,204,437,633 computed by Holger Stephan (2019) is exposed by signature conditions 14a and 14c in Adams & Shanks (1982, p. 257).
  23. ^ Kurtz, Shanks & Williams (1986, p. 697)
  24. ^ (sequence A332647 in the OEIS)
  25. ^ Adams & Shanks (1982, p. 280-283)

References