The Weekly Challenge ‐ Perl and Raku

CY's Take on The Weekly Challenge #155

If you want to challenge yourself on programming, especially on Perl and/or Raku, go to https://theweeklychallenge.org, code the latest challenges, submit codes on-time (by GitHub or email).

Do tell me, if I am wrong or you strongly oppose my statements!

It's time for challenges in Week #155 !

Two number theory tasks.

Task 1: Fortunate Numbers

(From Wikipedia) One of unsolved problems in mathematics:
Fortune's conjecture: Are any Fortunate numbers composite?

I have never thought that things like primorial can go with such an interesting and natural conjecture.

Back to our task. Looking at the wikipedia page linked from the task statement, we know that there are TWO sequences of Fortunate numbers:

The first one is a1(n) is the smallest integer m so that m + primorial(n) is prime.

The second one is sort {$a<=>$b} uniqint a1 (or equivalently, uniqint sort {$a<=>$b} a1 ^[1]).

The task requires us to generate the second sequence.

A concern is, a1(x) can be smaller than a1(y) when x > y. So how do we know when to stop generating terms of a1?

Keep reading on the Wikipedia gives us a hint: "The Fortunate number for primorial(n) is always above p_n."

So after we get a length-$N list of non-duplicate values of a1 we grasp its largest element $E, and check whether there are Fortunate numbers smaller than $E. The check will be stopped when we reach j, where j is determined from the inequality p_j ≤ $E < p_j+1 .

The worry has been cleared and we should start to code:

The first part, to generate a1 up to a certain "large-enough" value:

# first part

{
    my $k = scalar @primorials - 1;
    my $tmp_int = $primes[$k];
    do {
        $tmp_int++;
    } while (!is_prime( $primorials[-1] + $tmp_int ));
    push @fort_unsort, $tmp_int;
    iter_primorial();    # generate the next primorial
    @fort_sorted = uniqint sort {$a<=>$b} @fort_unsort;
}

Terms in @primorials grow fast and I set them to be Math::BigInt objects. This is the reason that I generate it one-by-one ‐ avoid eating up memory.

The second part, to check whether there are any small Fortunate numbers not yet included in @fort_unsort:

# second part

{
    # ... 
    my $k = scalar @primorials - 1;
    my $tmp_int = $primes[$k];
    do {
        $tmp_int++;
    } while (!is_prime( $primorials[-1] + $tmp_int ));
    push @fort_unsort, $tmp_int;
    @fort_sorted = uniqint sort {$a<=>$b} @fort_unsort;
}

I try to keep the size of @primes and @primorials same. And at certain moment I felt the name of @primes confusing, because it is not a full list of primes. Then I changed it into @ch_primes.

Here I reveal the loop conditional for the second part:

while ($fort_sorted[$N-1] > $ch_primes[-1])

The loop conditional for the first part, actually, can be quite arbitary. Each time the block runs, @fort_sorted is either increased by 1 or remain the same. One may casually use while (scalar @fort_sorted != 2*$N ) or while (scalar @fort_sorted < 2*$N ) , and s/he could even skip the second part ... IF LUCKY ENOUGH?

To go on a steady pace, I use

while (scalar @fort_sorted < $N) 

Since the first part and the second part look so similar, I combine them by a short-circuit or; here is the code:

my $N = $ARGV[0] || 8;

my @ch_primes = (2);
my @primorials = ( Math::BigInt->new(2) );
my @fort_unsort = ();
my @fort_sorted = ();


while ( !defined($fort_sorted[$N-1])
               ||
        $fort_sorted[$N-1] > $ch_primes[-1]
) {
    my $k = scalar @fort_unsort;
    my $tmp_int = $ch_primes[$k];
    do {
        $tmp_int++;
    } while (!is_prime( $primorials[$k] + $tmp_int ));
    push @fort_unsort, $tmp_int;
    iter_primorial() while scalar @primorials <= scalar @fort_unsort;
    iter_prime() while scalar @ch_primes <= scalar @primorials;
    @fort_sorted = uniqint sort {$a<=>$b} @fort_unsort;
}

say "Answer: \n", join ", ", @fort_sorted[0..$N-1];

while the supporting subroutines &iter_prime and &iter_primorial are

sub iter_prime {
    my $new_int = $ch_primes[-1];
    do {
        $new_int++;
    } while (!is_prime($new_int));
    push @ch_primes, $new_int;
}

sub iter_primorial {
    iter_prime() while scalar @ch_primes <= scalar @primorials; 
    push @primorials,
      $primorials[-1]*$ch_primes[ scalar @primorials ];
}

Sample run:

$ perl ch-1.pl 10
Answer: 
3, 5, 7, 13, 17, 19, 23, 37, 47, 59

---

Is the heuristic while (scalar @fort_sorted < 2*$N ) good? Let add a line say scalar @fort_unsort; and see:

---

I have drawn an alternative version of the code after reading polettix's blogpost on the task. We need not store the values of many primorials; we can just store one and proceed.

my @ch_primes = (2);
my $my_primorial = Math::BigInt->new(2);
my $ind_primorial = 1;
my @fort_unsort = ();
my @fort_sorted = ();


while ( !defined($fort_sorted[$N-1])
               ||
        $fort_sorted[$N-1] > $ch_primes[-1]
) {
    my $k = scalar @fort_unsort;
    my $tmp_int = $ch_primes[$k];
    do {
        $tmp_int++;
    } while (!is_prime( $my_primorial + $tmp_int ));
    push @fort_unsort, $tmp_int;
    iter_primorial() while $ind_primorial <= scalar @fort_unsort;
    iter_prime() while scalar @ch_primes <= $ind_primorial;
    @fort_sorted = uniqint sort {$a<=>$b} @fort_unsort;
}

say "Answer: \n", join ", ", @fort_sorted[0..$N-1];


sub iter_prime {
    my $new_int = $ch_primes[-1];
    do {
        $new_int++;
    } while (!is_prime($new_int));
    push @ch_primes, $new_int;
}

sub iter_primorial {
    iter_prime() while scalar @ch_primes <= $ind_primorial; 
    $my_primorial *= $ch_primes[ $ind_primorial ];
    $ind_primorial++;
}

An anecdote: we never assume the Fortunate conjecture true all along the code process, though, it's quite cool.

Task 2: Pisano Period

As usual I try to generalize. This time I generalize it into arbitary length of "order" of linear recurrence.

https://oeis.org/wiki/Index_to_OEIS:_Section_Rec

ch-2.pl $N 2   1 1       0 1   # Fibonacci numbers

ch-2.pl $N 2   1 1       2 1   # Lucas numbers    #OEIS:A000032

ch-2.pl $N 2   1 2       0 1   # Pell numbers     #OEIS:A000129

ch-2.pl $N 3   1 1 0  1 1 1 # Padovan numbers     #OEIS:A000931

ch-2.pl $N 3   1 1 1  0 0 1 # Tribonacci numbers  #OEIS:A000073

ch-2.pl $N 3   1 0 1  1 1 1 # Narayana's cows sequence #A000930

ch-2.pl $N 6  -1 -1 0 -1 2 1   1 2 2 4 5 9             #A001224

ch-2.pl $N 6   1 -1 -1 0 1 1   1 1 2 3 4 5             #A001399

So, how does the script work?

  1 sub pisano_period {
  2     my ($N, $t, $rec, $seq) = @_;
  3     die "(Some of) Parameters are too large.\n"
  4         if $N**$t + $t - 1 > 8_000_000;
  5 
  6     @$seq = map {$_ % $N} @$seq;
  7     my $ori_seqstate = [@$seq];
  8     my $new_seqstate = [@$ori_seqstate];
  9     my $count = 0;
 10     do {
 11         my $new_val = sum map {$rec->[$_]*$new_seqstate->[$_]} (0..$t-1);
 12         $new_val = $new_val % $N;
 13         push @{$seq}, $new_val;
 14         shift @{$new_seqstate};
 15         push @{$new_seqstate}, $new_val;
 16         $count++;
 17         die "Patterns not found\n" if $count > $N**$t + $t;
 18     } while (!cmp_num_arr($new_seqstate, $ori_seqstate));
 19 
 20     return [@$seq[0..$count-1]];
 21 }

line 02: Let's explain it using examples.

For the input ch-2.pl 7 3   1 0 1   1 1 1, the desired output is the info of the periodic behavior of the Narayana's cows sequence modulo 7.

$rec will be [1,0,1] (the sequence is defined by a(n-1)+a(n-3) for n ≥ 3) and $ori_seqstate will be [1,1,1] (a(0)=a(1)=a(2)=1).

For the input ch-2.pl 17 2   1 2   0 1, the desired output is the info of the periodic behavior of the Pell numbers modulo 17.

$rec will be [1,2] (the sequence is defined by 2 × a(n-1)+a(n-2) for n ≥ 2) and $ori_seqstate will be [0,1] (a(0)=0, a(1)=1).

line 03-04, line 17: $N**$t + $t - 1 is a lazy upper bound. By the pigeonhole principle in discrete math, there are $N**$t possible length-$t-sequence (while the possible values of members of the length-$t-sequence are 0 to $N-1).

line 11: I am curious if there is a more compact way to write it. Note that $rec->@* and $new_seqstate->@* have same number of elements.

line 14-15: They can be replaced by $new_seqstate = [ @$new_seqstate[1..$t-1] , $new_val ]; .

line 18: &cmp_num_arr is a rough supportive subroutine checking whether two lists are the same.

line 17: I have to admit that I have given up. Depending on the initial condition, $ori_seqstate may appear once; for instance, for a recurrence sequence specified by a(0)=a(1)=1 and a(n) = 4 a(n-1) + 4 a(n-2), then the code here will return a "Patterns not found" message for many $N. I have not done the details math of the necessary and sufficient condition for a linear recurrence sequence which surely generates a Pisano-like period.

Sample run:

$ perl ch-2.pl 16 2  1 1  0 1 # Fibonacci numbers 
1..4
length: 24
terms: 0, 1, 1, 2, 3, 5, 8, 13, 5, 2, 7, 9, 0, 9, 9, 2, 11, 13, 8, 5, 13, 2, 15, 1
ok 1 - Required Case (Period of Fibonacci numbers mod 3)
ok 2 - Pisano period for mod 2
ok 3 - Pisano-like period for Padovan numbers mod 3
ok 4 - Pisano-like period for Padovan numbers mod 7

Anyway, have fun exploring the periodic behavior of linear recurrence sequences.

Well, there are many fun discussions on the properties of Fibonacci sequences on Mathematics Stack Exchange, some related to the Pisano Period.

Unfortunately, I have been infected by COVID-19 Omicron this week ‐ tested with rapid tests. Three shots of vaccine were injected before. And three exams in April. I better speak less and take more rest. Everyone take care.

Stay alert and healthy everyone! Hope for peace and justice on East Europe! □

Remarks:

1: From my testing up to the highest member of the sequence is 1103, the performance between the different order of occurring of the two array operations seems to be very close.

Except from images and codes from other personnels, the content of this blogpost is released under a copyleft spirit. One may share (full or partial) content of this blogpost on other platform if you share it under the free and open content spirit.

link for CY's full codes: ch-1.pl, ch-2.pl

Contact on twitter: @e7_87.

Discuss via GitHub issues: here.

Email: fungcheokyin at gmail.com

Created Date: 10th March, 2022; typos and grammars fixed on 14th March.