Doing math with Perl 6

Different mathematical paradigms and how they are implemented in this language

Sets

Perl 6 includes the Set data type, as well as support for most set operations. Union and intersection are not only native operations, they use their natural symbols, ∩ and ∪. For instance, this code would check the fundamental laws of the arithmetic of sets for a limited number of sets:

my @arbitrary-numbers = ^100;
my \U = @arbitrary-numbers.Set;
 
my @sets;
 
@sets.push: Set.new@arbitrary-numbers.pick@arbitrary-numbers.elems.rand)) for @arbitrary-numbers;
 
my (@union@intersection);
 
for @sets -> $set {
    @union.push: $set  $set === $set;
    @intersection.push: $set  $set === $set;
}
 
say "Idempotent union is "so @union.all;
# OUTPUT: «Idempotent union is True» 
say "Idempotent intersection is "so @intersection.all;
# OUTPUT: «Idempotent intersection is True» 
my (@universe@empty-set@id-universe@id-empty);
 
for @sets -> \A {
    @universe.push: A  U === U;
    @id-universe.push: A  U === A;
    @empty-set.push: A  ∅ === ∅;
    @id-empty.push: A  ∅ === A;
}
 
say "Universe dominates "so @universe.all;    # OUTPUT: «Universe dominates True» 
say "Empty set dominates "so @empty-set.all;  # OUTPUT: «Empty set dominates True» 
 
say "Identity with U "so @id-universe.all;    # OUTPUT: «Identity with U True» 
say "Identity with ∅ "so @id-empty.all;       # OUTPUT: «Identity with ∅ True» 

In this code, which uses the empty set which is already defined by Perl 6, not only do we check if the equalities in the algebra of sets hold, we also use, via sigilless variables and the Unicode form of the set operators, expressions that are as close as possible to the original form; A ∪ U === U, for example, except for the use of the value identity operator <===> is very close to the actual mathematical expression in the Wikipedia entry.

We can even test De Morgan's law, as in the code below:

my @alphabet = 'a'..'z';
my \U = @alphabet.Set;
sub postfix:<>(Set $a{ U  $a }
my @sets;
@sets.push: Set.new@alphabet.pick@alphabet.elems.rand)) for @alphabet;
my ($de-Morgan1,$de-Morgan2= (True,True);
for @sets X @sets -> (\A, \B){
    $de-Morgan1 &&= (A  B)⁻  === A⁻  B⁻;
    $de-Morgan2 &&= (A  B)⁻  === A⁻  B⁻;
}
say "1st De Morgan is "$de-Morgan1;
say "2nd De Morgan is "$de-Morgan2;

We declare as the complement operation, which computes the symmetrical difference ⊖ between the Universal set U and our set. Once that is declared, it is relatively easy to express operations such as the complementary of the union of A and B, (A ∪ B)⁻, with a notation that is very close to the original mathematical notation.

Arithmetic

Perl 6 can do arithmetic using different data types. Num, Rat and Complex can all operate as a field under the operations of addition, subtraction, multiplication and division (technically, it should be noted that data types dealing with floating point number representations are not a field in the mathematical sense due to the inherent imprecisions of their arithmetic. However, they constitute an approximate enough, computer friendly version of such mathematical objects for most of the cases). The equivalent mathematical fields are:

Perl 6 class Field
Rat
Num
Complex

The Ints or ℤ, as they're usually called in mathematics, are not a mathematical field but rather a ring, since they are not closed under multiplicative inverses. However, if the integer division div is used, their operations will always yield other integers; if / is used, on the other hand, in general the result will be a Rat.

Besides, Int can do infinite-precision arithmetic (or at least infinite as memory allows; Numeric overflow can still occur), without falling back to Num if the number is too big:

my @powers = 22 ** * ... Infsay @powers[4].chars# OUTPUT: «19729␤» 

Also strictly speaking, the Rational class that behaves like a mathematical field is FatRat. For efficiency reasons, operating with Rats will fall back to Num when the numbers are big enough or when there is a big difference between numerator and denominator. FatRat can work with arbitrary precision, the same as the default Int class.

Some modules in the ecosystem can work with additional data types mathematically:

Numbers are duck-typed automatically to the numeric class they actually represent:

.^name.say for (4, ⅗, 1e-93+.1i); # OUTPUT: «Int␤Rat␤Num␤Complex␤» 

Which makes also arithmetic operations the most adequate for the particular type

say .33-.22-.11 == 0# OUTPUT: «True␤» 

In this case, all numbers are interpreted as Rats, which makes the operation exact. In general, most languages would interpret them as floating point numbers,

say .33.Num -.22.Num - .11.Num# OUTPUT: «1.3877787807814457e-17␤» 

For cases such as this, Perl 6 also includes an approximately equal operator,

say .33.Num -.22.Num - .11.Num  0# OUTPUT: «True␤» 

Sequences

A sequence is an enumerated collection of objects in which repetitions are allowed, and also a first-class data type in Perl 6 called Seq. Seq is able to represent infinite sequences, like the natural numbers:

my \𝕟 = 1,2 … ∞;
say 𝕟[3];# OUTPUT: «4␤» 

Infinite sequences use ∞, Inf or * (Whatever) as terminator. is the list generator, which in fact can understand arithmetic and geometric progression sequences as long as you insert the first numbers:

say 1,5,9 … * > 100;
# OUTPUT: «(1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77 81 85 89 93 97 101)␤» 
say 1,3,9 … * > 337# OUTPUT: «(1 3 9 27 81 243 729)␤» 

The first sequence will be terminated when the generated number is bigger than 100; the second sequence, which is a geometric progression, when it is bigger than 337.

The fact that an arbitrary generator can be used makes easy to generate sequences such as Fibonacci's:

say 1,1* + * … * > 50;#  OUTPUT: «(1 1 2 3 5 8 13 21 34 55)␤» 

We can, in fact, compute the approximation to the golden ratio this way:

my @phis = (2.FatRat1 + 1 / * ... *);
my @otherphi = (1 - @phis[200], 1 + 1 / * ... *);
say @otherphi[^10|(2030 ... 100)];# OUTPUT: 
# «((-0.61803398874989484820458683436563811772030918 
# -0.61803398874989484820458683436563811772030918 
# -0.61803398874989484820458683436563811772030918 
# -0.61803398874989484820458683436563811772030918 
# -0.61803398874989484820458683436563811772030918 
# -0.618033…» 

The Math::Sequences module includes many mathematical sequences, already defined for you. It defines many sequences from the encyclopedia, some of them with their original name, such as ℤ.

Some set operators also operate on sequences, and they can be used to find out if an object is part of it:

say 876  (7,14 … * > 1000) ; # OUTPUT: «False␤» 

In this particular case, we can find out if 876 is a multiple of 7 straight away, but the same principle holds for other sequences using complicated generators. And we can use set inclusion operators too:

say (55,89).Set  (1,1* + * … * > 200); # OUTPUT: «True␤» 

although it does not take into account if it is effectively a subsequence, just the presence of the two elements here. Sets have no order, and even if you don't explicitly cast the subsequence into a Set or explicitly cast it into a Seq it will be coerced into such for the application of the inclusion operator.

Mathematical constants

Perl 6 includes already a set of mathematical constants as part of the core.

say π# OUTPUT: «3.141592653589793» 
say τ# Equivalent to 2π; OUTPUT: «6.283185307179586» 
say 𝑒; # OUTPUT: «2.718281828459045␤» 

which are available also by their Latin name, e, pi and tau, with the same value (although 𝑒 is not available outside the MoarVM).

The Math::Constants module includes an additional series of physical and mathematical constants such as the previously mentioned golden ratio φ or the Planck's constant ℎ.

Since Perl 6 allows for definition of variables that use Unicode graphemes, and also variable and constant names without any kind of sigil, it is considered a good practice to use the actual mathematical name of concepts to denominate them wherever possible.

Numerical integration of ordinary differential equations

Perl 6 is an amazing programming language, and of course, you can do a lot of cool maths with it. A great amount of work during an applied mathematician's work is to simulate the models they create. For this reason, in every coding language, a numerical integrator is a must-have. Learning about how to do this in Perl 6 can easily be a good way to have a nice goal while learning a new coding language.

Requirements

In Perl 6 there are some modules in the ecosystem that can do the job we want in a very easy way :

For this example we are going to use Math::Model because its useful sintaxis, but remember that this module requires Math::RungeKutta as well. Simply install them with zef before using these examples. In the future it should be a nice idea to add other modules.

Malthus model

Let's start with the 'Hello World' of mathematical Ecology: Malthusian growth model. A Malthusian growth model, sometimes called a simple exponential growth model, is essentially exponential growth based on the idea of the function being proportional to the speed to which the function grows. The equation, then, looks like this:

dx/dt = g*x

x(0) = x_0

Where g is the population growth rate, sometimes called Malthusian parameter.

How can we translate that into Perl 6? Well Math::Model brings some help with a very understandable way to do that--let's see it:

use Math::Model;
 
my $m = Math::Model.new(
    derivatives => {
        velocity     => 'x',
    },
    variables   => {
        velocity           => { $:growth_constant * $:x},
        growth_constant    => { 1 },   # basal growth rate 
    },
    initials    => {
        x       => 3,
    },
    captures    => ('x'),
);
 
$m.integrate(:from(0), :to(8), :min-resolution(0.5));
$m.render-svg('population growth malthus.svg':title('population growth'));

To fully understand what is going on, let me explain it step by step.

Step by step explanation.

At this point our model is set. We need to run the simulation and render a cool plot about our results:

$m.integrate(:from(0), :to(8), :min-resolution(0.5));
$m.render-svg('population growth malthus.svg':title('population growth'));

There, we select our time limits and the resolution. Next, we generate a plot. All understood? Well, let's see our result!

link to the image

Looks great! But to be honest, it is quit unrepresentative. Let's explore other examples from more complex situations!

Logistic model

Resources aren't infinite and our population is not going to grow forever. P-F Verhulst thought the same thing, so he presented the logistic model. This model is a common model of population growth, where the rate of reproduction is proportional to both the existing population and the amount of available resources, all else being equal. It looks like this:

dx/dt = g*x*(1-x/k) x(0)=x_0

where the constant g defines the growth rate and k is the carrying capacity. Modifying the above code we can simulate its behaviour in time:

use Math::Model;
 
my $m = Math::Model.new(
    derivatives => {
        velocity     => 'x',
    },
    variables   => {
        velocity           => { $:growth_constant * $:x - $:growth_constant * $:x * $:x / $:k },
        growth_constant    => { 1 },   # basal growth rate 
        k                  => { 100 }# carrying capacity 
    },
    initials    => {
        x       => 3,
    },
    captures    => ('x'),
);
 
$m.integrate(:from(0), :to(8), :min-resolution(0.5));
$m.render-svg('population growth logistic.svg':title('population growth'));

Let's look at our cool plot: link to the image

As you can see population growths till a maximum.

Strong Allee Effect

Interesting, isn't it? Even if these equations seem basic they are linked to a lot of behaviours in our world, like tumor growth. But, before end, let me show you a curious case. Logistic model could be accurate but... What happens when, from a certain threshold, the population size is so small that the survival rate and / or the reproductive rate drops due to the individual's inability to find other ones?

Well, this is an interesting phenomenon described by W.C.Allee, usually known as Allee effect. A simple way to obtain different behaviours associated with this effect is to use the logistic model as a departure, add some terms, and get a cubic growth model like this one:

dx/dt=r*x*(x/a-1)*(1-x/k)

where all the constants are the same as before and A is called critical point.

Our code would be:

use Math::Model;
 
my $m = Math::Model.new(
    derivatives => {
        velocity_x    => 'x',
    },
    variables   => {
        velocity_x           => { $:growth_constant * $:x *($:x/$:a -1)*(1$:x/$:k},
        growth_constant    => { 0.7 },   # basal growth rate 
        k                  => { 100 }# carrying capacity 
        a                  => { 15 },  # critical point 
    },
    initials    => {
        x       => 15,
    },
    captures    => ('x'),
);
 
$m.integrate(:from(0), :to(100), :min-resolution(0.5));
$m.render-svg('population growth alle.svg':title('population growth'));

Try to execute this one yourself to see the different behaviours that arise when you change the initial condition around the critical point!

Weak Allee Effect

Whilst the strong Allee effect is a demographic Allee effect with a critical population size or density, the weak Allee effect is a demographic Allee effect without a critical population size or density, i.e., a population exhibiting a weak Allee effect will possess a reduced per capita growth rate at lower population density or size. However, even at this low population size or density, the population will always exhibit a positive per capita growth rate. This model differs lightly from the strong one, getting a new formula:

dx/dt=r*x*(1-x/k)*(x/a)**n, with n>0

Our code would be:

use Math::Model;
 
my $m = Math::Model.new(
    derivatives => {
        velocity_x    => 'x',
    },
    variables   => {
        velocity_x           => { $:growth_constant * $:x *(1$:x/$:k)*($:x/$:a)**$:n },
        growth_constant    => { 0.7 },   # basal growth rate 
        k                  => { 100 }# carrying capacity 
        a                  => { 15 },  # critical point 
        n                  => { 4  }
    },
    initials    => {
        x       => 15,
    },
    captures    => ('x'),
);
 
$m.integrate(:from(0), :to(100), :min-resolution(0.5));
$m.render-svg('population growth alle.svg':title('population growth'));

Extra info.

Do you like physics? Check the original post by the creator of the modules that have been used, Moritz Lenz: https://perlgeek.de/blog-en/perl-6/physical-modelling.html.