Finding Perrin Pseudo Primes_Part 1 


My latest posts can be found here: Previous blog posts:
Additionally, some earlier writings: 
Finding Perrin Pseudoprimes:

This isn't an exact quotation, but it's close enough. 
counterexample. 
So I started on that. I figured Andrew had already done the basic work of looking for a small counterexample, but I also figured that I would need to get warmed up and start to get a feel for the problem.
Before that, I need to put this into some historical context. This was when, if you wanted to do arithmetic on more than 32 bit numbers, you had to write your own arbitrary arithmetic package. It was before the WWW as we know it, although not preinternet per se. So some of what follows will seem unnecessary. If you think that, I invite you to find what is common for the first 1000 counterexamples, except for a small number of them, and tell me which ones do not follow the pattern.
And so ...

We can work out the ratio between successive terms and see that it starts to converge quite rapidly to about 1.324736... In fact we can show that:
$$k(n)=\alpha^n+\beta^n+\gamma^n$$
where $\alpha,$ $\beta,$ and $\gamma$ are solutions to the equation:
$$x^3=x+1$$
(Usually we need constants in front of the powers of the roots, but in this case the constants all turn out to be 1).
No? Can't see that? Well the definition is: $$k(n)=a.k(n1)+b.k(n2)+c.k(n3)$$ where a=0, b=1, and c=1. Then the characteristic equation is $x^3=a.x^2+b.x+c.$ Anyway ... 
The roots to that cubic are:
Now the number of digits in a number m (base 10) is given by log(m) (base 10), so that means the number of digits in k(n) is log(k(n)) which is roughly $\log(\alpha^n),$ which in turn is $n.\log(\alpha).$
Now $\log(\alpha)$ is about 0.122..., so by the time n is 100, k(n) will have 12 digits (actually 13, because we should round up), and it will only get worse.
So we can't go far with this by hand. But that's OK, because what we want to know is this:
the Perrin Sequence, so I call this the "Perrin Test." 
#!/usr/bin/python from time import clock for limit_power in [0, 1, 2, 3, 4]: n = 2 limit = 4**limit_power t_start = clock() while clock()t_start < limit: a,b,c = 3,0,2 for i in xrange(n): a,b,c = b,c,(a+b) % n n += 1 print limit, n, float(n)/2**limit_power 
So here we have some Python code to compute k(n) mod n for as long as allowed by some time limit, and then we print how far we got. We're not actually doing any of the searching for a counterexample in this code, since we're not actually seeing whether the numbers are predicted to be prime, nor are we seeing if they actually are prime. This is just showing that it's going to be a slow process if we use this method.

#!/usr/bin/python from time import clock from sys import argv def is_prime_q(n): if n in [2, 3, 5, 7]: return True if n < 11 : return False if n % 2 == 0 : return False if pow(2,n1,n) != 1: return False if n % 3 == 0 : return False if pow(3,n1,n) != 1: return False if n % 5 == 0 : return False if pow(5,n1,n) != 1: return False d = 7 while d*d <= n: if (n % d) == 0: return False d += 2 return True def is_perrin_q(n): if n in [2, 3]: return True a0, a1, a2 = 3, 0, 2 i = n while i > 0: a0, a1, a2 = a1, a2, (a0+a1) % n i = 1 return a0 == 0 # ====================================== n = 2 t = 1 while clock() < int(argv[1]): is_prime = is_prime_q(n) is_perrin = is_perrin_q(n) if is_prime != is_perrin: print n, is_prime, is_perrin if t < clock(): print 'Tested to %d in %d secs' % (n,t) t += 1 n += 1 
$$n\;\approx\;3150\sqrt(t)$$
So here we have the code to start the testing process. We have a simple test for primality, first testing for being a small prime, then using Fermat's Little Theorem, and finally performing trial division up to the square root of the number being tested.
Then we have the routine to test the Perrin condition. We start with the first few terms in the Perrin sequence, then advance up the sequence, working modulo n to avoid having the numbers get stupidly large. In the end, we return True if the relevant term in the sequence is zero, indicating that n divides k(n).
The main loop starts at n=2 and then tests to see if the number is prime, and whether it passes the Perrin test. If these answers don't match, then the number is printed, and is a counterexample.
Specific measurements give:
$$n\;\approx\;4575\sqrt{t}$$
This is all very rough and ready, and will fall apart as we go significantly further, but it's enough to get an idea of what's going on. We can see that getting to a million will take about 13 hours, so we really, really want to speed this up.
So what do we do? How can we get further, faster?
In Part 2 we'll look at algorithmic improvements.
The Unwise Update  :  FindingPerrinPseudoPrimes Part2 ... 
You should follow me on twitter 
I've decided no longer to include comments directly via the Disqus (or any other) system. Instead, I'd be more than delighted to get emails from people who wish to make comments or engage in discussion. Comments will then be integrated into the page as and when they are appropriate.
If the number of emails/comments gets too large to handle then I might return to a semiautomated system. That's looking increasingly unlikely.
Quotation from Tim BernersLee 