Integer Triangle Traits
A Julia Implementation
Part II
Outline
- Integer triangle trait cards
- Binomial transform
- Implementing the binomial transform
- Transformed sequences as traits of triangles.
- The traits N0TS and NATS.
- The inverse binomial transform
- The Fubini/Worpitzky Triangles
- A prototypical triangle module
- Julia's Romeo is Nemo
- Triangles as skeletons of polynomials
- Generating polynomials
- Evaluating polynomials.
- Two examples: evaluations at ± 1/2.
- The traits poshalf and neghalf.
- Bivariate exp. generating functions
- The Lah and the Bernoulli polynomials
- Reinterpretation as a rectangle
- Flattening a sequence of triangles
- The polynomial array
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.
- In the first part, we put together some basic terminology and concepts.
- In the second part we examine the connection with polynomial sequences.
- The third part deals with general algorithms for constructing 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.
Integer triangle trait cards
Let us begin the second part by re-focusing on the goal of our efforts: the creation of "Integer Triangle Trait Cards"™. This is a compilation of the essential characteristics of a triangle, whereby we understand the characteristics of a triangle to be sequences that can be obtained from the triangle by elementary transformations.
The definition of most traits were given in the last chapter. The missing traits, PosHalf, NegHalf, N0TS and NATS will be explained below.
A first example is the famous triangle of binomial coefficients.
🔶 Pascal's Δ
[0] 1 [1] 1 1 [2] 1 2 1 [3] 1 3 3 1 [4] 1 4 6 4 1 [5] 1 5 10 10 5 1 [6] 1 6 15 20 15 6 1 [7] 1 7 21 35 35 21 7 1 [8] 1 8 28 56 70 56 28 8 1 [9] 1 9 36 84 126 126 84 36 9 1 Trait Seq Anum Author Year Δ Binomial: 1, 1, 1, 1, 2, 1, 1, 3, ... A007318 Sloane 1994 Δ InvBinom: 1, -1, 1, 1, -2, 1, -1, ... A130595 Deléham 2007 Δ DiagBinom: 1, 1, 1, 1, 1, 2, 1, 3, ... A011973 Sloane 1996 Δ PolyBinom: 1, 1, 1, 1, 2, 1, 1, 4 ... A009998 Sloane 1996 -- Sum: 1, 2, 4, 8, 16, 32, 64, ... A000079 Sloane 1973 EvenSum: 1, 1, 2, 4, 8, 16, 32, ... A011782 Killough 1996 OddSum: 0, 1, 2, 4, 8, 16, 32, ... A131577 Curtz 2007 AltSum: 1, 0, 0, 0, 0, 0, 0, ... A000007 Sloane 1994 DiagSum: 1, 1, 2, 3, 5, 8, 13, ... A000045* Sloane 1964 Central: 1, 2, 6, 20, 70, 252, ... A000984 Sloane 1973 LeftSide: 1, 1, 1, 1, 1, 1, 1, ... A000012 Sloane 1994 RightSide: 1, 1, 1, 1, 1, 1, 1, ... A000012 Sloane 1994 PosHalf: 1, 3, 9, 27, 81, 243, ... A000244 Sloane 1973 NegHalf: 1, -1, 1, -1, 1, -1, ... A033999 Danilov 1998 N0TS: 0, 1, 4, 12, 32, 80, ... A001787 Sloane 1973 NATS: 1, 3, 8, 20, 48, 112, ... A001792 Sloane 1973 -- A000045*(n) = A000045(n+1)
The date gives the year the sequence was entered into the EIS,
where EIS stands for Handbook or Encyclopedia or OEIS.
The sequence A000045*(n) = A000045(n+1)
does not have an independent entry. Judging from the entries in A000045
and by the fact, that DiagSum(Pascal) = Fibonacci is a
perfect definition for this sequence, this arguably would have been justified.
Binomial transform
Given a sequence s and a function T defined for integers 0 ≤ k ≤ n with T(n,n) != 0 ∀ n, and 0 otherwise. Then we call the sequence Ts the T-transform of s if Ts(n) = Sum(0 ≤ k ≤ n) T(n,k)*s(k) (∀ n ≥ 0).
iT is the inverse transform of T if the matrix iT(n,k) is the inverse matrix of the matrix T(n,k).
If T(n,k) = binomial(n,k) then T is the binomial transform and (n,k) ↦ iT(n,k) = (-1)^(n,k)*binomial(n,k) the inverse binomial transform. By the binomial theorem, binomial(n,k) and (-1)^(n,k)*binomial(n,k) are the connection coefficients between the bases {x^n} and {(x-1)^n}.
The definitions given here are consistent with their use in the OEIS, especially with the definition in Sloane's 'Transforms'. Unfortunately, there are Internet sources that deviate from this definition. Note that Wikipedia is unreliable in this respect and gives a wrong definition! A good reference is, for example, "Combinatorial Theory" by Martin Aigner.
Implementing the binomial transform
The conventional way to implement the binomial transform looks similar to the function below:
???: function BinomialTransform(seq)
sz = size(seq, 1)
T = ℤSequence(sz)
T[1] = seq[1]
for k ∈ 1:sz-1
next = ℤ(0)
for j ∈ 0:k
next += binomial(k, j)*seq[j+1]
end
T[k+1] = next
end
T
end
Our approach is different. To keep it general, we first extract the content of the two loops: they are nothing more than a matrix times vector multiplication, as well known from linear algebra. In fact, the transformations we are looking at are linear.
LinMap(F::Function,V::ℤSeq,n) = sum(F(n-1)[k]*V[k] for k ∈ 1:n) LinMap(F::Function,V::ℤSeq) = LinMap(F, V, length(V)) LinMap(M::ℤTri,V::ℤSeq,n) = sum(M[n][k] * V[k] for k ∈ 1:n) LinMap(M::ℤTri,V::ℤSeq) = (n -> LinMap(M, V, n)).(1:length(V))
The call parameters differ depending on whether the transformation kernel is given as a function or as a triangle; the result is the same in both cases.
[Homework] If K(n) = [ℤ(1) for _ ∈ 0:n] and L(n) = [ℤ(k) for k ∈ 0:n] what is [LinMap(K, K(n)) for n ∈ 0:9] respectively [LinMap(L, L(n)) for n ∈ 0:9]?
Now to come back to our example:
First we assume that we have two incarnations of the binomial function (called methods in Julia), Binomial(n) which returns the n-th row of Pascal's triangle, and Binomial(n,k) which returns the k-th element of the n-th row. We have already seen this interface in part 1 for the Lah numbers.
Our implementation are three one-liners and has an entirely different look and feel than the above implementation.
DEF: Telescope(V::ℤSequence) = (V[1:k] for k ∈ 1:size(V,1))
Binomial(V::ℤSequence) = LinMap(Binomial, V)
BinomialTransform(V::ℤSequence) = Binomial.(Telescope(V))
The Binomial Transform
So what is happening here? Let's start with the second line. Here the binomial function is extended to Vectors and defined as the image of the vector V under the linear transform associated with the binomial matrix.
In the third line, the binomial function is applied to the Telescope of the sequence V. The telescope of V is the expansion of V into a sequence of sequences (i.e., into a ℤTriangle), where the n-th row is the initial segment of V of length n+1.
At this point, one should note a small but important syntactical difference. We have defined Telescope(V) = (V[1:k] ...) and not [V[1:k] ...]. In other words, we used round brackets instead of square brackets. In consequence, Telescope(V) is not an array but an iterator, and the transformation is more efficient.
Much of the joy of coding with Julia derives from the fact that Julia makes it easy to formulate the notions on which the computations are based so easily and succinctly. This makes it very different from 'old style' coding, which is chasing indices through nested loops.
But this is by far not the end of the story. By purely formal substitution of the name in the last two lines (Binomial → LahNumbers), one gets the Lah transform, a transform Sloane added to his Transform library in 2003 at the suggestion of Vladeta Jovovic.
Def: LahNumbers(A::ℤSequence) = LinMap(LahNumbers, A, length(A))
LahTransform(A) = LahNumbers.(Telescope(A))
And this procedure works for any integer triangle Δ! All it needs are two functions with the same interface as the LahNumbers and the Binomial coefficients above: Δ(n) and Δ(n,k).
Def: Δ(A::ℤSequence) = LinMap(Δ, A, length(A))
ΔTransform(A::ℤSequence) = Δ.(Telescope(A))
In part 4 of this series we will discuss Julia's macro language and take things a step further (into the realm of metaprogramming) and show how one can generate all these functions with a macro for a given triangle.
Transformed sequences as traits of triangles.
Let's look at an example, where we also use the Motzkin transform derived from the Motzkin triangle defined in the next chapter:
IN: U = [ℤ(1) for _ ∈ 0:8] BinomialTransform(U) |> println LahTransform(U) |> println MotzkinTransform(U) |> println OUT: [1, 2, 4, 8, 16, 32, 64, 128, 256] [1, 1, 3, 13, 73, 501, 4051, 37633, 394353] [1, 2, 5, 13, 35, 96, 267, 750, 2123]
If we look at the transformed sequences, we see that they connect to things that are known to be closely related to the notions behind the respective triangles:
Binomial ↦ A000079; Lah ↦ A000262; Motzkin ↦ A005773. In other words, they show characteristic traits of the triangles.
The traits N0TS and NATS.
We are now adding two more to our collection of traits of a triangle, the sequences arising from the transformation of the fundamental sequences:
N0TS = ΔTransform([0,1,2,3, ...])NATS = ΔTransform([1,2,3,4, ...])
We call the transform of the positive natural numbers NATS and the transform of the non-negative natural numbers N0TS. Perhaps this is not the very summit of scientific terminology, but quite suitable for distinguishing these two cases in daily use.
Note that the important transform
UNOS = ΔTransform([1,1,1,1, ...])
is already in our list of traits: it generates the row sums of triangles.
DEFS: N0TS(T::Function) = T.(0:n-1)
NATS(T::Function) = T.(1:n)
In our context, we simply assume that the function T has been defined as a transformation in the sense of the last section.
The inverse binomial transform
In mathematics, inverse transformations are often more important than direct ones. And often more difficult to calculate. Fortunately, this is not the case here, we can reuse our scheme introduced above with a slight variation.
Again, the construction is based on the concept of the telescope of a sequence. So we assign a triangle to a sequence that consists of the initial segments of the sequence, only that this time the terms are signed.
DEF: SignedTelescope(A::ℤSequence) =
([(-1)^(k-j)*A[j] for j ∈ 1:k] for k ∈ 1:size(A,1))
SignedTelescope(len::Int, f::Function) =
([(-1)^(k-j)*ZZ(f(j-1)) for j ∈ 1:k] for k ∈ 1:len)
The first of the two incarnations corresponds exactly to the above definition of the telescope except for the factor (-1)^(k-j). The second variant returns the same thing, except that the first len terms of the sequence are calculated first with the function f. This form is sometimes more flexible in application.
The inverse binomial transformation is now defined as the binomial applied to the signed telescope of the sequence.
DEF:
BinomialInvTransform(A::ℤSequence) = Binomial.(SignedTelescope(A))
BinomialInvTransform(len::Int, f::Function) =
Binomial.(SignedTelescope(len, f))
The Inverse Binomial Transform
It is obvious that in this construction the binomial function has only an exemplary character. Just as described in the case of the direct transformation, we can use this setup as a blueprint for the definition of any inverse transformation which we will consider.
ΔInvTransform(A::ℤSequence) = Δ.(SignedTelescope(A)) ΔInvTransform(len::Int, f::Function) = Δ.(SignedTelescope(len,f))
This suggests to take the final step and remove the type of the transformation from the name and pass it in as parameter instead.
InvTransform(Δ, A::ℤSequence) = Δ.(SignedTelescope(A))
InvTransform(Δ, len::Int, f::Function) =
Δ.(SignedTelescope(len, f))
The Fubini/Worpitzky Triangles
As a first application we consider the calculation of the Fubini triangle A131689.
The question is: Consider the power functions (x,n) -> x^n.
What is the inverse binomial transform of this sequence?
Since the functions depend on a parameter n, let's look at the entire sequence
of the transforms and put them row by row into a triangle. One possible implementation is:
function SomeTriangle(dim::Int, f::Function)
T = ℤTriangle(dim)
for n ∈ 0:dim-1
T[n+1] = BinomialInvTransform(n+1, x -> f(x, n))
end
T
end
After a brief consideration, it is clear that such a function could also be useful for transformations other than the inverse binomial transform. Also, we do not like for-loops. In fact we can avoid the loop here if we remember the constructors of triangles. This brings us to the following form:
TransformTriangle(Trans::Function, dim::Int, f::Function) =
ℤTriangle(dim, n -> Trans(n+1, x -> f(x, n)))
A hat tip goes to Haskell Curry. Let's look at the result, the Fubini triangle, which is the inverse binomial transform of the power functions:
IN: FubiniTriangle(dim) = TransformTriangle(BinomialInvTransform, dim, (x,n) -> x^n) FubiniTriangle(8) OUT: 8-element Array{Array{fmpz,1},1}: [1] [0, 1] [0, 1, 2] [0, 1, 6, 6] [0, 1, 14, 36, 24] [0, 1, 30, 150, 240, 120] [0, 1, 62, 540, 1560, 1800, 720] [0, 1, 126, 1806, 8400, 16800, 15120, 5040]
But the story doesn't end here. The Fubini triangle has an elder brother, the Worpitzky triangle, from which it differs only in the signs. If one forms the row sums of the Worpitzky triangle in the form
Sum(k=0..n) T(n,k)/(k+1), then one gets the Bernoulli numbers,
as Worpitzky showed.
But the following is even more descriptive: If the terms of the Worpitzky triangle are interpreted
as the coefficients of polynomials Fn(x), then
Integral{x=0..1} Fn(x) = Bn(1),
where Bn(x) are the Bernoulli polynomials. In other words, the integrals equal the Bernoulli numbers. See also the illustration in A278075.
A prototypical triangle module
Before we go any further, we want to show how what has been said so far can be put together in a Julia module. As an example, we consider the Laguerre triangle.
module LaguerreTriangles
using IntegerTriangles
export Laguerre, LaguerreTriangle
export LaguerreTransform, LaguerreInvTransform
const CacheLaguerre = Dict{Int, ℤSequence}(0 => [ℤInt(1)])
function Laguerre(n::Int)
haskey(CacheLaguerre, n) && return CacheLaguerre[n]
prevrow = Laguerre(n-1)
row = ℤSequence(n+1)
row[n+1] = ℤInt(1)
for k ∈ 1:n
row[k] = ( get(prevrow, k-1, 0)
+ get(prevrow, k, 0)*(2*k-1)
+ get(prevrow, k+1, 0)*k^2
)
end
CacheLaguerre[n] = row
end
Laguerre(n::Int, k::Int) = Laguerre(n)[k+1]
Laguerre(A::ℤSequence) = LinMap(Laguerre, A)
function LaguerreTriangle(size::Int)
length(CacheLaguerre) < size && Laguerre(size)
[CacheLaguerre[n] for n ∈ 0:size-1]
end
LaguerreTransform(A::ℤSequence) = Laguerre.(Telescope(A))
LaguerreInvTransform(A::ℤSequence) =
Laguerre.(SignedTelescope(A))
end
The functions LinMap and Telescope (and of course the types ℤInt and ℤSequence) are imported from our package IntegerTriangles. We leave it up to the reader to look up other details in the Julia documentation (especially the use of the get function, which takes care that the values outside the range of prevrow are handled correctly.)
And here is the trait card of the Laguerre triangle:
🔶 Laguerre Δ
[0] 1 [1] 1 1 [2] 2 4 1 [3] 6 18 9 1 [4] 24 96 72 16 1 [5] 120 600 600 200 25 1 [6] 720 4320 5400 2400 450 36 1 [7] 5040 35280 52920 29400 7350 882 49 1 [8] 40320 322560 564480 376320 117600 18816 1568 64 1 [9] 362880 3265920 6531840 5080320 1905120 381024 42336 2592 81 1 Trait Seq Anum Author Year Δ Triangle 1, 1, 1, 2, 4, 1, 6, 18, 9, ... A021009* Sloane 1999 Δ Inverse -1, 1, 2, -4, 1, -6, 18, -9, .. A021009# Sloane 1999 Δ Diagonal 1, 1, 2, 1, 6, 4, 24, 18, 1, .. A084950 Adamson 2003 Δ Polynom 1, 1, 1, 2, 2, 1, 6, 7, 3, .. A289192 Heinz 2017 -- Sum 1, 2, 7, 34, 209, 1546, ... A002720 Sloane 1973 EvenSum 1, 1, 3, 15, 97, 745, 6571, ... A331325 Luschny 2020 OddSum 0, 1, 4, 19, 112, 801, 6756, .. A331326 Luschny 2020 AltSum 1, 0, -1, -4, -15, -56, ... A009940 Meeussen 1999 DiagSum 1, 1, 3, 10, 43, 225, 1393, ... A001040* Sloane 1973 Central 1, 4, 72, 2400, 117600, ... A295383* Gutkovskiy 2017 LeftSide 1, 1, 2, 6, 24, 120, 720, ... A000142 Sloane 1973 RightSide 1, 1, 1, 1, 1, 1, 1, 1, 1, ... A000012 Sloane 1994 PosHalf 1, 3, 17, 139, 1473, 19091, ... A025167 Meeussen 1999 NegHalf 1, -1, 1, 7, -127, 1711, ... A025166* Meeussen 1999 N0TS 0, 1, 6, 39, 292, 2505, ... A103194 Jovovic 2005 NATS 1, 3, 13, 73, 501, 4051, ... A000262* Sloane 1973 --- A021009* = unsigned A021009. A021009# = A021009 with different signs. A001040* = A001040 shifted left once. A295383* = abs(A295383). A025166* = A025166 with different signs. A000262* = A000262 shifted left once.
Julia's Romeo is Nemo
Julia is not a CAS. Julia is a computing language that was designed aiming at numerical computing and current methods of data modeling. So one cannot expect the availability of symbolic mathematics or implementations of algebraic structures like polynomial rings out of the box. But fortunately, some excellent packages exist which provide much of this. Particularly noteworthy is the package Nemo which is used in this article. It provides Julia bindings for various mathematical libraries (including flint, a fast library for number theory). William Hart and others develop it.
IN: using Nemo
In the following sections, we will make use of polynomial rings and power series imported from Nemo. We introduce some type aliases and constructors.
DEFS:
const RationalSequence = Array{fmpq,1}
const RationalTriangle = Array{Array{fmpq,1},1}
ℚSequence(len) = Vector{fmpq}(undef, len)
ℚTriangle(dim) = Vector{RationalSequence}(undef, dim)
const ℤPolySeq = Array{fmpz_poly,1}
const ℚPolySeq = Array{fmpq_poly,1}
ℤPolynomials(x) = PolynomialRing(ZZ, x)
ℚPolynomials(x) = PolynomialRing(QQ, x)
Triangles as skeletons of polynomials
There is a one-to-one correspondence between regular integer triangles that do not vanish on the right side (main diagonal) and ℤ-polynomial sequences. A polynomial sequence is a sequence of polynomials such that degree(pn(x)) = n (∀ n ≥ 0). The importance of polynomial sequences derives from the fact that they form a basis. Examples are the standard polynomials xn, the falling factorials and the rising factorials.
To illustrate this further, we have put together a vocabulary which translates the concepts introduced in Part 1 into the language of the polynomials.
| row | ⇄ | polynomial |
| rowsum | ⇄ | pn(1) |
| altsum | ⇄ | pn(-1) |
| leftside | ⇄ | pn(0) |
| rightside | ⇄ | lead(pn) |
We used this concepts to build a profile of an integer triangle. Now we add two more traits to the profile.
| poshalf | ⇄ | 2n pn(1/2) |
| neghalf | ⇄ | (-2)n pn(-1/2) |
We see that it would be difficult to express these two traits in the language of integer triangles. First, we have to learn how to handle polynomials in Julia/Nemo.
Generating polynomials
To warm-up, we show how to convert an integer triangle to a polynomial sequence and vice versa.
DEF: function Polynomial(s::ℤSequence) R, x = ℤPolynomials("x") sum(c*x^(k-1) for (k,c) ∈ enumerate(s)) end Polynomial(T::ℤTriangle) = Polynomial.(T) IN: T = [[1], [0,1], [0,2,1], [0,6,6,1], [0,24,36,12,1]] Polynomial(T) OUT: 5-element Array{fmpz_poly,1}: 1 x x^2 + 2*x x^3 + 6*x^2 + 6*x x^4 + 12*x^3 + 36*x^2 + 24*x
Unfortunately, the polynomials in the output are sorted by descending powers, whereas in the triangle, the exponent is the index in the row, read from left to right. But we can live with that.
What hurts more is x^(k-1)
in the definition of the polynomials. The -1 in the exponent is solely due to Julia's decision to have arrays begin with 1. The result is an unnatural representation of polynomials in formulas and, of course, encourages off-by-one errors ...
More pleasing is how elegantly the concept of the coefficient expands. In the first definition, it applies to a polynomial and in the second definition to a sequence of polynomials. Here again, the use of the "broadcasting operator" is the method of choice.
DEF: coefficients(p) = coeff.(p, 0:degree(p)) coefficients(P::ℤPolySeq) = coefficients.(P) IN: R, x = ℤPolynomials("x") P = [R(1),x,x^2+2*x,x^3+6*x^2+6*x,x^4+12*x^3+36*x^2+24*x] coefficients(P) OUT: 5-element Array{Array{fmpz,1},1}: [1] [0, 1] [0, 2, 1] [0, 6, 6, 1] [0, 24, 36, 12, 1]
Using this method we avoided all for-loop syntax in our definitions. Note the input of the polynomial sequence in the last example: [R(1),x,...]. Here we had to use the constructor for the constant polynomial 1 and had to write R(1); otherwise Julia would had interpreted 1 as an integer which would result in an error.
Evaluating polynomials.
Polynomials are no functions. But we can evaluate them by substitution. Here we are interested in the values at +1/2 and -1/2. This forces us to consider rational polynomials (i.e., polynomials whose coefficients are in ℚ).
DEF: function Evaluate(p, x) Q, q = PolynomialRing(QQ, "q") subst(Q(p), x) end DEF: Evaluate(P::ℤPolySeq, x) = Evaluate.(P, x) Evaluate(T::ℤTriangle, x) = Evaluate.(ℤPolySeq(T), x)
To do this, we have to convert the integer polynomial into a rational polynomial before we evaluate it. In the second definition, we apply the
evaluation to a sequence of integer polynomials, but the result is a ℚSequence. If we pay attention to type stability, we have to convert this sequence back to a ℤSequence.
The numerator function can be used for this, as we will see below.
Two examples: evaluations at ± 1/2.
The evaluation of an integer polynomial at 1/2:
IN: R, x = ℤPolynomials("x") p = x^4+12*x^3+36*x^2+24*x Evaluate(p, 1//2) OUT: 361//16
The evaluation of an integer triangle at 1/2 (via the associated polynomials):
IN: T = [[1], [0,1], [0,2,1], [0,6,6,1], [0,24,36,12,1]] Evaluate(ℤTriangle(T), 1//2) OUT: 5-element Array{fmpq,1}: 1, 1//2, 5//4, 37//8, 361//16
The triangle T as written has
the type IntTriangle because integer literals are of
type Int. Therefore T has to be converted to a
ℤTriangle before applying the evaluation function.
The traits poshalf and neghalf.
Now we can finally introduce the announced traits poshalf and neghalf of an integer triangle.
DEFS:
poshalf(p) = numerator(2^degree(p)*Evaluate(p,1//2))
neghalf(p) = numerator((-2)^degree(p)*Evaluate(p,-1//2))
poshalf(T::ℤTriangle) = poshalf.(Polynomial(T))
neghalf(T::ℤTriangle) = neghalf.(Polynomial(T))
poshalf(P::ℤPolySeq) = poshalf.(P)
neghalf(P::ℤPolySeq) = neghalf.(P)
This simply means: interpret the rows of an integer triangle as the coefficients of a polynomial, evaluate this polynomial at the x = 1/2 respectively x = -1/2 and normalize the result with 2^n respectively (-2)^n, where n
is the degree of the polynomial. Finally, and only to resurrect the type ℤInt, apply the numerator function to the resulting value.
As a test case let's look at the Laguerre triangle.
IN: T = LaguerreTriangle(8) poshalf(T) |> println neghalf(T) |> println OUT: [1, 3, 17, 139, 1473, 19091, 291793, 5129307] [1, -1, 1, 7, -127, 1711, -23231, 334391]
These are, apart from the sign-pattern, sequences A025167 and A025166 introduced by Wouter Meeussen in 1999.
Bivariate exp. generating functions
The next function returns the coefficients of a sequence of polynomials generated by a bivariate exponential generating function.
DEF: Laplace(s, k) = factorial(k) * coeff(s, k) function EgfExpansion(prec, gf::Function, coeff=true) R, x = QQPolynomials("x") S, t = PowerSeriesRing(R, prec+1, "t") ser = gf(x, t) P = Laplace.(ser, 0:prec-1) coeff ? coefficients.(P) : P end EgfExpansionCoeff(prec, gf) = EgfExpansion(prec, gf, true) EgfExpansionPoly(prec, gf) = EgfExpansion(prec, gf, false)
Note how much more precise the specifications are compared to the symbolic approach often taken with CASs like Maple or Mathematica. Nemo is more similar in its architecture to programs like Magma or SageMath. We will not go into the details and refer to the documents of the package Nemo.jl, easily to find on GitHub.
Lah and Bernoulli polynomials
We will use the above functions to demonstrate the construction of the Lah triangle from its generating function exp(t*(x/(1-t)).
DEF: G271703(x,t) = exp(t*divexact(x, 1 - t))
T271703(dim) = EgfExpansionCoeff(dim, G271703)
P271703(dim) = EgfExpansionPoly(dim, G271703)
The naming scheme used is obvious: G271703 denotes the generating function of A271703, T271703 is the associated triangle and P271703 the associated sequence of polynomials.
IN: T271703(5) OUT: 5-element Array{Array{fmpq,1},1}: [1] [0, 1] [0, 2, 1] [0, 6, 6, 1] [0, 24, 36, 12, 1] IN: P271703(5) OUT: 5-element Array{fmpq_poly,1}: 1 x x^2+2*x x^3+6*x^2+6*x x^4+12*x^3+36*x^2+24*x
What is easy to miss is the small 'q' in the output: the returned array is a rational triangle, not an integer triangle! Therefore we introduced the type alias ℚTriangle for Array{Array{ fmpq,1},1}. But we deliberately did not include the numerator function in EgfExpansion to reduce the risk that the denominators will be cut off accidentally from fractions whose denominator is not equal to 1.
Instead, we lift the numerator function from the rational numbers two levels up to rational triangles.
DEF: import Nemo.numerator
numerator(T::RationalTriangle) = [numerator.(t) for t ∈ T]
Now we can call numerator(T271703(7)) to get the Lah numbers as an integer triangle.
For a second example consider the computation of the Bernoulli polynomials.
IN: egfBernoulli(x, t) = divexact(t*exp(x*t), exp(t)-1) BernoulliPolynomials(n) = EgfExpansionPoly(n, egfBernoulli) BernoulliPolynomials(6) OUT: 6-element Array{fmpq_poly,1}: 1 x - 1//2 x^2 - x + 1//6 x^3 - 3//2*x^2 + 1//2*x x^4 - 2*x^3 + x^2 - 1//30 x^5 - 5//2*x^4 + 5//3*x^3 - 1//6*x
The Bernoulli polynomials.
In Julia rationals are constructed using the // operator. This operator has precedence over multiplication and addition. No need for brackets, no ambiguity here.
As a third example we look at the Worpitzky polynomials, which we have already considered above. Their coefficients can also be calculated from their exponential generating function like this:
IN: G278075(x, t) = divexact(1, 1 - x*(1 - exp(-t))) WorpitzkyTriangle(n) = EgfExpansionCoeff(n, G278075) WorpitzkyTriangle(7) OUT: [1] [0, 1] [0, -1, 2] [0, 1, -6, 6] [0, -1, 14, -36, 24] [0, 1, -30, 150, -240, 120] [0, -1, 62, -540, 1560, -1800, 720]
The coefficients of the Worpitzky polynomials.
Reinterpretation as a rectangle
Before we go any further, we need some more infrastructure.
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||
Informally we can say that pushing the k-th column of a triangle
up k cells transforms the triangle into the rectangle. Formally the indices are transformed (n,k) → (n+k,k) when proceeding from the triangle form to the rectangular form.
We call this a reinterpretation and not a transformation because the data does not change in this process. Reading the triangle by rows is the same as reading the rectangle by ascending antidiagonals. However, we get a different view of the data which might lead to further insights into the structure of the generating function.
Note that the first column of the triangle is identical to the first column of the rectangle. That is the supporting pillar in this reinterpretation. Make sure that this pillar remains in place also when the reading direction is reversed (rectangle ↦ triangle). In particular, this means: avoid reading the rectangle by descending antidiagonals.
This is in line with our definition of the left side of a sequence or triangle and it conforms to the western way to read text from left to right.
Flattening a sequence of triangles
How to lay out a sequence of triangles in linear storage? Let's assume a sequence of triangles, A, B, C, D, ... given as flattened triangles in a rectangular array. Thus T(n) here means the n-th row of the triangle T, in the same way as we used it in the implementations above. Now using the same reinterpretation as in the last section, this time reading from right to left, we arrive at a triangle of rows.
|
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Reading the right triangle by rows we finally have the linear layout of the sequence of triangles we looked for. From the mathematical point of view we defined an indexing of ℕ x ℕ where ℕ = {0, 1, 2, ...}. See A332662 for further remarks.
The polynomial array
The trait cards summarize 12 sequences and 4 triangles. Their names in the captions of the cards are: Δ Triangle, Δ Inverse, Δ Diagonal, Δ Polynom, Sum, EvenSum, OddSum, AltSum, DiagSum, Central, LeftSide, RightSide, PosHalf, NegHalf, N0TS, NATS.
We have now discussed all but one triangle: the polynomial triangle, which, however, is better presented as a rectangular array. It is defined as the evaluation of the polynomial associated with the n-th row on the nonnegative natural numbers. This means that the n-th row of the array is generated as indicated in the scheme below.
T(n,0), T(n,1), ..., T(n,n)
↦ p(x) = Sum_{k=0..n} T(n,k)*x^k
↦ p(0), p(1), p(2), ....
For instance, the first three lines of the Laguerre polynomial array are computed as follows:
[1] ↦ 1 ↦ 1, 1, 1, 1, ...
[1, 1] ↦ 1 + x ↦ 1, 2, 3, 4, ...
[2, 4, 1] ↦ 2 + 4*x + x^2 ↦ 2, 7, 14, 23, ...
DEFS:
function PolyArray(T::ℤTriangle)
P = Polynomial(T)
dim = length(T)
U = ℤTriangle(dim)
for n = 1:dim
p = P[n]
eva = [Evaluate(p, k) for k ∈ 0:dim-1]
U[n] = numerator.(eva)
end
U
end
function PolyTriangle(T::ℤTriangle)
A = PolyArray(T)
U = ℤTriangle(length(T))
for n = 1:length(T)
U[n] = [A[n-k][k+1] for k ∈ 0:n-1]
end
U
end
Applied to the Laguerre triangle these functions return the following triangular, respectively rectangular, arrays:
IN: T = LaguerreTriangle(5) Println.(PolyTriangle(T)) Println.(PolyArray(T)) OUT: [1] [1, 1] [2, 2, 1] [6, 7, 3, 1] [24, 34, 14, 4, 1] [1, 1, 1, 1, 1] [1, 2, 3, 4, 5] [2, 7, 14, 23, 34] [6, 34, 86, 168, 286] [24, 209, 648, 1473, 2840]
If the coefficients list of the polynomials is reversed, one gets an entirely different array whose interpretation is not clear. A setup where the lead coefficient is 1 is always preferable.
1, x+1, x^2+4*x+2, x^3+9*x^2+18*x+6, ...
1, x+1, 2*x^2+4*x+1, 6*x^3+18*x^2+9*x+1, ... (reversed)
# The reversed Laguerre array.
[1, 1, 1, 1, 1]
[1, 2, 3, 4, 5]
[1, 7, 17, 31, 49]
[1, 34, 139, 352, 709]
[1, 209, 1473, 5233, 13505]
Finally the question: How efficient is our general setup in the present case?
@time PolyArray(LaguerreTriangle(100)) 1.256 seconds (1.88 M allocations: 57.7 MiB, 8.2% gc time)
In our opinion, this result fulfills the promise: "Julia was designed to help researchers write high-level code in an intuitive syntax and produce code with the speed of production programming languages." Let's compare it to the method proposed in A289192 using Maple.
A := (n, k) -> n!*add(binomial(n, i)*k^i/i!, i=0..n): t := time(): seq(seq(A(n,k), k=0..99),n=0..99): time()-t; 4.141 seconds
However, I am not sure whether this is a fair comparison. After all, "Maple Academic Single-User Edition" currently costs US$ 995, Julia is free and open. In addition, Maple is promoted to introduce the concept of "clickable math". Julia, in contrast, does not have much to offer about clicks.
Get started!
In the third part of this blog post we dive deeper 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.
If a reader can't wait to try out Julia Language now, here are the next steps:
- Installation of Julia. Shown in this video by Carsten Bauer for Windows, and similar on other platforms.
- Download the notebook IntegerTriangles.ipynb from GitHub as a starting point.
- Recommended introduction: Ben Lauwens & Allen Downey Think Julia: How to Think Like a Computer Scientist, a free Creative Commons 3.0 online version on GitHub.
- Indispensable: Julia documentation.
- If you have questions, go to the very helpful Julia community on discourse.