Integer Triangle Traits
A Julia Implementation
Part I

Outline

In three parts, we want to consider how the programming language Julia can be applied in the study of integer sequences, in particular of integer triangles.

The article is written in the hope of making it easier for friends of number sequences (SeqFans) to get started with the Julia language. And conversely, that it might animate one or the other member of the Julia community to write contributions for the OEIS. It can also be read as an introductory tutorial to the registered package "IntegerTriangle".

We choose Julia because it is currently (2020) the language best suited for many tasks in scientific computing and, equally important, is an elegant language that makes coding fun. As one co-creator of Julia puts it: "Julia was designed to help researchers write high-level code in an intuitive syntax and produce code with the speed of production programming languages."

Julia is free and open-source. We use version 1.4 which was released in 2020. This article and the code is available under the Creative Commons Attribution 3.0 license.

Triangles as sequences of sequences

Here we want to look at a simple example, the integer triangles. The conventional view of an integer triangle is that of a matrix split into two parts at the main diagonal of which the lower part is retained. This triangle is indexed by an index 'n' indicating the rows and an index 'k' specifying the position of an element within a row. The typical signature of a function computing the entries is T(n, k).

However, this will not be our view of the underlying data structure. Our generalization of the 1-dimensional sequence is a 1-dimensional sequence of sequences. This is a 'rectangle' (also called an 'array' in the OEIS) if the sequences are infinite, or a 'triangle' if the sequences are finite. The typical signature of a function computing the entries is T(n).

That means when going from the one to the two dimensional case we do not change the type of the container (sequence ↦ matrix) but the type of the elements (integer ↦ sequence).

This approach characterizes the entire project. In short, we always introduce first a function for the case of a sequence, and then extend it to the case of the triangle, usually by broadcasting.

Another immediate consequence is that both types of tables in the OEIS, the "regular table (tabl)", where the length of rows increases by one, and the "irregular table (tabf)" are triangles in our sense.

Definition: An integer triangle is a sequence of sequences whose terms are integers.

Integers, BigInts and fmpz

We will show how this translates into the terminology of Julia. First we have to say what we mean by an integer. It is the type fmpz with the constructor ℤ. These integers are provided by the package Nemo, on which we will rely for many algebraic and number-theoretical calculations. Just regard fmpz as another name for Julia's BigInt. It provides integers which are unrestricted in size (mp stands for 'multiple precision'). The name of the constructor, ZZ, imitates the blackboard ℤ, which in turn derives from the first letter of the German word Zahl. (Don't worry if you see the characters 'ℤ' and '∈' below but only have an old editor: you can write 'ZZ' and 'in' instead.) 

IN:   using Nemo

      [[ℤ(k) for k ∈ 0:n] for n ∈ 0:6]
OUT:  7-element Array{Array{fmpz,1},1}:
      [0]
      [0, 1]
      [0, 1, 2]
      [0, 1, 2, 3]
      [0, 1, 2, 3, 4]
      [0, 1, 2, 3, 4, 5]
      [0, 1, 2, 3, 4, 5, 6]

Julia is helpful and tells us exactly how she understands our input: as an 7-element Array{Array{fmpz,1},1}. So an array of arrays of fmpz integers is in Julia parlance precisely what we will call an integer triangle. The type names are not inviting. We will improve that next.

Types and constructors

Three types and four constructors are the base of our development. The first three definitions are type aliases whose sole purpose is to provide more intuitive names for the Julia types. They are marked as constant.

DEF:  const ℤInt      = fmpz
      const ℤSequence = Array{ℤInt,1}
      const ℤTriangle = Array{ℤSequence,1} 
DEF:  ℤSequence(len) = ℤSequence(undef, len)

      function ℤTriangle(dim::Int; reg=false)
          reg ? ℤSequence.(1:dim) : ℤTriangle(undef,dim)
      end

ℤSequence(len) and ℤTriangle(dim) are the basic constructors for integer sequences and integer triangles. They return uninitialized instances of their respective types. These structures provide the support to which the integer triangles will be applied.

IN:   S = ℤSequence(4)
OUT:  4-element Array{fmpz,1}:
       #undef
       #undef
       #undef
       #undef

 IN:   T = ℤTriangle(4)
 OUT:  4-element Array{Array{fmpz,1},1}:
        #undef
        #undef
        #undef
        #undef

If in the case of triangles the optional argument reg is set true then a regular triangular array with dim rows and uninitialized cells is provided.

IN:   T = ℤTriangle(6, reg=true)
OUT:  6-element Array{Array{fmpz,1},1}:
      [#undef]
      [#undef, #undef]
      [#undef, #undef, #undef]
      [#undef, #undef, #undef, #undef]
      [#undef, #undef, #undef, #undef, #undef]
      [#undef, #undef, #undef, #undef, #undef, #undef]

With the queries isa(S,ℤSequence) and isa(T, ℤTriangle) one can convince oneself that the returned objects are indeed of type ℤSequence respectively ℤTriangle.

Irregular triangles and empty rows

The next two constructors create the objects by a generating function. The given function is applied at all integers in the interval 0:dim-1.

DEF:  ℤSequence(dim, f::Function) = f.(0:dim-1)
      ℤTriangle(dim, f::Function) = f.(0:dim-1)

Equivalently we could have written [f(n) for n ∈ 0:dim-1] for the right hand side of the definitions. We prefer the terse for-free constructor form over the comprehension form. The index-free dot syntax is called broadcasting in Julia. It applies the function f elementwise to the index array 0:dim-1 = [0,1,2,3...,dim-1]. This can be seen as a vectorized form of function evaluation.

To demonstrate such a construction of a ℤTriangle by a generating function we use the prime divisors of an integer. 

      # Return the list of prime divisors of n.

DEF:  function PrimeDivisors(n)
          n < 2 && return ℤInt[]
          sort!([p for (p, e) ∈ factor(ℤ(n))])
      end
IN:   ℤTriangle(9, PrimeDivisors)

OUT:  9-element Array{Array{fmpz,1},1}:
        []
        []
        [2]
        [3]
        [2]
        [5]
        [2, 3]
        [7]
        [2]

Julia returns an ℤTriangle, even if it may not visually meet expectations. The first two entries are empty lists. This is explained on the one hand by the fact that we always evaluate a generating function starting at 0, and on the other hand because according to the fundamental theorem of arithmetic 0 and 1 do not have a prime factorization.

Such a situation occurs quite often in number theory, where many sequences start at 1 (the arithmetic functions) or even at 2. In many similar cases, the underlying reason can be traced back to the fact just given.

To handle empty rows in a ℤTriangle is no problem, after all an item of a ℤTriangle is defined to be a list, and this includes the case of an empty list. And we will uniformly treat the case 'undefined' like the case 'nonexistent' by assigning them the empty list as their extension.

A point of confusion arises when the two-dimensional ℤTriangle has to be transformed into an one-dimensional ℤSequence. This transformation is called flattening. This situation arises, for instance, when you want to submit an integer triangle to the OEIS.

Then one has to decide how to deal with empty lists. One way is to circumvent the problem and to introduce an 'offset', start the sequence at the first non-empty row and set its index to this offset. This approach is a handball trick which often leads to all sorts of inconsistencies in the long run. And it forces the bookkeeping of another quantity, the offset — a nightmare for the editors of the OEIS.

The second possibility is to explicitly define a value which takes the place of the undefined or nonexistent value (the typical phrase in definitions in the OEIS is: "... or 0 if there is none.") This may be unsatisfactory at times, but it minimizes the problems. In any case, that's the approach we choose.

In fact, this approach can also be defined in an abstract way by introducing "pointed sets" as is done in universal algebra (see Wikipedia). But that would be far beyond our goal here. For our purposes it is sufficient to define the mapping 'flattening' ℤTriangle ↦ ℤSequence as follows: 

    DEF:  function flat(T::ℤTriangle)
             Empty(s) = isempty(s) ? [ℤ(0)] : s
             [z for t ∈ T for z ∈ Empty(t)]
          end

This function is analogous to the sum function which uses the convention that the empty sum is assigned the value 0. Indeed, the problem is more general: How to collect empty collections? There is no general valid answer. One has to give an explicite definition or follow a convention.

In all cases, however, it turns out that it is more useful to set a value for the empty set than to keep it undefined. Even if the original interpretation of the sequence does not directly extend to this case. This is sometimes summed up with the phrase: "Defend the formula against interpretation." In particular, because there may be other interpretations that fit the formula better.

Before we proceed we summarize the six basic types/constructors:

    DEF:  
    const ℤSequence = Array{ℤInt,1}
    const ℤTriangle = Array{ℤSequence,1}

    ℤSequence(len) = ℤSequence(undef,len)
    ℤTriangle(dim; reg=false) = reg ?
                     ℤSequence.(1:dim) : ℤTriangle(undef,dim)

    ℤSequence(len, f::Function) = f.(0:dim-1)
    ℤTriangle(dim, f::Function) = f.(0:dim-1)

The generating function f is expected to be defined for all integers ≥ 0, in the case of a ℤSequence to return an integer of type ℤInteger, and in the case of a ℤTriangle to return a ℤSequence.

The indexed Lah triangle

The next example is a classic and important number triangle, the Lah numbers. To warm-up, we first show the implementation of which we said above that we do not want, the indexed triangle. 

DEF:  function LahIndexed(n, k)
         function recLah(n, k)
             k < 0 && return ℤ(0)
             k == n && return ℤ(1)
             recLah(n-1, k-1) + recLah(n-1, k)*(n+k-1)
         end
         recLah(n, k)
      end
IN:   [[LahIndexed(n, k) for k ∈ 0:n] for n ∈ 0:6]

OUT:  7-element Array{Array{fmpz,1},1}:
      [1]
      [0, 1]
      [0, 2, 1]
      [0, 6, 6, 1]
      [0, 24, 36, 12, 1]
      [0, 120, 240, 120, 20, 1]
      [0, 720, 1800, 1200, 300, 30, 1]

The cached Lah triangle

Before we get to our actual implementation, we have to discuss a technical matter. We use recursion. That is good because it makes the algorithm easy to understand. That is bad because it makes the execution for large terms slow. But this can be overcome if we use a cache. We use a dictionary for it. The declaration looks in Julia like this: 

const CacheLah = Dict{Int, ℤSequence}([0 => [ℤ(1)]])

We do not want to explain all the syntactical details. Julia is well documented. But one should understand that here an integer key is associated with a ℤSequence, namely the row index n to the n-th row of the triangle, and that indexing starts at 0. We also put the base case of recursion into the cache: [0 => [ℤ(1)]] as the top of the triangle is the row [ℤ(1)].

All we have to do is to weave that cache into the recursion given above. This is done in the next function below. In particular, note that the call uses only one parameter, the row index. 

DEF:  function LahNumbers(n)
          haskey(CacheLah, n) && return CacheLah[n]
          prevrow = LahNumbers(n-1)
          row = ℤSequence(n+1)
          row[1] = 0; row[n+1] = 1
          for k ∈ 2:n
              row[k] = prevrow[k-1] + prevrow[k]*(n+k-2)
          end
          CacheLah[n] = row
      end

Let's test this function:

IN:   LahNumbers(6)

OUT:  7-element Array{fmpz,1}:
      0, 720, 1800, 1200, 300, 30, 1

We derive two further functions from this basic function. The first one is: 

DEF:  LahNumbers(n, k) = LahNumbers(n)[k+1]

With this function, we can retrieve the elements of the triangle as defined in mathematics. The irritating +1 on the right side of the definition is due to the unfortunate fact that Julia starts indexing arrays at 1. (Perhaps the only really unattractive feature of this otherwise well-thought-out language).

Note how intuitive the provided notation now is. LahNumbers(n) denotes the n-th row and LahNumbers(n,k) the k-th element of the n-th row. If we were to access the elements via a Julia matrix, we would always have to write L[n+1, k+1], which would be unnatural and error-prone.

The second function uses the cache directly and returns the triangle in all its glory:

DEF:  function LahTriangle(size)
          length(CacheLah) < size && LahNumbers(size)
          [CacheLah[n] for n ∈ 0:size-1]
      end
IN:   LahTriangle(7)

OUT:  7-element Array{Array{fmpz,1},1}:
      [1]
      [0, 1]
      [0, 2, 1]
      [0, 6, 6, 1]
      [0, 24, 36, 12, 1]
      [0, 120, 240, 120, 20, 1]
      [0, 720, 1800, 1200, 300, 30, 1]

The Lah Triangle.

Let's take a quick look at the performance. 

IN:   @time LahTriangle(100);
OUT:  0.004907 seconds (39.79 k allocations: 947.781 KiB)

That's satisfactory. In a further use of the triangle, the computing time vanishes almost completely thanks to caching. A second call of the timing macro confirms this: 

IN:   @time LahTriangle(100);
OUT:  0.000009 seconds (5 allocations: 1.031 KiB)

This finishes our example. But we still want to gather the code in a small package to send it to our friends. It will essentially look like this: 

DEF:  module Lah
          export LahNumbers, LahTriangle
          const CacheLah = Dict{Int, ...
          function LahNumbers(n) ...
          LahNumbers(n, k) = LahNumbers(n+1)[k+1]
          function LahTriangle(size) ...
      end # module

The module provides (via export) two methods for the generic function LahNumbers and one method for the generic function LahTriangle. In the second part of this tutorial another important element will be added to the export list: the LahTransform.

Some characteristics of a triangle

The two most straightforward characteristics of an integer triangle are the sum of the rows and the alternating sum of the rows. 

IN:   [sum(row) for row ∈ LahTriangle(9)]

OUT:  9-element Array{fmpz,1}:
      1,  1, 3, 13, 73, 501, 4051, 37633, 394353

For this we could introduce a function rowsum like this: 

DEF:  rowsum(T) = [sum(row) for row ∈ T]

But one can also take a different view. One may ask: "What is the sum of a finite sequence?" A: "It is the sum of the terms of the sequence." Q: "What is the sum of an integer triangle?" A: "It is the (sum of the terms) of the terms of the integer triangle." Thus we might interpret the sequence of the row sums as the sum of an integer triangle. This way of thinking does not introduce a new concept, like 'rowsum', but broadens the scope of the notion 'sum.'

This alternative way of thinking, and this is perhaps surprising, fits nicely with Julia. (Unfortunately, the idea comes under a frightening name: 'multiple dispatch' or 'dynamically dispatched polymorphism'.)

So what does this look like in concrete terms? We broaden the use of sum to ℤTriangles introducing the following definition, that reflects precisely what we just described: 

DEF:  sum(T::ℤTriangle) = sum.(T)

IN:   sum(LahTriangle(9))
OUT:  9-element Array{fmpz,1}:
      1, 1, 3, 13, 73, 501, 4051, 37633, 394353

(Footnote: We must mention a technicality: we have first to hide "import Base.sum" somewhere in the header of our script before using this definition, but then we can forget about this detail.)

Before we go further, we have to clarify one more thing. What to do with an empty sum? We already discussed this above in connection with empty lists. The generally accepted convention is to return 0 in this case. This is also what Julia does if the sum runs over integers. It also works with fmpz's if you use the Nemo library with version ≥ 0.16.1. Otherwise, you have to add the line: 

IN:   Base.zero(::Type{fmpz}) = ℤ(0)

Now we can also compute the sum of the irregular triangle of the prime divisors of n in this simple, concise form:

IN:   sum(PrimeDivisors.(0:7))

OUT:  8-element Array{fmpz,1}:
      0, 0, 2, 3, 2, 5, 5, 7

Apart from the fact that our sequence starts at offset 0 (and has therefore an additional 0), it agrees with A008472.

Next, we introduce the alternating row sum. Again, we first define the alternating sum for arrays and then declare the alternating row sum of an ℤTriangle as the alternating row sums of its components (i.e., of the rows). Observe how nicely the broadcast notation reflects this relation avoiding irrelevant syntactical fluff (indices). 

DEF:  evensum(A) = sum(A[1:2:end])
      oddsum(A)  = sum(A[2:2:end])
      altsum(A)  = evensum(A) - oddsum(A)

      evensum(T::ℤTriangle) = evensum.(T)
      oddsum( T::ℤTriangle) = oddsum.(T)
      altsum( T::ℤTriangle) = evensum(T) - oddsum(T)

Let's see how this works.

IN:   T = LahTriangle(10)
      println.([sum(T), evensum(T), oddsum(T), altsum(T)]);

OUT:  fmpz[1, 1, 3, 13, 73, 501, 4051, 37633, 394353, 4596553]
      fmpz[1, 0, 1, 6, 37, 260, 2101, 19362, 201097, 2326536]
      fmpz[0, 1, 2, 7, 36, 241, 1950, 18271, 193256, 2270017]
      fmpz[1, -1, -1, -1, 1, 19, 151, 1091, 7841, 56519]

The prefix in the output fmpz[...] is a bit annoying. In the accompanying notebook we have therefore defined a function Println that suppresses these type prefixes.

The diagonalized triangle

There is still an important sum missing: the sum of the antidiagonals which we will denote by diagsum(T). It is defined as the sum of the diagonalized triangle T. So diagsum(T) = sum ⚬ diag(T) where diag is a triangle to triangle transformation mapping T to its diagonal form. We will assume that T is a regular triangle. The result, however, is an irregular triangle. This may be a psychological reason why this transformation is so rarely used. 

DEF:  function DiagonalTriangle(T::ℤTriangle)
          dim = length(T)
          U = ℤTriangle(dim)
          for n ∈ 1:dim
              R = ℤSequence(div(n+1, 2))
              for k ∈ 0:div(n-1, 2)
                  R[k+1] = T[n-k][k+1]
              end
              U[n] = R
          end
          U
      end

Applied to the Lah triangle we get DiagonalTriangle( LahTriangle) = A330609 (see also A221913). 

OUT:  [1]
      [0]
      [0, 1]
      [0, 2]
      [0, 6,    1]
      [0, 24,   6]
      [0, 120,  36,   1]
      [0, 720,  240,  12]
      [0, 5040, 1800, 120, 1]

DEF:  diagsum(T::ℤTriangle) = sum(DiagonalTriangle(T))
IN:   diagsum(LahTriangle(9)) # A001053 Sloane & Guy 1973

OUT:  9-element Array{fmpz,1}:
      1, 0, 1, 2, 7, 30, 157, 972, 6961

The comment "A001053 Sloane & Guy 1973" means that this was the earliest date the sequence appeared in the EIS. There is no indication that the connection with the Lah numbers was known to the authors.

The profile of an integer triangle

A study of an integer triangle should start with a look at the profile of the triangle, i.e., a look at the most simple characteristics of the triangle. In addition to the collection of sum functions, we add three others: the left and the right side and the central axis of the triangle. This set of characteristics is distinguished by the fact that it depends only on the row index, but not on the column index. 

DEF:  central(A) = A[div(end+1, 2)]
      central(T::ℤTriangle) = central.(T[1:2:end])

IN:   central(LahTriangle(16)) |> println

OUT:  [1, 2, 36, 1200, 58800, 3810240, 307359360, 29682132480]
DEF:  leftside(A)  = A[begin]
      rightside(A) = A[end]

DEF:  leftside( T::ℤTriangle) = leftside.(T)
      rightside(T::ℤTriangle) = rightside.(T)

The left and right sides of an integer triangle are introduced here as a generalization of the start and endpoint of an integer sequence. In the spirit of our multidispatch philosophy, a common name is used for both cases. The generalization lifts the meaning of the components on the collection via broadcasting.

This is in clear contrast to the usual names in OEIS: there the left side is called "first column", the right side is called "main diagonal" and the beginning of a sequence is called offset. Nobody seems to mind that a triangle has no main diagonal at all.

Let's put things together. A call to the profile function below provides an overview of sequences derived from the triangle, closely connected to the structure of the triangle. 

DEF:  function TraitCard(T::ℤTriangle)
          for row ∈ T println(row) end; println()
          print("Sum:      "); sum(T)       |> println
          print("EvenSum:  "); evensum(T)   |> println
          print("OddSum:   "); oddsum(T)    |> println
          print("AltSum:   "); altsum(T)    |> println
          print("DiagSum:  "); diagsum(T)   |> println
          print("Central:  "); central(T)   |> println
          print("LeftSide: "); leftside(T)  |> println
          print("RightSide:"); rightside(T) |> println
       end

This function allows for creating surveys in no time, which one wishes the OEIS would provide by higher-level search functions. Much time would be saved in the study of triangles.

🔶 Lah Δ — Trait Card 

Triangle:  1, 0, 1, 0, 2, 1, 0, 6, ...  A271703 Luschny    2016
Sum:       1, 1, 3, 13, 73, 501, ...    A000262 Sloane     1973
EvenSum:   1, 0, 1, 6, 37, 260, ...     A088312 Jovovic    2003
OddSum:    0, 1, 2, 7, 36, 241, ...     A088313 Jovovic    2003
AltSum:    1, -1, -1, -1, 1, 19, ...    A293116 Manyama    2017
DiagSum:   1, 0, 1, 2, 7, 30, 157, ...  A001053 Sloane     1973
Central:   1, 2, 36, 1200, 58800, ...   A187535 Munarini   2011
LeftSide:  1, 0, 0, 0, 0, 0, 0, 0, ...  A000007 Sloane     1994
RightSide: 1, 1, 1, 1, 1, 1, 1, 1, ...  A000012 Sloane     1994
PosHalf:   1, 1, 5, 37, 361, 4361, ...  A025168 Meeussen   1999
NegHalf:   1, 1, -3, 13, -71, 441, ...  A318223 Gutkovskiy 2018

The two traits PosHalf and NegHalf will be introduced in the second part of this article, along with two others.

The date gives the year the sequence was entered into the EIS, where EIS stands for Handbook or Encyclopedia or OEIS. We added author and date because this illustrates the mosaic-like character of the OEIS.

Each of the sequences listed here is rightly part of the OEIS and is appealing for multiple reasons. Sloane even wrote an article about one of these sequences. But it took 45 years to collect them all. If the systematic approach proposed here had already been used in 1973, it would not have taken more than half an hour to obtain all the sequences, even without the help of Julia.

Any study of a number triangle should begin with an overview like the one given in the table. We call such a table an "Integer Triangle Trait Card." In the following sections, we will create another dozen such cards. To be able to create these cards by calling a single function is the fruit of our programming efforts.

Relative offsets

We want to emphasize one thing: at no point, we used an offset when creating the above table. We defined the integer triangle we started from to be this:

       [[1], [0,1], [0,2,1], ...],
and only this; there is no 'offset' attached.

Well, it is true that we use an offset internally, for example when setting up the cache of a triangle. But we see indices as an implementation detail a user should normally not have to deal with.

We know that all derived sequences (the traits of the triangle) have, by construction, the same start as the triangle. It is this relative offset, the offset of the sequences in relation to the triangle, which is essential and which allows the interpretation of one sequence to be transferred directly to the others. This transition from a fixed offset to a relative one feels similar to the transition from Euclidean/Cartesian geometry to affine geometry.

We think it would be a good idea to give some core triangles a normative power and align the OEIS offsets of the derived sequences to them. Such a systematic is, of course, far from what we see in our beloved haystack. Such consistency is currently not even given with the two Stirling triangles. One should not underestimate how irritating this can be for some users of the database.

Three more traits

We add three more traits to the above list. They are triangles: the inverse Lah triangle, the triangle of the antidiagonals we already introduced, and the polynomial values of the Lah triangle, which, however, is best seen as an array. This array, the PolyLah array, will be defined in the next part of the article when we discuss the interpretation of integer triangles as sequences of polynomials. 

InvLah:  1, 0, 1, 0, -2, 1, 0, 6, -6,  ... A111596 Lang    2005
DiagLah: 1, 0, 0, 1, 0, 2, 0, 6, 1, 0, ... A330609 Luschny 2019
PolyLah: 1, 0, 1, 0, 1, 1, 0, 3, 2, 1, ... A253286 Luschny 2015

In the second part of this blog post, we will introduce the missing traits, discuss the connection with polynomials before we dive deeper in part three into the mathematical side of the story and look at several methods to generate interesting integer triangles, the Schröder triangles, the Motzkin triangle, the normal and the extended Catalan triangle, the Jacobsthal triangle and the Delannoy triangle.

Get started!

If a reader can't wait to try out Julia Language now, here are the next steps:

Part 2 of the tutorial       Part 3 of the tutorial

↻