The Weekly Challenge ‐ Perl and Raku

CY's Take on The Weekly Challenge #169 ‐ Number Theoretic Quiz

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 #169 !

Both tasks of this week ask for the first 20 terms of a certain integer sequence; I ran an extra mile to make it suitable for some larger values. However my running performance has been so-so.

Task 1: Brilliant Numbers

This is a simple revision of my code.

#!/usr/bin/perl
# The Weekly Challenge 169
# Task 1 Brilliant Numbers # CY's solution version 2.0
use v5.24.0;
use warnings;

my $req = $ARGV[0] || 20;

my @prime = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,
             41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97);

die "You're asking too many, or I am too little!\n"
    if $req > 241;  # 21*20/2 +  21 + 4*3/2 + 4


my $pt = ceil( ( sqrt(8*$req-79) +7 )/2 ) - 1;

my @brilliant = (-1, 4, 6, 9, 10, 14, 15, 21, 25, 35, 49);

my @temp;


for my $i (4..$pt) {
    for my $j ($i..24) {
        my $product = $prime[$i]*$prime[$j];
        push @temp, $product;
    }
}

push @brilliant, sort {$a<=>$b} @temp;

say join ", ", @brilliant[1..$req];

Sample Run:

$ ch-1.pl 30
4, 6, 9, 10, 14, 15, 21, 25, 35, 49, 121, 143, 169, 187, 209, 221, 
247, 253, 289, 299, 319, 323, 341, 361, 377, 391, 403, 407, 437, 451

(The 20th term is bold.)

Improvisation

FULL CODE

Declaring myself as an ambitious person, here I extend the methods up to "brilliant products" by 4-digit primes.

Number of brilliant numbers from x-digit primes is equal to the number of prime squares and the number of the pairs, for instance, 4² + 4*3/2 = 10, 21² + 21*20/2 = 231.

There are many pre-processing for this version.

sub bn {
    my $req = $_[0];
    my $show_all_values = $req <= 50;
    # show all values if the required number is small

    my %prime = (1 => [2, 3, 5, 7],
                 2 => [11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
                       47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97],
                 3 => [],
                 4 => [],
                );

    my $digit_len;

    given ($req) {
        $digit_len = 1 when $_ <= 10;
        $digit_len = 2 when $_ <= 241;
        $digit_len = 3 when $_ <= 10537;
        $digit_len = 4 when $_ <= 573928;
        default {die "You're asking too many!\n";}
    }

    eval { forprimes {push $prime{3}->@*, $_}  100, 999; } if $digit_len >= 3;
    eval { forprimes {push $prime{4}->@*, $_} 1000,9999; } if $digit_len >= 4;

    my $ind = $req;
    $ind = $req-10 if $digit_len == 2;
    $ind = $req-241 if $digit_len == 3;
    $ind = $req-10537 if $digit_len == 4;
    my $c = 3;
    $c = ceil( ( sqrt(8*$ind+1) - 1 ) / 2 ) - 1 if $digit_len >= 2;

    my @temp;

    for my $i (0..$c) {
        my $I = $prime{$digit_len}[$i];
        my $_j = scalar $prime{$digit_len}->@* - 1;
        for my $J ($prime{$digit_len}->@[$i..$_j]) {
            my $product = $I*$J;
            push @temp, $product;
            last if $product > a_function_of_$c_or_$prime{$digit_len}[$c];
            # It would be great if I could find such a function!
        }
    }
    @temp = sort {$a<=>$b} @temp;
    if ($show_all_values) {
        my @brilliant = (-1, 4, 6, 9, 10, 14, 15, 21, 25, 35, 49);
        push @brilliant, @temp;
        say join "\n", @brilliant[1..$req];
        return;
    }
    return $temp[$ind-1];
}

( sqrt(8*$ind+1) - 1 ) / 2 comes from solving quadratic inequality. To act as a writer of math textbook, I leave the derivation as an exercise... Huh!?

OK. I am not good at acting. I explain. See my $pt = ceil( ( sqrt(8*$req-79) +7 )/2 ) - 1; from the simple version. It comes from 10 + ($pt-4)*($pt-5)/2+$pt-4 >= $req.

( sqrt(8*$ind+1) - 1 ) / 2 comes from the same principle, just we do not care how many brilliant numbers are from products of primes with less digits, because we set $ind equal to $req minus number_of_brilliant_numbers_from_products_of_primes_with_less_digits.

Grade 10 math. :)

Finally it is the testing part:

use Test::More tests=> 14;
ok bn(60) == 841;
ok bn(80) == 1079;
ok bn(100) == 1369;
ok bn(120) == 1763;
ok bn(140) == 2231;
ok bn(160) == 2809;
ok bn(241) == 9409;
ok bn(5000) == 197503;
ok bn(10000) == 696191;
ok bn(10416) == 851927;
ok bn(10537) == 994009;  # (OEIS limit)
ok bn(10538) == 1009*1009;
ok bn(573927) == 9967*9973;
ok bn(573928) == 9973*9973;

# ref: OEIS:A078972

All are done! Sample Run:

$ ch-1.pl 2022
1..14
85207
ok 1
ok 2
ok 3
ok 4
ok 5
ok 6
ok 7
ok 8
ok 9
ok 10
ok 11
ok 12
ok 13
ok 14

Lessons Learnt in Program Design

The Version 1 is here. The difference between v1 and v2 is essentially only the for my $i (0..XX). In v1, it is a constant 24; in v2, the value of $pt is determined by the quadratic formula. When I started to code this task, the tricky point for me is how NOT to generate all brilliant products of x-digit primes. Then I divided it into two parts to be tackled: (1) find a prime P so that "(almostly-)enough" brilliant products are generated, and then (2) compute the brilliant products of "small enough" x-digit primes and large x-digit primes in order to include some small brilliant numbers from their products; combine the result, sort it, and locate the $req-th brilliant number.

The red part of the code posted made me headache. Though I haven't got a criterion for last, the idea of (2) is still valuable. This is the second version of this blogpost as after I read my own blogpost, I cleaned up my mind and realized I should not abandon to reduce the number of computations just because I could not figure out a good last statement. "Don't throw the baby out with the bathwater."

Lessons Learnt in Perl Syntax

The rewriting has forced me to get more familiar with many stuff:

• Perl's given-when, which is a bit different from switch in many other programming languages;
• @prime{$digit_len}[$i..$_j] is invalid... $prime{$digit_len}->@[$i..$_j] is valid;
• forprimes in Math::Prime::Util; usually only using subroutines in modules, I seldom use this kind of contruct (loop construct) from module.

The current code is a bit messy due to a desire to fulfill task requirements (show all first 20 values) and my ambition of getting n-th brilliant number for large n.

Anyway. Good task to be practised!

Task 2: Achilles Numbers

I did not investigate deeply and have just followed the mathemational definition to check out the required number of Achilles numbers.

-1+72*($prime[-1]+2 on the first line of the outermost for-loop is a bold stroke; originally I set the upper bound to 100_000 and set the terminate condition as $req > 21.

A short note:

Though using !($dividend % $divisor) looks cooler than using ($dividend % $divisor == 0), I have chosen the latter for readability and maintainability.

Code logic:

Preload a reasonable size of list of prime numbers from 2. The script checks the positive integers one by one. Firstly it checks whether the integer is "powerful" (from the task statement: "...for every prime factor p of n, p² is also a divisor."), then it collects the indices of the prime factors in the prime factorization formula of the integer. Finally it checks whether the integer is "imperfect", i.e. the indices do NOT have a common factor larger than 1.

"Short-circuits" to a next number, when a hint that the current integer is not an Achilles number, disperse throughout the code lines.

my $req = $ARGV[0] || 20;

my @prime = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,
             41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97);

my @achilles = (-1);


my $found = 0;

for my $i (2..-1+72*($prime[-1]+2)) {
    my $powerful = 1;
    my $imperfect;
    my @prime_factors;
    for (@prime) {
        if ($i % $_ == 0) {
            push @prime_factors, $_;
            $powerful = 0 if ($i % ($_*$_) != 0);
        }
        last if !$powerful;
    }
    next if !$powerful;
    next if scalar @prime_factors == 0;    # a large prime
    next if scalar @prime_factors == 1;
    # the above: made up of only one prime in @prime and a large prime
    my @ppower;
    for my $pf (@prime_factors) {
        my $ind = 2;
        my $pp = $pf*$pf*$pf;
        while ($i % $pp == 0) {
            $ind++;
            $pp = $pp*$pf;
        }
        push @ppower, $ind;
    }
    $imperfect = gcd(@ppower);

    if ($imperfect == 1) {
        push @achilles, $i;
        $found++;
    }
    last if $found == $req;
}

say join ", ", @achilles[1..$found];

die "$found terms are shown. You're asking too many, or I am too little!\n"
    if $req > $found;

sub gcd {
    # skip details... Euclidean gcd algorithm
}

sub gcd_pair {
    # skip details... Euclidean gcd algorithm
}

Sample Run:

$ ch-2.pl 30
72, 108, 200, 288, 392, 432, 500, 648, 675, 800, 864, 968, 972, 1125,
1152, 1323, 1352, 1372, 1568, 1800, 1944, 2000, 2312, 2592, 2700, 2888,
3087, 3200, 3267, 3456

# checked against OEIS:A052486

Stay alert and healthy! □

The image of hexagonal cross section cut from a cube is a public domain picture made by CY 16 years ago.

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: 18th June, 2022.

Last Edited Date Time: 20th June, 2022 AM06:48:23 HKT