Find all Perrin pseudoprimes less than 100 million

Question

Simply put, the Perrin sequence has the following recurrence relation:

P(n) = P(n-2) + P(n-3);

P(0) = 3; P(1) = 0; P(2) = 2;

But wait!! There's more!!

The Perrin sequence has a very bizarre property:

If N is prime, then P(N) mod N = 0.

Furthermore, there are very few composite numbers which have this property. This can serve as a simple primality test which gives very few false positives.

There are 5761462 numbers under 100 million which satisfy this congruence. Out of those, 5761455 (99.9999%) are prime numbers.

Your mission is to compute those 7 numbers which pass the test but which are not prime. To be specific, your program must output these numbers in a reasonable format:

271441, 904631, 16532714, 24658561, 27422714, 27664033, 46672291

_{You can find more perrin pseudo primes (all of them up to 10^16) with factorization here. For this challenge only the first 7 are important, since they are the only one smaller than 100 million.}

(Your program must compute these numbers without having them pre-made, or accessing the internet, or similar stuff.)

Also, your program should not produce integer overflows. This may at first seem difficult because the 100 millionth Perrin number is really, really big. Algorithms do exist, however, which would not cause overflows.

This is code golf, shorter solutions are better.

Answer 1

Ruby, 183

require'Prime'
f=->a{a=a.flat_map{|n|String===n ? n : n<3?%w[3 0 2][n]:[n-2,n-3]} until a.all?{|x|String===x}
a.map(&:to_i).inject :+}
n=1
p n if f[[n+=1]]%n==0&&!n.prime? while n<1e8

Technically works but runs in an unreasonable amount of time. Try setting the lower bound to 271435 and upper bound to 271450.

Answer 2

PARI: 148

The oeis webpage on the Perrin number contains a PARI code written by Joerg Arndt that can easily be modified (a number & spacing) to fit the criteria

N=10^8; 
default(primelimit, N); 
M=[0,1,0;0,0,1;1,1,0]; 
a(n)=lift(trace(Mod(M,n)^n));
{for(n=1,N,if(isprime(n),next());if(a(n)==0,print1(n,",");););}

It takes quite a while, but it definitely works. I will work on translating this to Fortran for speed computation, but it definitely will not be small (especially if I put the prime test in the code & not as a use <module> declaration).

Answer 3

GolfScript, 53 50 characters

10 8?,{3 0 2{@2$+3$%}4$2-*])!*},{.,2>{1$\%!},*},n*

Do not try to run this version - it'll take very long time to complete (if ever). The code creates a list of all candidates and then consists of two filter steps, first to select all assumed primes and then filter all non-primes from the remaining list. In order to see the code work you can try the following snippet:

30,{271430+}%
.p
{3 0 2{@2$+3$%}4$2-*])!*},
.p
{.,2>{1$\%!},*},
p

Answer 4

Python (127)

a=i=3
b,c=0,2
while i<10**8:
    c,b,a=a+b,c,b
    j=2
    if c%i==0:
        while j*j<=i:
            if i%j==0:print i;break
            j+=1
    i+=1

Unnecessary usage of memory avoided but the solution tends to be slow due to big integer operations. Primality testing is done with a slow algo.

Answer 5

AutoHotkey 132

o:={-2:1,-1:-1,0:3}
loop % 10**8
{
p:=0
if !Mod(o[(l:=A_index)-2]+o[l-3],l)
loop % l
if !Mod(l,A_index)
p++
if p>2
s.=l ","
}
msgbox % s  

Not possible as Integer limit is surpassed.

Answer 6

Husk, 22 bytes

↑7fo¬ṗtW¬z%Ntƒ(+d302Ẋ+

(Don't bother to) Try it online!: as many other answers this is extremely slow, so you won't see anything in the 60s TIO allows. Try this online instead! for a modified version that finds the first Perrin pseudoprime after 271440 (yes, it's 271441).

Explanation

↑7fo¬ṗtW¬z%Ntƒ(+d302Ẋ+
             ƒ(+d302Ẋ+    Recursive definition of the Perrin sequence:
                d302       Start with the digits of 302
               +           and append after them  
                    Ẋ+     the pairwise sums
             ƒ             of the sequence itself
            t             Drop the first value (we like counting from 1 in Husk)
         z%N              Compute each value modulo its index
       W¬                 Find the indices of 0s
      t                   Drop the fist one (0 mod 1 is 0, but for some reason
                           1 is not one of the required outputs)
  f                       Keep only those values
   o¬                     that are not
     ṗ                    primes
↑7                        Take the first 7 results

Answer 7

Haskell (169)

p=3:0:2:zipWith(+)p(tail p)
s(x:y)=x:s[a|a<-y,a`rem`x/=0]
s[]=[]
w=s$2:[3,5..10^8]
k(y,x)z|x`rem`y==0&&(y`notElem`w)=y:z|1<2=z
main=print.foldr k[]$zip[2..10^8].drop 2$p

p is an infite list producing the Perrin numbers, s is a simple sieve function, w contains all primes up to 10^8. k checks whether the pair (y,x) (where y is the Perrin number's index and x is the number itself) fulfils the condition and adds them to the list z. main folds k along the zip of the list [2..10^8] and the Perrin numbers (except for N=0 and N=1).

Don't run this program. It will generate the prime numbers first, which takes very, very, very long due to the inefficient sieve. Faster versions would use a wheel sieve and a binary search instead of the current linear search. Also, starting at around 10^5, garbage collection will eat up almost all productivity.

Answer 8

Perl 68

use ntheory":all";forcomposites{say if is_perrin_pseudoprime($_)}1e8

Takes 66 seconds on my macbook with the latest version, 6 minutes with previous one (that had no prefilters -- basically a sped up version of Arndt's Pari code).

Of course this is cheesy by using a module that does all the work in C. 88 lines of C with prefilters for faster operation, 28 lines if removing the filters and an optimization to avoid mulmod for small inputs.

Answer 9

Jelly, 30 bytes

Ø.’ịṖS;@
Ø.Ḥ3;Ç⁸¡ị@‘
Ẓ=Ç%$Ɗ8#Ḋ

Try it online!

Wildly inefficient, times out on TIO before printing anything.

How it works

Ø.’ịṖS;@    - Helper link. Takes a list l on the left
            - l is the list of previous Perrin numbers and the link appends the next Perrin number
Ø.’         - Yield [-1, 0]
    Ṗ       - Remove the last element of l
   ị        - Index into the cropped version of l
            - This gives the second and third elements from the end of l
     S      - Take the sum of these two
      ;@    - Append to l

Ø.Ḥ3;Ç⁸¡ị@‘ - Helper link. Takes an integer k on the left
            - This link takes an integer k and returns the kth Perrin number
Ø.          - Yield [0, 1]
  Ḥ         - Double to [0, 2]
   3;       - Prepend 3
      ⁸¡    - Do the following k times, operating on the list [3, 0 2]:
     Ç      -   Call the above helper link
          ‘ - Yield k+1
        ị@  - Take the 1-indexed (k+1)th index of this list

Ẓ=Ç%$Ɗ8#Ḋ   - Main link. Takes no arguments
     Ɗ8#    - Do the following over n = 0, 1, 2, ... until 8 values return True:
Ẓ           -   Is n prime?
    $       -   Group the previous two links and run over n:
  Ç         -     Get the nth Perrin number
   %        -     Modulo it by n
 =          -   Are the two values equal?
            -   The two values will either be equal if:
            -     n is prime and P(n) mod n = 1 or
            -     n is composite and P(n) mod n = 0 or
            -     n = 1
            -   The first is never true, therefore this only returns True if n is composite and P(n) mod n = 0 or n = 1
            - This returns [1, 271441, 904631, 16532714, 24658561, 27422714, 27664033, 46672291]
        Ḋ   - Remove 1 from this list