# Find all Perrin pseudoprimes less than 100 million

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.

• What is P(2)...? Dec 30, 2013 at 19:04
• @DoorknobofSnow Simple mistake. Dec 30, 2013 at 19:06
• Is having p(0),p(1), & p(2) stored as a single global array okay? Dec 30, 2013 at 19:24
• Do you think this would have made a better fastest-code competition? Dec 30, 2013 at 19:48
• Do you want us to consider up to the 100 millionth Perrin number, or numbers under 100 million? Dec 30, 2013 at 22:02

# 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.

## 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).

• With newer Pari/GP this can be reduced. Get rid of the first two lines, change the last to forcomposite(n=1,10^8,a(n)||print1(n",")) Jun 5, 2015 at 19:06

### 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

• Can be shortened a bit if both tests are packed into one but then it'll be even slower. Dec 30, 2013 at 19:54
• Unless I've missed something, the +,( could just be ,*. Dec 30, 2013 at 23:11
• And similarly, I think -1=! could be )!*. Dec 30, 2013 at 23:30
• @PeterTaylor Of course you're right. I was a little bit too anxious to leave anything on the stack and forgot that you can use the array as boolean also. Dec 31, 2013 at 8:20

## 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.

• You can use 1e8 instead of 10**8, I think. May 8, 2015 at 14:04

## 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.

• According to the AutoHotkey documentation, integers are signed 64-bit. The largest number required for a direct approach like this is about 40000000 bits. Are you sure this doesn't overflow and give nonsense results? Jan 2, 2014 at 8:41
• @PeterTaylor According to doc "AutoHotkey supports 64-bit signed integers, which range from -9223372036854775808 (-0x8000000000000000) to 9223372036854775807 (0x7FFFFFFFFFFFFFFF)." 10**8 is less than 9223372036854775807 and that is the reason I think it's okay.
– Avi
Jan 2, 2014 at 11:49
• 108 is fine, but you're computing 1.3...**(108), which is much much larger. Jan 2, 2014 at 14:11
• @PeterTaylor I didn't thought P(n) for say n=100 can get that big. Yes, you are right. Ahk won't compute that. I will edit that.
– Avi
Jan 3, 2014 at 3:40

# 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:
+           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


p=3:0:2:zipWith(+)p(tail p)
s(x:y)=x:s[a|a<-y,aremx/=0]
s[]=[]
w=s$2:[3,5..10^8] k(y,x)z|xremy==0&&(ynotElemw)=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. # 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.

# 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