Suggest a change

These fields are all optional and need only
be supplied if you would like a direct reply.

Subject
Your email address
Your real name

You must answer this!

If you don't, my spam filtering will
ensure that I never see your email.

What's 8 plus five (in digits only)?

Please make your changes here and then
Editing tips and layout rules.

File: FindingPerrinPseudoPrimes_Part2

''' <link rel="alternate" type="application/rss+xml"
'''       href="/rss.xml" title="RSS Feed">

********> width="25%"
|>>
''' <a title="Subscribe to my feed"
'''    rel="alternate"
'''    href="https://www.solipsys.co.uk/rss.xml">
'''    <img style="border-width: 0px;"
'''         src="https://www.feedburner.com/fb/images/pub/feed-icon32x32.png"
'''         align="middle"
'''         alt="" />Subscribe!</a>
_
''' <a href="https://twitter.com/ColinTheMathmo">
''' <img src="https://www.solipsys.co.uk/new/images/TwitterButton.png"
'''      title="By: TwitterButtons.net"
'''      width="212" height="69"
'''      alt="@ColinTheMathmo"
'''     /></a>
<<|
----
My latest posts can be
found here:
  * ColinsBlog
----
Previous blog posts:
  * FindingPerrinPseudoPrimes_Part1
  * TheUnwiseUpdate
  * MilesPerGallon
  * TrackingAnItemOnHackerNews
  * HackerNewsUserAges
  * PokingTheDustyCorners
  * ThereIsNoTimeForThis
  * PublicallySharingLinks
  * LearningTimesTables
  * GracefulDegradation
  * DiagrammingMathsTopics
  * OnTheRack
  * SquareRootByLongDivision
  * BeyondTheBoundary
  * FillInTheGaps
  * SoftwareChecklist
  * NASASpaceCrews
  * TheBirthdayParadox
  * TheTrapeziumConundrum
  * RevisitingTheAnt
  * TheAntAndTheRubberBand
  * IrrationalsExist
  * MultipleChoiceProbabilityPuzzle
  * RandomEratosthenes
  * WrappingUpSquareDissection
  * DissectingASquarePart2
  * DissectingACircle
  * DissectingASquare
  * AnOddityInTennis
  * DecisionTreeForTennis
  * DecisionTreesInGames
  * AMatterOfConvention
  * DoYouNourishOrTarnish
  * BinarySearchReconsidered
  * TwoEqualsFour
  * TheLostPropertyOffice
  * TheForgivingUserInterface
  * SettingUpRSS
  * WithdrawingFromHackerNews
----
Additionally, some earlier writings:
  * RandomWritings.
  * ColinsBlog2010
  * ColinsBlog2009
  * ColinsBlog2008
  * ColinsBlog2007
  * ColinsBlogBefore2007
''' <img src="/cgi-bin/CountHits.py?FPPP_2" alt="" />
********
!!! Finding Perrin Pseudo-primes: _ Part 2 - 2015/05/17

[[[>
{{{
#!/usr/bin/python

n = 2
t = 10

t_prime = 0.0
t_perrin = 0.0

while clock() < int(argv[1]):

t_prime -= clock()
    is_prime = is_prime_q(n)
    t_prime += clock()
    t_perrin -= clock()
    is_perrin = is_perrin_q(n)
    t_perrin += clock()
    if is_prime != is_perrin:
        print n, is_prime, is_perrin

if t < clock():
        print 'Tested to %d in %d secs' % (n,t)
        t += 10

n += 1

print 'Time in prime :', t_prime
print 'Time in Perrin:', t_perrin

# ================
# Edited output:
#
# Tested to 42763 in 100 secs
#
# Time in prime :  0.11
# Time in Perrin: 99.72
#                 =====
#                 99.83
}}}
]]]
So now we've got the scaffolding of a program to find these
Perrin Pseudo-Primes. Here is the main loop of the code, with
some simplistic timing added to it.  (Note that this code is
incomplete and won't run as is).

The output shows that when given 100 seconds to run it gets as
far as /n=42763,/ but more importantly, the timing shows that
overwhelmingly it spends its time in the routine to test whether
or not a number passes the "Perrin Test."  So there are a few
things we need to do:

* Make fewer calls to the Perrin Test
  * Make the Perrin Test run faster.

In this post we'll look at a simple technique to accomplish
the first of these.

There's a maxim in programming:

|>> [[[ |>>
You can't make a program _
run faster, you can only _
make it do less.
<<| ]]] <<|

One way of making the program do less is to pre-compute things
and then use the result multiple times.  Another way is to do
a quick test to avoid a longer one.  We'll do both of those.

Here's the critical observation:

* If /m/ divides /n/,
    * and /n/ divides /k(n),/
      * then /m/ divides /k(n)./

The second part of the observation is that the Perrin sequence
modulo /m/ has to repeat itself after at most EQN:m^3 steps.
There are only EQN:m^3 possible triples, and the sequence is
reversable, so eventually we get the first triple appearing
again.  Here are the sequences for 2, 3, and 5, up to the point
where they start to repeat

* mod 2 : 1 0 0 1 0 1 1 1 0 0 ...
    * repeats after 7 terms
  * mod 3 : 0 0 2 0 2 2 2 1 1 0 2 1 2 0 0 2 ...
    * repeats after 13 terms
  * mod 5 : 3 0 2 3 2 0 0 2 0 2 2 2 4 4 1 3 0 4 3 4 2 2 1 4 3 0 2 ...
    * repeats after 24 terms
[[[>
{{{
#!/usr/bin/python

def compute_perrin_number(n):

if n in [2, 3]: return 0
    a0, a1, a2 = 3, 0, 2
    i = n
    while i > 0:
        a0, a1, a2 = a1, a2, (a0+a1) % n
        i -= 1

return a0

def make_pre_filter(m):

rtn = [ 3 % m, 0, 2 % m ]
    while True:
        rtn.append( (rtn[-3] + rtn[-2]) % m )
        if rtn[:3] == rtn[-3:]:
            return rtn[:-3]

pre_filters = {}
pre_filter_elements = range(2,104)
for m in pre_filter_elements:
    pre_filters[m] = make_pre_filter(m)

def is_perrin_q(n):

for m in pre_filter_elements:
        if n % m == 0:
            pre_filter_to_use = pre_filters[m]
            ndx = n % len(pre_filter_to_use)
            if 0 != pre_filter_to_use[ndx]:
                return False

return compute_perrin_number(n) == 0
}}}
]]]
So what we can do is this.  For each of several primes (and maybe
even some composites) we precompute the Perrin sequence modulo
that number.  Then when we are looking at a number /n/ we look at
all the precomputed sequences that divide /n/ and ask if the
sequence has a 0 in the right place.  If not, then /n/ will not
divide /k(n)./  Of course, if it does have a 0 in the right place
for every precomputed /m/ that divides /n/ then we will still have
to do the computation, but with any luck that will be a lot less.

So here is the result.

* Without the sieving:
    * Tested to 45761 in 100 secs
    * Perrin computations: 45760

* With sieving from m=2 to 103:
    * Tested to 109159 in 100 secs
    * Perrin computations: 19004

As you can see, we get more than twice as far, so we're clearly
doing a good thing.  It's obvious why when we look at the number
of Perrin computations we do.  Without the sieving every number
has the full calculation, but with sieving only 19 in every 109
gets the calculation, or roughly 17.4%, a saving of 82.6%.

But there's still a problem that we can see clearly when we look
at the timing:

* Time in prime : 0.25 secs
  * Time in Perrin : 99.47 secs

Despite saving over 80% we're still spending all our time in the
Perrin test routine. In the next installment we'll see how we can
speed that up further by using matrices and an adaptation of the
Peasant Multiplication Algorithm.
----

|>>
| |>> <<<< Prev <<<< ---- FindingPerrinPseudoPrimes_Part1 <<| | : | |>> >>>> Next >>>> ---- RussianPeasantMultiplication <<| |
----
********>
''' <a href="https://twitter.com/ColinTheMathmo">You should follow me on twitter</a>
********
''' <a href="https://twitter.com/ColinTheMathmo">
''' <img src="https://www.solipsys.co.uk/new/images/TwitterButton.png"
'''      title="By: TwitterButtons.net"
'''      width="212" height="69"
'''      alt="@ColinTheMathmo"
'''     /></a>
********<
<<|
----
!! 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.

********<