PWC 168 › Perrin Primes

This post is part of a series on Mohammad Anwar’s excellent Weekly Challenge, where hackers submit solutions in Perl, Raku, or any other language, to two different challenges every week. (It’s a lot of fun, if you’re into that sort of thing.)

Task #1 this week asks us to find the first 13 Perrin primes. “What’s a Perrin prime,” I can’t hear you asking? To answer that, we first have to look at the Perrin sequence, as described in OEIS A001608. It’s easy to generate:

Starting with [3, 0, 2], each new term is determined by adding the 2^nd and 3^rd last terms. So, the 4^th number is 3 + 0 = 3, giving us [3, 0, 2, 3]. The 5^th number is 0 + 2 = 2, and so on.

Perrin primes are simply the elements of the Perrin sequence that also happen to be prime.

Normally (and per the example output in the task) we are to find the unique Perrin primes, in order. So, we’ll just seed the first prime in our @r results, and then rely on the fact that the sequence is strictly increasing after the first five terms.

Building the sequence is very simple from there:

use Math::Prime::Util qw< is_prime >;

# Calculate the first $n Perrin primes
sub perrin {
    my ($n2, $n1, $n0) = (3, 0, 2); # Sliding window
    my @r;

    while (@r < $_[0]) {
        push @r, $n0 if $n0 > $r[-1] and is_prime($n0);
        ($n2, $n1, $n0) = ($n1, $n0, $n2 + $n1);
    }

    @r;
}

You’ll note I do not store the entire Perrin sequence—just the three most recent terms. When calculating higher Perrin primes, the numbers get large very quickly, and memory usage becomes a real problem, and that slows down the execution markedly (assuming it doesn’t run out of memory!) Storing the Perrin primes (@r) is not a significant issue, but those could certainly be printed and discarded as we go, instead.

ch-1.pl source code

The Raku code is very similar to its Perl counterpart, here:

sub perrin(Int $max) {
    my ($n2, $n1, $n0) = (3, 0, 2);
    my @r = (2);

    while (@r.elems < $max) {
        @r.push($n0) if $n0 > @r.tail and $n0.is-prime;
        ($n2, $n1, $n0) = ($n1, $n0, $n2 + $n1);
    }

    @r;
}

ch-1.raku source code

Efficiency

I’m already performing a constant number of calculations at each step, so there isn’t anything to optimize that will lower the complexity class for our task of returning the first n Perrin primes. I’m leaning on Math::Prime::Util for the well-implemented is_prime() function, but our best methods are still polynomial.

Optimizing memory usage actually turns out to be more important, which is why I’m not storing previous Perrin terms.

Growth

The 13^th Perrin prime is already 792.6 trillion, and the 14^th is another million times larger than that, which already requires bigint in order to calculate. If we look only at the Perrin numbers themselves (not Perrin primes), we can get an idea of how quickly they grow.

Whenever you are adding previous terms, you will get exponential growth. The Fibonacci sequence is a classic example of this—each term is bigger than the previous term by a factor of the Golden ratio, or \(\varphi = \frac{ 1 + \sqrt{5} }{2} \approx 1.618\,033 \dots \)).

Nailing down the exact growth (α, giving αⁿ = size of n^th Perrin number) is similar. As a sanity check, we can see \(1 < \alpha < 2 \). Each new term is less than twice the previous, since we’re adding the 2^nd and 3^rd last terms). However, we can do better than those bounds. The math is a bit beyond the scope of this post, but as it turns out, the growth is equal to:

\(\alpha = \displaystyle{ \frac{ \sqrt[3]{108+12 \sqrt{69}} + \sqrt[3]{108 – 12 \sqrt{69}} }{6} } \\
\alpha \approx 1.324\,717\,957\,244\,746\,025\,960 \dots \)

Even though 1.32ⁿ isn’t a terribly large exponential, by the 50^th term, we are already well over a million.

It’s nice when a little bit of research and math can definitively reveal whether it’s a good time to interrupt my program or not!

The Sequence Itself

Here are as many terms of this sequence as my computer could generate while I wrote this blog post, and had lunch:

perrin(   3) is prime# 1 =         2 = 2
perrin(   4) is prime# 2 =         3 = 3
perrin(   6) is prime# 3 =         5 = 5
perrin(   8) is prime# 4 =         7 = 7
perrin(  11) is prime# 5 =        17 = 17
perrin(  13) is prime# 6 =        29 = 29
perrin(  21) is prime# 7 =       277 = 277
perrin(  22) is prime# 8 =       367 = 367
perrin(  25) is prime# 9 =       853 = 853
perrin(  35) is prime#10 =  1.42e+04 = 14197
perrin(  39) is prime#11 =  4.37e+04 = 43721
perrin(  76) is prime#12 =  1.44e+09 = 1442968193
perrin( 123) is prime#13 =  7.93e+14 = 792606555396977
perrin( 167) is prime#14 =  1.87e+20 = 187278659180417234321
perrin( 237) is prime#15 =  6.62e+28 = 66241160488780141071579864797
perrin( 356) is prime#16 =  2.26e+43 = 22584751787583336797527561822649328254745329
perrin( 357) is prime#17 =  2.99e+43 = 29918426252927024136988188355201180399482197
perrin( 931) is prime#18 = 3.76e+113 = 375650352810749628391658393147651164149079195002314045738061982119710039976648976965060598469931973177804611901813
perrin(1043) is prime#19 = 1.79e+127 = 17889871724792219792241402014701050416254403054909819082963323121939408639274412767017724313639101409409795922319558694157739957
perrin(1215) is prime#20 = 1.81e+148 = 18106564606349058350871445556416183706383627605153862231876341960946635847221883756661544295450957270512362655785866338801138896957806303459431839801
perrin(1462) is prime#21 = 2.64e+178 = 26443665126671039192963010650954408309392693422822076064578125303560832561672888722088906692033449248344378194605701099265071815485284432217750405098433434144179485693746031340517
perrin(1623) is prime#22 = 1.21e+198 = 1213927704065079865017068478668276043626477148780514011765015731886286159454983721480068033892046357328417429372450987777059793416910075913180181245051185193201551033755831307534780351082477949347441
perrin(4431) is prime#23 = 1.02e+541 = 10157009252817374963867100949951608928714862242745008993013668540854184107874570804968501397570379041274564782116665400515007185872727535557633044532545504298417513460010708859590624519737577132699128528530905048976280744785692707368299964909111445284217209819026508401682969213029773502708330316828337469827393737449858748826727773566201071908906567992775961863663250545199702810339801764180200620104056601639153965055826816646056412891949330608030933040756303987388596508709113305229398404925505186056022798817893091541647706591557044644581

The largest Perrin prime I could find reference to is Perrin(581133), which, even if I forego the is_prime() check, takes a long time to compute.