laurent_r

These are some answers to the Week 8 of the Perl Weekly Challenge organized by Mohammad S. Anwar.

Challenge # 1: Perfect Numbers

Write a script that computes the first five perfect numbers. A perfect number is an integer that is the sum of its positive proper divisors (all divisors except itself). Please check Wiki for more information. This challenge was proposed by Laurent Rosenfeld.

Since Mohammad Anwar revealed that I suggested this challenge, let me add a few things about it. The idea of this challenge comes from a very long memory: the requirement is almost the same as a programming assignment I had during the first year (actually the first few months) of my CS study back in 1991, except for the fact that we were asked to write a program (in Pascal at the time) to find the perfect numbers below 10,000 (i.e. 4 perfect numbers). If a student can do it in the very first months of his or her CS curriculum, then it would certainly be too easy for a challenge aimed at experienced CS professional. This the reason why my suggestion was to find the first 5 perfect numbers (instead of 4) as this makes a real difference: as we will see, finding the fifth perfect number is significantly more difficult that finding the first 4 ones.

Let me also be clear on something: 1 is not a perfect number, because the definition says the sum of all its divisors except itself. The sum of all divisors except itself of 1 is 0; thus, 1 is not a perfect number. The first perfect number is 6 ( = 1 + 2 + 3).

Perfect Numbers in Perl 5

Brute Force Approach

Let's try brute force approach to grasp the problem. I'll start with my assignment of 1991, i.e. finding the perfect numbers below 10,000:

#!/usr/bin/perl
use strict;
use warnings;
use feature "say";

sub divisors {
    my $num = shift;
    my $limit = 1 + int $num / 2;
    my @divisors = (1);
    for my $test_div (2..$limit) {
        push @divisors, $test_div if $num % $test_div == 0;
    }
    return @divisors;
}

for my $num (2..10000) {
    my $sum_div = 0;
    my @divisors = divisors $num;
    for my $div (@divisors) {
        $sum_div += $div;
    }
    say " $num : @divisors" if $sum_div == $num ;
}

Easy, nothing complicated. I made a separate divisors subroutine because we may later want to reuse it. If we run that program, we get the following output:

$ time perl perfect1.pl
6 : 1 2 3
28 : 1 2 4 7 14
496 : 1 2 4 8 16 31 62 124 248
8128 : 1 2 4 8 16 32 64 127 254 508 1016 2032 4064

real    0m2.421s
user    0m2.343s
sys     0m0.062s

Two and a half seconds to find the first four perfect numbers. It probably took at least ten or twenty times that duration when I completed this assignment with the hardware of 1991 (my computer at the time didn't even have a hard disk, just two 5.5-inch floppy drives), but you would still get the result in a relatively reasonable time.

So, perhaps we can try with numbers up to 100,000, maybe we'll find the next perfect number. Well, that does not work. Not only we do not find any new perfect number with this new limit, but with 10 times more numbers to visit, the program takes almost 4 minutes to execute:

$ time perl perfect1.pl
6 : 1 2 3
28 : 1 2 4 7 14
496 : 1 2 4 8 16 31 62 124 248
8128 : 1 2 4 8 16 32 64 127 254 508 1016 2032 4064

real    3m50.514s
user    3m48.265s
sys     0m0.140s

This program really does not scale well: with ten times more numbers to study, it takes more or less 100 more time to run. The problem stems from the fact that we have two nested loops (even if it's a bit hidden by the fact that one of the loops is in the subroutine and the other calling that subroutine). Note that I very well knew that something like this was going to happen, but nonetheless made the test to get actual timings to illustrate the problem. We are not going to continue in this direction and try with a limit of one million, as this would obviously take hours to run (with no guarantee to find the next perfect number).

Triangular Numbers

When we were given this assignment back in 1991, not only was the hardware much slower than today, but the Web did not exist (well, it had just been created but was still an internal tool of the CERN, it went public only in 1993) and there certainly wasn't something like Wikipedia (which was created in 2001) to find additional information about perfect numbers. Our professor had been quite wise to ask only for the first four perfect numbers.

Looking at the Wikipedia page on perfect numbers, we can find a lot on interesting information on perfect numbers. First, all known perfect numbers are even (but it is not known whether there are odd perfect numbers, it can only be said that none has been found so far). Assuming that we are looking only for even numbers will cut the duration by half, but that's far from sufficient to solve the performance problem.

The article provides an additional piece of information: all even perfect numbers are triangular numbers. This is a list of triangular numbers:

0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91, 105, ..., 1225, 1275, 1326, 1378,

You start with 0. Add 1 to it to get the next one (1). Add 2 to second one to get the third one (3). Add 3 to the third one (6). And so on. This is interesting because there becoming increasingly sparce and we'll need to check the divisors of much fewer numbers.

So, the forloop at the end of the previous program is now changed as follows:

my $triangular_num = 1;
for my $num (2..10000) {
    $triangular_num += $num;
    my $sum_div = 0;
    my @divisors = divisors $triangular_num;
    for my $div (@divisors) {
        $sum_div += $div;
    }
    say " $triangular_num : @divisors" if $sum_div == $triangular_num;
}

This is indeed much faster than before, and I was able to get the fifth perfect number after about 45 minutes:

$ time perl perfect.pl
6 : 1 2 3
28 : 1 2 4 7 14
496 : 1 2 4 8 16 31 62 124 248
8128 : 1 2 4 8 16 32 64 127 254 508 1016 2032 4064
33550336 : 1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8191 16382 32764 65528 131056 262112 524224 1048448 2096896 4193792 8387584 16775168

and it took more than two hours to complete to the last value of the for loop (without finding any additional perfect number).

So, mission accomplished? Well, yes, in a way, we've got our first five perfect numbers, but I'm not really satisfied with this solution: this is way too long.

Mersenne Numbers

The Mersenne number are powers of 2 minus one. More precisely, for every integer number n, the corresponding Mersenne number is the number (2 ** n) -1. Their name is associated with Marin Mersenne, a French monk living in the seventeenth century. These numbers had been studied long before, notably in antiquity by Greek mathematician Euclid of Alexandria, but Mersenne studied them in detail in the context of primality of large numbers, and he may have demonstrated that a Mersenne number (2 ** n) -1 can be prime only if n is also prime. The opposite is not true: for example, 11 is prime, but (2 ** 11) -1 is not prime. Mersenne prime numbers are still very important in the search of very large prime numbers: as of 2018, the seven largest known prime numbers are Mersenne primes.

Euclid discovered and demonstrated an interesting property about the Mersenne primes (although they were obviously not called like that at the time) and perfect numbers: if a Mersenne number (2 ** p) -1 is prime (which requires p to be also prime), then the number ((2 ** p) -1) * 2 ** (p - 1) is a perfect number. In the eighteenth century, Swiss mathematician Leonhard Euler proved conversely that every even perfect number has this form. This is sometimes known as the Euclid Euler theorem.

This can give us a much faster way to identify perfect numbers: compute Mersenne numbers (2 ** p) -1 with p prime, check if such Mersenne number is itself prime, and if it is, that we know that ((2 ** p) -1) * (2 ** (p - 1)) is a perfect number.

For this, we'll write two suboutines, find_primes, which will populate an array of prime numbers, and is_prime, which will determine if its argument is a prime numbers. Contrary to what some readers might expect, is_prime is not going to be used to populate the list of prime numbers, but things will happen the other way around: is_prime will use the list of primes found by the first subroutine. This may seem counter-intuitive, but how do you find out if a number is prime? The basic test is to try to divide it by all numbers smaller than it and if it is not evenly divided by any of these smaller numbers, then it is prime. For a large number, this can take quite a bit of time. This basic idea can be improved very significantly in two ways. First, we don't need to test even division by all integers smaller than the target number being tested, but we can limit the test to numbers smaller than the square root of the target number. So, for example, if you want to know whether 10,001 is prime, you don't need to test divisibility by all numbers smaller than 10,001, but only by numbers smaller than 100, a gain of a factor of 100. Second, you don't need to test all numbers smaller than the square root of the target number. For example, if you start by trying to divide it by 2, you'll find that 10,001 is not evenly divided by 2. This also means that the target number is not evenly divided by any multiple of 2, so when you continue you can skip the test for all even numbers, 2, 4, 6, 8, 10, etc. Similarly, the next test would be to try to divide by 3. Once you've found that 10,001 is not divisible by 3, you can also skip test with multiples of 3, i.e. 6, 9, 12, 15, etc. This method is called the sieve of Eratosthenes. Pushing the argument further, what you actually need is to test divisibility by prime numbers smaller than the square root of the target number.

So, for our is_prime subroutine to be efficient with a large number, we need to have first a list of all prime numbers smaller than the square root of that larger number we're likely to investigate. Creating this list will take some time, but it will make the is_prime subroutine much faster, and this will be a net gain if we're going to test primality of a relatively large quantity of large numbers. The is_prime subroutine is doing just that: testing divisibility by all primes smaller than the square root of the target number. Note that there are some more efficient primality tests, such as the Miller-Rabin test (or, in the case of Mersenne numbers, the Lucas–Lehmer primality test), but implementing such tests is beyond the scope of this blog post.

The find_primes subroutine uses the same principle, with an added small tweak: it populates the @primes array and use it at the same time, i.e. it uses the list of primes already generated to find out whether a new number is a prime. This can be a bit tricky for very small numbers (especially in view of the fact that 2 is a very special case, as it is the only number that is both even and prime), so we pre-populate the @primes array with the first three primes (2, 3 and 5) and start from there. This subroutine has been previously used and described in a previous Perl Weekly challenge, so you can check there if you want more information on it.

Once we have these two subroutines, we just need to generate Mersenne numbers in the form (2 ** p) -1, with p prime, check whether said Merssenne number is prime, and, if such is the case, generate the ((2 ** p) -1) * 2 ** (p - 1) perfect number associated with it.

The following program finds the first eitght perfect numbers in about 0.3 second:

#!/usr/bin/perl
use strict;
use warnings;
use feature "say";
use constant largest_num => 100000;

sub find_primes {
    my $num = 5;
    my @primes = (2, 3, 5);
    while (1) {
        $num += 2;     # check only odd numbers
        last if $num > largest_num;
        my $limit = int $num ** 0.5;
        my $num_is_prime = 1;
        for my $prime (@primes) {
            last if $prime > $limit;
            if ($num % $prime == 0) {
                # $num evenly divided by $prime, exit the for loop
                $num_is_prime = 0;
                last;
            }
        }
        push @primes, $num if $num_is_prime; #  Found a new prime, add it to the array of primes
    }
    return @primes;
}

my @prime_numbers = find_primes;  
# print "@prime_numbers \n";

sub is_prime {
    my $num = shift;
    my $limit = 1 + int $num ** 0.5;
    for my $p (@prime_numbers) {
        return 1 if $p > $limit;
        return 0 if $num % $p == 0;
    }
    warn "Something got wrong (primes list too small)\n";
    return 0; # If we've reach this point, then our list of primes is 
              # too small, we don't now if the argument is prime, return
              # a false value to be on the safe side of things
}

my @perfect_nums;
for my $prime (@prime_numbers) {
    my $mersenne = 2 ** $prime - 1;
    if (is_prime $mersenne) {
        say "$prime $mersenne ", $mersenne * (2 ** ($prime - 1) );
        push @perfect_nums, $mersenne * (2 ** ($prime - 1) );
        last if scalar @perfect_nums >= 8;
    }
}

say "Perfect numbers: @perfect_nums";

The program displays the following information:

$ time perl perfect.pl
2 3 6
3 7 28
5 31 496
7 127 8128
13 8191 33550336
17 131071 8589869056
19 524287 137438691328
31 2147483647 2305843008139952128
Perfect numbers: 6 28 496 8128 33550336 8589869056 137438691328 2305843008139952128

real    0m0.312s
user    0m0.171s
sys     0m0.077s

The first lines display the p number used for generating the Mersenne number, the Mersenne number itself, and the perfect number corresponding to it. We can see that we checked Mersenne numbers until prime p equal to 31, that Mersenne numbers for primes 11, 23 and 29 have not been retained (they are not prime), and that the eighth perfect number has 19 digits. And the program ran in less than a third of a second.

I have the feeling that it will probably not possible to go further, as the next perfect number has 37 digits, which is probably beyond the integer precision in Perl. We could do it with the Math::BigInt or bigint module (but we might encounter again some long run times). We will try to do that in Perl 6.

Perfect Numbers in Perl 6

The First Four Perfect Numbers in Perl 6

Even though we know by now that the brute force approach will easily yield the first four perfect numbers, but not the fifth one, let's implement our old 1991 CS course assignment in Perl 6, if only for the sake of comparing with Perl 5.

use v6;

sub divisors (Int $num) {
    return 1, | grep { $num %% $_ }, 2 .. (1 + ($num / 2).Int);
}

for 2..10000 -> $num {
    my @divisors = divisors $num;
    say "$num: ", @divisors if $num == [+] @divisors;
}

which produces the following output:

6: [1 2 3]
28: [1 2 4 7 14]
496: [1 2 4 8 16 31 62 124 248]
8128: [1 2 4 8 16 32 64 127 254 508 1016 2032 4064]

Notice how much shorter the P6 code is compared to P5. In fact, the divisors subroutine is now so short that we could inline it in the main code:

use v6;

for 2..10000 -> $num {
    my @divisors = 1, | grep { $num %% $_ }, 2 .. (1 + ($num / 2).Int);
    say "$num: ", @divisors if $num == [+] @divisors;
}

and make it even shorter. Well, if we just want the perfect numbers without willing to print the list of divisors, we could even get rid of the @divisors temporary array and make it short enough for a one-liner:

perl6 -e 'for 2..10000 -> $num { $num.say if $num == [+] (1, | grep { $num %% $_ }, 2 .. (1 + ($num / 2).Int))};'

So, the code is much shorter, but I regret to say that both these Perl 6 versions are also more than 10 times slower that the Perl 5 counterpart.

A Better Algorithm to Speed up Things

Of course, a better algorithm solves this. For example, we may reuse the idea of the triangular numbers:

use v6;

sub divisors (Int $num) {
    return 1, | grep { $num %% $_ }, 2 .. (1 + ($num / 2).Int);
}
my $triangular-num = 1;
for 2..200 -> $num {
    $triangular-num += $num;
    my $sum-div = [+] divisors $triangular-num;
    say " $triangular-num : @divisors[]" if $sum-div == $triangular-num;
    last if $triangular-num > 10000;
}
say now -  INIT now;

Now, finding the first four perfect numbers takes about a third of a second:

~ perl6 triangular_perfect.p6
6
28
496
8128
0.3385362

So, perhaps, rather than complaining about Perl 6 being slow as I just did above, I should have started with a better algorithm.

Using Lazy Infinite Sequences

One of the difficulties in the code above is to determine the required range for $num. The math is not too difficult to compute that the upper bound should be close the square root of twice the ceiling for perfect numbers (i.e. about 141 for a ceiling of 10,000), but I used an upper bound of 200 to be on the safe side of things. The best would be not to have to compute that upper bound. Here come to the rescue lazy infinite lists. For example, we can generate an infinite list @nums of consecutive integers for $num, and Perl 6 will compute them as and when until it reaches the limit for $triangular-num:

use v6;

sub divisors (Int $num) {
    return 1, | grep { $num %% $_ }, 2 .. (1 + ($num / 2).Int);
}
my $triangular-num = 1;
my @nums = 2 ... *;   # Infinite sequence
for @nums -> $num {
    $triangular-num += $num;
    my $sum-div = [+] divisors $triangular-num;
    say $triangular-num if $sum-div == $triangular-num;
    last if $triangular-num > 10000;
}
say now -  INIT now;

It is more concise and probably more idiomatic to generate directly an infinite sequence of triangular numbers:

use v6;

sub divisors (Int $num) {
    return 1, | grep { $num %% $_ }, 2 .. (1 + ($num / 2).Int);
}   
my $num = 1;
my @triangular-numbers = 1,  * + ++$num ... *;
for @triangular-numbers -> $triangular-num {
    last if $triangular-num > 10000;
    say $triangular-num if $triangular-num == [+] divisors $triangular-num;
}
say now -  INIT now;

The key code line here is where the @triangular-numbers sequence is defined. It is an infinite sequence using an explicit generator: each value is created by adding the previous value and a number, $num, which is itself incremented at each step through the process.

Rather than building an infinite list, we can build a list of the triangular numbers less than 10,000 (thereby making the last statement in the loop unnecessary) by adding a code block saying where the sequence should stop:

use v6;

sub divisors (Int $num) {
    return 1, | grep { $num %% $_ }, 2 .. (1 + ($num / 2).Int);
}   
my $num = 1;
my @triangular-numbers = 1,  * + ++$num ...^ * > 10000 ;
.say if $_ == [+] divisors $_ for @triangular-numbers;
say now - INIT now;

This program displays the same four perfect numbers:

perl6 triangular_perfect.p6
1
6
28
496
8128
0.2695492

You might object that the requirement was to display the first five perfect numbers, and we haven't done that yet in Perl 6. Yes, indeed. Even if our program now runs pretty fast, we know from the experience in Perl 5 that this won't be fast enough for the fifth perfect number. So, just as for P5, we'll use the Mersenne Primes in P6.

Mersenne Prime Numbers in Perl 6

In our Perl 5 program using Mersenne numbers, we spent quite a bit of energy trying to make an efficient is_prime subroutine.

We don't need that in Perl 6, since it has a built-in function, is-prime, which uses the very fast Miller-Rabin algorithm to find out whether a number is prime. To tell the truth, the Miller-Rabin algorithm is a probabilistic test, which may in theory report that a composite number is prime, but the probability of occurrence of such an event is so abysmally low that we can just dismiss it for practical purposes.

The other advantage of P6 is its built-in ability to do integer calculations with arbitrary precision.

This a Perl 6 program displaying the first 12 perfect numbers in about half a tenth of a second:

use v6;

for grep { is-prime $_ }, 0..300 -> $prime {
    my $mersenne =  2 ** $prime - 1;
    say $mersenne * 2 ** ($prime - 1) if is-prime $mersenne;
}
say "time taken ", now - INIT now;

This displays the following:

~ perl6 mersenne_perfect.p6
6
28
496
8128
33550336
8589869056
137438691328
2305843008139952128
2658455991569831744654692615953842176
191561942608236107294793378084303638130997321548169216
13164036458569648337239753460458722910223472318386943117783728128
14474011154664524427946373126085988481573677491474835889066354349131199152128
time taken: 0.0533907

The last number displayed has 77 digits, and calculating these 12 perfect numbers took less than 0.06 second.

Computing the first 15 perfect numbers took 0.64 second on my old inefficient laptop, and the fifteenth perfect number (built with prime number 1279) has 770 digits.

Centering Text

Write a function, 'center', whose argument is a list of strings, which will be lines of text. The function should insert spaces at the beginning of the lines of text so that if they were printed, the text would be centered, and return the modified lines.

For example,

center("This", "is", "a test of the", "center function");

should return the list:

"     This", "      is", " a test of the", "center function"

because if these lines were printed, they would look like:

       This
        is
   a test of the
  center function

The Center Subroutine in Perl 5

This challenge is fairly simple. This is a possible implementation in Perl 5:

#!/usr/bin/perl
use strict;
use warnings;
use feature qw / say /;

sub center {
    my $max_size = 0;
    for my $str (@_) {
        my $length = length $str;
        $max_size = $length if $length > $max_size;
    }
    map { " " x (int ($max_size - length) / 2) . $_} @_;
}   

my @strings = ("This", "is", "a test of the", "center function");
say for center @strings;

This displays the following output:

$ perl center.pl
     This
      is
 a test of the
center function

We could also use the max function of the core List::Util module to make the center subroutine simpler:

#!/usr/bin/perl
use strict;
use warnings;
use feature qw / say /;
use List::Util qw / max /;

sub center {
    my $max_size =  max map length, @_;
    map { " " x (int ($max_size - length) / 2) . $_} @_;
}

my @strings = ("This", "is", "a test of the", "center function");
say for center @strings;

This script displays the same output as before.

Note that Perl 5 has the format functionality which could be used for centering pieces of text. I have barely used it, perhaps 2 or 3 times in the last 15 years or so, and I believe that it has become somewhat obsolete. Besides, I don't think using it here would really satisfy the requirement (and it would be overkill).

The Center Subroutine in Perl 6

Perl 6 has a built-in max function that we can use to compute the size of the longest string much in the same way as in our second P5 implementation:

use v6;

sub center (@str) {
    my $max-size = max map { .chars }, @str;
    return map { " " x (($max-size - .chars) / 2).Int ~ $_}, @str;
}

my @strings = "This", "is", "a test of the", "center function";
.say for center @strings;

This script displays the following output:

~ perl6 center.pl6
     This
      is
 a test of the
center function

Although I personally prefer the functional style when using map, we could also do it with a chained method-invocation syntax:

sub center (@str) {
    my $max-size = @str.map({.chars}).max;
    @str.map({" " x (($max-size - .chars) / 2).Int ~ $_})
}

my @strings = "This", "is", "a test of the", "center function";
.say for center @strings;

Wrapping up

There was a third challenge this week: Create a simple Perl client for @mailgun API. I'm not sure I even understand the requirement, but I guess this probably has to do with creating a Web interface. Since I know next to nothing about this kind of topic, I won't undertake anything on that subject and even less blog about it

The next week Perl Weekly Challenge is due to start very soon. If you're interested in participating in this challenge, please check https://perlweeklychallenge.org/ and make sure you answer the challenge before 6 p.m. BST (British summer time) on next Sunday, May 26. And, please, also spread the word about the Perl Weekly Challenge if you can.