From Babylon to Pascal with Raku

by Arne Sommer

From Babylon to Pascal with Raku

[6] Published 14. April 2019.

Perl 6 → Raku

This article has been moved from «perl6.eu» and updated to reflect the language rename in 2019.

This is my response to the Perl Weekly Challenge #3.

Challenge #3.1

Create a script to generate 5-smooth numbers, whose prime divisors are less or equal to 5. They are also called Hamming/Regular/Ugly numbers. For more information, please check this wikipedia.

The task is to multiply the prime numbers 1, 2, 3 and 5 (and any amount of them) to get a list of numbers. As a mathematical function, the solution looks like 2a * 3b * 5c, where «a», «b» and «c» are integers from 0 and upwards.

The first five values are:

Value   Expression
1 20 * 30 * 50     (1 * 1 * 1)
2 21 * 30 * 50(2 * 1 * 1)
3 20 * 31 * 50(1 * 3 * 1)
4 23 * 30 * 50(1 * 3 * 1)
5 20 * 30 * 51(1 * 1 * 5)

First Try

There are two possible approaches, specifying the upper limit for a, b and c, or giving a limit for the regular number. It is easier (for me as a programmer) to go for the first one, so I do that. We can let the user decide the upper limit for the values a, b and c. For simplicity's sake we use the same value for all of them.

Then we just set up three loops inside each other, one for each number (2, 3 and 5). The values will not come in the right order, so we can either put them in a list and sort it before printing (after removing any duplicate values, if any), or we can use as hash. I have chosen to use a hash, or rather the «SetHash» variant which only accepts «True» or «False» as values. (The value «False» actually means that the value isn't present in the «SetHash».)

See docs.raku.org/type/SetHash for more information about «SetHash».

File: regular
sub MAIN (Int $limit where $limit > 0)
{
  my SetHash $solution = SetHash;

  for 0 .. $limit -> $a
  {
    for 0 .. $limit -> $b
    {
      for 0 .. $limit -> $c
      {
         $solution{ 2 ** $a * 3 ** $b * 5 ** $c} = True;
      }
    }
  }
  
  say $solution.keys.sort.join(" ");
}

Running it (with newlines added to make it fit the page):

$ raku regular 1
1 2 3 5 6 10 15 30

$ raku regular 2
1 2 3 4 5 6 9 10 12 15 18 20 25 30 36 45 50 60 75 90 100 150 180 225 300 450 \
900

$ raku regular 3
1 2 3 4 5 6 8 9 10 12 15 18 20 24 25 27 30 36 40 45 50 54 60 72 75 90 100 \
108 120 125 135 150 180 200 216 225 250 270 300 360 375 450 500 540 600 675 \
750 900 1000 1080 1125 1350 1500 1800 2250 2700 3000 3375 4500 5400 6750 \
9000 13500 27000

$ raku regular 4
1 2 3 4 5 6 8 9 10 12 15 16 18 ...

The number 16 appears in the last one, and not when we used «3» as argument:

Value   Expression
1624 * 30 * 50     (16 * 1 * 1)

So this approach isn't working.

Second Try

This time we'll take the second approach, where we specify the upper limit for the regular number. The loop order is swapped, so that the lowest prime (2) is in the inner most loop.

File: regular2
sub MAIN (Int $limit where $limit > 0)  # [1]
{
  my SetHash $solution = SetHash;

  for 0 .. Inf -> $c                    # [2]
  {
    last if 5 ** $c > $limit;           # [3]
    
    for 0 .. Inf -> $b                  # [4]
    {
      last if 3 ** $b > $limit;         # [5]

      for 0 .. Inf -> $a                            # [6]
      {
         my $value =  2 ** $a * 3 ** $b * 5 ** $c;
         last if $value > $limit;                   # [7]
         $solution{$value} = True;                  # [8]
      }
    }
  }
   
  say $solution.keys.sort.join(" ");                # [9]
}

[1] [3] [5] [7] The limit argument is now a true limit, the highest number to consider.

[2] [4] [6] We have swapped the order of the loops, starting with 5 (the highest prime), so that the increase in values by the inner most loop will be as low as possible. They go from 0 to Infinity.

[8] Register the value.

[9] Display the list, sorted and separated by spaces.

Running it:

$ raku regular2 1
1

$ raku regular2 2
1 2

$ raku regular2 3
1 2 3

$ raku regular2 4
1 2 3 4

$ raku regular2 5
1 2 3 4 5

$ raku regular2 6
1 2 3 4 5 6

$ raku regular2 7
1 2 3 4 5 6

$ raku regular2 8
1 2 3 4 5 6 8

$ raku regular2 9
1 2 3 4 5 6 8 9

$ raku regular2 16
1 2 3 4 5 6 8 9 10 12 15 16

$ raku regular2 22
1 2 3 4 5 6 8 9 10 12 15 16 18 20

$ raku regular2 100
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 40 45 48 50 54 60 64 \
72 75 80 81 90 96 100

A walk-through, with 16 as limit:

$a   $b   $c   ExpressionValue   Comment
00020 * 30 * 50   (1 * 1 * 1) 1OK
10021 * 30 * 50(2 * 1 * 1) 2OK
20022 * 30 * 50(4 * 1 * 1) 4OK
30023 * 30 * 50(8 * 1 * 1) 8OK
40024 * 30 * 50(16 * 1 * 1) 16OK
50025 * 30 * 50(32 * 1 * 1)   32Not OK, [7]
01020 * 31 * 50(1 * 3 * 1) 3OK
11021 * 31 * 50(2 * 3 * 1) 6OK
21022 * 31 * 50(4 * 3 * 1) 12OK
31023 * 31 * 50 (8 * 3 * 1) 24Not OK, [7]
02020 * 32 * 50(1 * 9 * 1) 9OK
12021 * 32 * 50(2 * 9 * 1) 18Not OK, [7]
03033(27)27Not OK, [5]
00120 * 30 * 51(1 * 1 * 5) 5OK
10121 * 30 * 51(2 * 1 * 5) 10OK
20122 * 30 * 51(4 * 1 * 5) 20Not OK, [7]
01120 * 31 * 51(1 * 3 * 5) 15OK
11121 * 31 * 51(2 * 3 * 5) 30Not OK, [7]
02120 * 32 * 51(1 * 9 * 5) 45Not OK, [7]
03133(27) 27Not OK, [5]
00252(25)25Not OK, [3]

We can modify the program to print the walk-through when we give the «--verbose» argument:

File: regular2v:
sub MAIN (Int $limit where $limit > 0, :$verbose = False)
{
  my SetHash $solution = SetHash;

  for 0 .. Inf -> $c
  {
    if 5 ** $c > $limit
    {
      say "0 0 $c -> 5 ** $c -> { 5 ** $c } -> NOT OK" if $verbose;
      last;
    }
    
    for 0 .. Inf -> $b
    {
      if 3 ** $b > $limit
      {
        say "0 $b $c -> 3 ** $b -> { 3 ** $b} -> NOT OK" if $verbose;
        last;
      }

      for 0 .. Inf -> $a
      {
         my $value = 2 ** $a * 3 ** $b * 5 ** $c;
         if $value > $limit
	 {
           say "$a $b $c = $value -> NOT OK" if $verbose; 
	   last;
	 }
         say "$a $b $c = $value -> OK" if $verbose; 
         $solution{$value} = True;
      }
    }
  }
  print "\n" if $verbose;
  
  say $solution.keys.sort.join(" ");
}

Running it:

$ raku regular2v --verbose 100
0 0 0 = 1 -> OK
1 0 0 = 2 -> OK
2 0 0 = 4 -> OK
3 0 0 = 8 -> OK
4 0 0 = 16 -> OK
5 0 0 = 32 -> OK
6 0 0 = 64 -> OK
7 0 0 = 128 -> NOT OK
0 1 0 = 3 -> OK
1 1 0 = 6 -> OK
2 1 0 = 12 -> OK
3 1 0 = 24 -> OK
4 1 0 = 48 -> OK
5 1 0 = 96 -> OK
6 1 0 = 192 -> NOT OK
0 2 0 = 9 -> OK
1 2 0 = 18 -> OK
2 2 0 = 36 -> OK
3 2 0 = 72 -> OK
4 2 0 = 144 -> NOT OK
0 3 0 = 27 -> OK
1 3 0 = 54 -> OK
2 3 0 = 108 -> NOT OK
0 4 0 = 81 -> OK
1 4 0 = 162 -> NOT OK
0 5 0 -> 3 ** 5 -> 243 -> NOT OK
0 0 1 = 5 -> OK
1 0 1 = 10 -> OK
2 0 1 = 20 -> OK
3 0 1 = 40 -> OK
4 0 1 = 80 -> OK
5 0 1 = 160 -> NOT OK
0 1 1 = 15 -> OK
1 1 1 = 30 -> OK
2 1 1 = 60 -> OK
3 1 1 = 120 -> NOT OK
0 2 1 = 45 -> OK
1 2 1 = 90 -> OK
2 2 1 = 180 -> NOT OK
0 3 1 = 135 -> NOT OK
0 4 1 = 405 -> NOT OK
0 5 1 -> 3 ** 5 -> 243 -> NOT OK
0 0 2 = 25 -> OK
1 0 2 = 50 -> OK
2 0 2 = 100 -> OK
3 0 2 = 200 -> NOT OK
0 1 2 = 75 -> OK
1 1 2 = 150 -> NOT OK
0 2 2 = 225 -> NOT OK
0 3 2 = 675 -> NOT OK
0 4 2 = 2025 -> NOT OK
0 5 2 -> 3 ** 5 -> 243 -> NOT OK
0 0 3 -> 5 ** 3 -> 125 -> NOT OK

1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 40 45 48 50 54 60 64 \
72 75 80 81 90 96 100

Challenge #3.2

Create a script that generates Pascal Triangle. Accept number of rows from the command line. The Pascal Triangle should have at least 3 rows. For more information about Pascal Triangle, check this wikipedia page.

I'll jump right in and show the triangle as printed by the program, as it makes it easier to discuss the algorithm. Here it is with 10 lines, as specified as argument:

$ raku pascal 10
                             1
                          1     1
                       1     2     1
                    1     3     3     1
                 1     4     6     4     1
              1     5    10    10     5     1
           1     6    15    20    15     6     1
        1     7    21    35    35    21     7     1
     1     8    28    56    70    56    28     8     1
  1     9    36    84   126   126    84    36     9     1

Some observations:

  • The first row has 1 element, and the value is 1
  • The second row has 2 elements, the third has three elements, and so on. (The number of elements is the same as the row number.)
  • Any given value in the triangle is the sum of its to parents, the value above it to the left, and the one above it to the right. A missing parent value (outside of the triangle) is the same as zero. This rule applies to every value, except the one in the first row.

I have chosen to use a two dimentional array to store the values. The first dimension is the row, and the second is the index.

The value of any given cell (c) in a row (r) can be expressed like this in quasi-quode:

value[r][c] = value[r-1][c-1] + value[r-1][c+1] 

As negative indices in an array doesn't work, I chose to start with 11 (for the top node). That means that it is possible to compute a pyramid with 12 rows. If we try with 13 we get a run time error. This is esily fixed by using the specified number (actually +1) as the start index.

The index pyramid looks like this:

$ raku pascal-verbose --graphic  10
 1:                                      11    
 2:                                  10      12    
 3:                               9      11      13    
 4:                           8      10      12      14    
 5:                       7       9      11      13      15    
 6:                   6       8      10      12      14      16    
 7:               5       7       9      11      13      15      17    
 8:           4       6       8      10      12      14      16      18    
 9:       3       5       7       9      11      13      15      17      19    
10:   2       4       6       8      10      12      14      16      18      20

This is the code printing it:

File: pascal-verbose (partial)
my $top = 10;

sub MAIN (Int $size where $size > 0, :$graphic = False)
{
  $top = $size if $top < $size;               # [1]

  if $graphic
  {
    for 1 .. $size -> $level                  # [2]
    {
       print $level.fmt('%2d'), ":";          # [3]
       print " " x ($size - $level) * 4;      # [4]

       for 1 .. $level -> $current            # [5]
       {
          print ($top - $level + $current * 2).fmt('%4d'), "    ";  # [6]
       }
       print "\n";                            # [7]
    }
  }
}

[1] Avoid negative indices.

[2] One iteration for each line in the pyramid, from the top.

[3] Print the line number. fmt('%2d') means that the number will be printed with at least two characters. If the number only has one digit, it is prefixed by a space.

[4] Add indentation before the first value.

[5] One iteration for each value on the line. The count is the same as the line number.

[6] Print the ID, using four characters (adding a leading space if needed), and four spaces as padding before the next value, also in [6].

[7] Add a newline before starting with the next row, in [2].

The program (without the code printing the index pyramid):

File: pascal
my $top = 10;
 
sub MAIN (Int $size where $size > 0)
{
  $top = $size if $top < $size;         # [1]

  my @values;
  @values[1][$top +1] = 1;              # [2]

  for 2 .. $size -> $level              # [3]
  {    
     for 1 .. $level -> $current        # [4]
     {
        my $id = $top - $level + $current * 2;                     # [5] 
        @values[$level][$id] = (@values[$level -1][$id -1] // 0)   # [6]
	                     + (@values[$level -1][$id +1] // 0);  # [7]
     }
  }

  for 1 .. $size -> $level              # [8]
  {
     print " " x ($size - $level) * 3;  # [9]
     say @values[$level; *].grep( *.defined ).map( *.fmt('%3d') ).join("   ");
         # [10] ##########  # [11] #########  # [12] ###########  # [13] ####
  }
}

[1] Prevent a negative index (which would have terminated the program) from happening.

[2] The topmost value is 1, and we put that in manually.

[3] The first row has been set up, so we iterate over the second, third, and so on, rows.

[4] All the values in the row.

[5] Computing the ID, or index. Look at the index pyramid, and try to verify the values.

[6] The value of the left parent.

[7] The value of the right parent.

[8] Finally we'll print the pyramid, one row at the time.

[9] Starting with padding.

[10] This gives a one dimentional list with all the values for this row.

[11] Every element in the list with lower indices than the one we set exist, but the values are undefined:

> my @c;
> @c[1000] = 1;
> say @c.elems;  # -> 1001

So we use grep to only select the ones that have a defined value.

[12] We use map to print them as three characters (leading space added as needed).

[13] And finally joining the values, with four spaces between them.

Note that using arrays like this is a waste of memory (and some time as well), as we set aside space for a lot of unused positions - especially for large triangles. (You can try running the program with 1000 as argument. This gives a pyramid with 1000 rows that is impossible to read on a screen. Except for that, it works fine, and takes about 45 seconds to run on my computer.)

// is the «defined-or» operator. The value on the right side is used if the left hand side is undefined. In this case a normal || «or» could have been used, but «defined-or» is logically the correct choice. See docs.raku.org/routine/$SOLIDUS$SOLIDUS for more information.

Optimized Pascal

Take a look at the Pascal Triangle again:

$ raku pascal 9
                          1
                       1     1
                    1     2     1
                 1     3     3     1
              1     4     6     4     1
           1     5    10    10     5     1
        1     6    15    20    15     6     1
     1     7    21    35    35    21     7     1
  1     8    28    56    70    56    28     8     1

You may have noticed that it is symmetrical (on both sides of the middle column; with the values 1, 2, 6, 20 and 70)? This means that we compute a lot of values 2 times. That is a waste.

We can fix that. We give the top node the index zero this time. The new index pyramid looks like this:

$ raku pascal2-verbose  --graphic  10
 1:                                       0    
 2:                                  -1       1    
 3:                              -2       0       2    
 4:                          -3      -1       1       3    
 5:                      -4      -2       0       2       4    
 6:                  -5      -3      -1       1       3       5    
 7:              -6      -4      -2       0       2       4       6    
 8:          -7      -5      -3      -1       1       3       5       7    
 9:      -8      -6      -4      -2       0       2       4       6       8
10:  -9      -7      -5      -3      -1       1       3       5       7       9

Whenever the program encounters a negative index, it should simply use the positive one - as that gives the correct value.

The changes to the computation is minimal (and the index pyramid code doesn't need any changes).

File: pascal2-verbose (partial)
my $top = -1;
...
# $top = $size if $top < $size;        # Remove this line
...
my $id = $top - $level + $current * 2; # Existing line
next if $id < 0;                       # This line is new
@values[$level][$id] = (@values[$level -1][abs($id -1)] // 0)   # Note the "abs"
                     + (@values[$level -1][$id +1]      // 0);

It is probably possible to rewrite the inner loop so that we could avoid the next statement. But this way of doing it is probably easier to understand than a complex start value in the loop.

Note the abs on the index, to remove a potential negative sign.

The challenge actually said that «the Pascal Triangle should have at least 3 rows». We can enforce that by changing the signature:

#| size is an integer from 3 and upwards.
sub MAIN (Int $size where $size >= 3)

The special comment gives us a much nicer usage message if we run the program without arguments (and «pascal2» is the new version):

$ raku pascal
Usage:
  pascal <size>

$ raku pascal2
Usage:
  pascal2  -- size is an integer from 3 and upwards.

If we use a custom type for the argument, we get the usage message if the argument is of the wrong type (or value) as well:

subset Int3 of Int where * >= 3;

#| size is an integer from 3 and upwards.
sub MAIN (Int3 $size)

See docs.raku.org/language/typesystem#index-entry-subset-subset for more information about «subset».

$ raku pascal A
Type check failed in binding to parameter '<anon>'; expected Any but got Mu (Mu)
  in block >unit> at pascal line 3

$ raku pascal2 A
Usage:
  pascal2 <size> -- size is an integer from 3 and upwards.

The code printing the Triangle is more complicated this time:

File: pascal2-verbose (partial)
  for 1 .. $size -> $level
  {
     print " " x ($size - $level) * 3;                    # [1]
     my @partial = @values[$level; *].grep( *.defined );  # [2]

     my @row = @partial.reverse;                          # [3]

     @partial.pop if @values[$level;0].defined;           # [4]

     @row.append(@partial);                               # [5]

     say @row.map( *.fmt('%3d') ).join("   ");            # [6]
  }

[1] Indentation, as the last time.

[2] The defined values for the row. They are from index 0 and to the right (in the index pyramid)

[3] We reverse the values from [2], giving the left side of the pyramid (including the middle value; with index 0 - if any).

[4] If the row has a middle value (index zero), we remove it from the list we got in [2].

[5] We add the list from step [4] to the row, giving a full row.

[6] Printing the full row of values.

This version has no unused cells in the array, so there is no waste of memory as we had in the first version.

The complete program (without the code printing the index pyramid):

File: pascal2
my $top = -1;

subset Int3 of Int where * >= 3;

#| size is an integer from 3 and upwards.
sub MAIN (Int3 $size)
{
  my @values;
  @values[1][$top +1] = 1;

  for 2 .. $size -> $level
  {    
     for 1 .. $level -> $current
     {
        my $id = $top - $level + $current * 2;
	next if $id < 0;
	
        @values[$level][$id] = (@values[$level -1][abs($id -1)] // 0)
	                     + (@values[$level -1][$id +1]      // 0);
     }
  }

  for 1 .. $size -> $level
  {
     print " " x ($size - $level) * 3;
     my @partial = @values[$level; *].grep( *.defined );

     my @row = @partial.reverse;

     @partial.pop if @values[$level;0].defined;

     @row.append(@partial);

     say @row.map( *.fmt('%3d') ).join("   ");
  }
}

This version with 1000 as argument takes about 43 seconds, compared with 45 seconds for the first version. The speedup isn't very impressive, but most of the time (in both versions) is used by the printing of the values. (If we remove the printing code, the speedup gain is about 33%.)