 Subscribe!
My latest posts can be found here: Previous blog posts:
Additionally, some earlier writings: |
Finding Perrin Pseudo-primes:
|
|
This isn't an exact quotation, but it's close enough. |
- "A chap at work who does maths, but is not actually
a mathematician, asked me the following.
Consider the sequence $k(n)$ with- $k(1)=0,\;k(2)=2,\;k(3)=3,$ and
- $k(n)=k(n-2)+k(n-3).$
|
counter-example. |
So I started on that. I figured Andrew had already done the basic work of looking for a small counter-example, 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 pre-internet 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 counter-examples, 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(n-1)+b.k(n-2)+c.k(n-3)$$ 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:
- $\alpha\approx{}1.324736\ldots$
- $\beta\approx{}-0.66236-0.56228{}i$
- $\gamma\approx{}-0.66236+0.56228{}i$
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:
- When n is prime,
- does n divide k(n) ?
- does n divide k(n) ?
|
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 counter-example 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,n-1,n) != 1: return False
if n % 3 == 0 : return False
if pow(3,n-1,n) != 1: return False
if n % 5 == 0 : return False
if pow(5,n-1,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 counter-example.
Specific measurements give:
- Roughly 104 per second when n=100k
- Roughly 42 per second when n=250k
- Roughly 10 per second when n=1000k
$$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 |
|
Comments
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 semi-automated system. That's looking increasingly unlikely.
Suggest a change ( <--
What does this mean?) /
Send me email
