Fun With The Pascal Triangle 3

Posted by Cleve Moler,

The Wikipedia article on Pascal's Triangle has hundreds of properties of the triangle and there are dozens of other Web pages devoted to it. Here are a few facts that I find most interesting.

Contents

Blaise Pascal

Blaise Pascal (1623-1662) was a 17th century French mathematician, physicist, inventor and theologian. His Traité du triangle arithmétique (Treatise on Arithmetical Triangle) was published posthumously in 1665. But this was not the first publication about the triangle. Various versions appear in Indian, Chinese, Persian, Italian and other manuscripts centuries before Pascal.

Binomial Coefficients

The binomial coefficient usually denoted by ${n} \choose {k}$ is the number of ways of picking $k$ unordered outcomes from $n$ possibilities. These coefficients appear in the expansion of the binomial $(x+1)^n$. For example, when $n = 7$
  syms x
  n = 7;
  x7 = expand((x+1)^n)
x7 =
x^7 + 7*x^6 + 21*x^5 + 35*x^4 + 35*x^3 + 21*x^2 + 7*x + 1
Formally, the binomial coefficients are given by $${{n} \choose {k}} = \frac {n!} {k! (n-k)!}$$ But premature floating point overflow of the factorials makes this an unsatisfactory basis for computation. A better way employs the recursion $$ {{n} \choose {k}} = {{n-1} \choose {k}} + {{n-1} \choose {k-1}}$$ This is used by the MATLAB function nchoosek(n,k).

Pascal Matrices

MATLAB offers two Pascal matrices. One is symmetric, positive definite and has the binomial coefficients on the antidiagonals.
   P = pascal(7)
P =
     1     1     1     1     1     1     1
     1     2     3     4     5     6     7
     1     3     6    10    15    21    28
     1     4    10    20    35    56    84
     1     5    15    35    70   126   210
     1     6    21    56   126   252   462
     1     7    28    84   210   462   924
The other is lower triangular, with the binomial coefficients in the rows. (We will see why the even numbered columns have minus signs in a moment.)
   L = pascal(7,1)
L =
     1     0     0     0     0     0     0
     1    -1     0     0     0     0     0
     1    -2     1     0     0     0     0
     1    -3     3    -1     0     0     0
     1    -4     6    -4     1     0     0
     1    -5    10   -10     5    -1     0
     1    -6    15   -20    15    -6     1
The individual elements are
P(i,j) = P(j,i) = nchoosek(i+j-2,j-1)
And (temporarily ignoring the minus signs) for i $\ge$ j
L(i,j) = nchoosek(i-1,j-1)
The first fun fact is that L is the (lower) Cholesky factor of P.
   L = chol(P)'
L =
     1     0     0     0     0     0     0
     1     1     0     0     0     0     0
     1     2     1     0     0     0     0
     1     3     3     1     0     0     0
     1     4     6     4     1     0     0
     1     5    10    10     5     1     0
     1     6    15    20    15     6     1
So we can reconstruct P from L.
   P = L*L'
P =
     1     1     1     1     1     1     1
     1     2     3     4     5     6     7
     1     3     6    10    15    21    28
     1     4    10    20    35    56    84
     1     5    15    35    70   126   210
     1     6    21    56   126   252   462
     1     7    28    84   210   462   924

Pascal Triangle

The traditional Pascal triangle is obtained by rotating P clockwise 45 degrees, or by sliding the rows of L to the right in half increments. Each element of the resulting triangle is the sum of the two above it.
   triprint(L)
                         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

Square Root of Identity

When the even numbered columns of L are given minus signs the matrix becomes a square root of the identity.
   L = pascal(n,1)
   L_squared = L^2
L =
     1     0     0     0     0     0     0
     1    -1     0     0     0     0     0
     1    -2     1     0     0     0     0
     1    -3     3    -1     0     0     0
     1    -4     6    -4     1     0     0
     1    -5    10   -10     5    -1     0
     1    -6    15   -20    15    -6     1
L_squared =
     1     0     0     0     0     0     0
     0     1     0     0     0     0     0
     0     0     1     0     0     0     0
     0     0     0     1     0     0     0
     0     0     0     0     1     0     0
     0     0     0     0     0     1     0
     0     0     0     0     0     0     1
Here is an exercise for you. What is sqrt(eye(n))? Why isn't it L?

Cube Root of Identity

When I first saw this, I was amazed. Rotate L counterclockwise. The result is a cube root of the identity.
   X = rot90(L,-1)
   X_cubed = X^3
X =
     1     1     1     1     1     1     1
    -6    -5    -4    -3    -2    -1     0
    15    10     6     3     1     0     0
   -20   -10    -4    -1     0     0     0
    15     5     1     0     0     0     0
    -6    -1     0     0     0     0     0
     1     0     0     0     0     0     0
X_cubed =
     1     0     0     0     0     0     0
     0     1     0     0     0     0     0
     0     0     1     0     0     0     0
     0     0     0     1     0     0     0
     0     0     0     0     1     0     0
     0     0     0     0     0     1     0
     0     0     0     0     0     0     1

Sierpinski

Which binomial coefficients are odd? It's a fledgling fractal.
   odd = @(x) mod(x,2)==1;
   n = 56;
   L = abs(pascal(n,1));
   spy(odd(L))
   title('odd(L)')

Fibonacci

The sums of the antidiagonals of L are the Fibonacci numbers.
   n = 12;
   A = fliplr(abs(pascal(n,1)))
   for k = 1:n
       F(k) = sum(diag(A,n-k));
   end
   F
A =
     0     0     0     0     0     0     0     0     0     0     0     1
     0     0     0     0     0     0     0     0     0     0     1     1
     0     0     0     0     0     0     0     0     0     1     2     1
     0     0     0     0     0     0     0     0     1     3     3     1
     0     0     0     0     0     0     0     1     4     6     4     1
     0     0     0     0     0     0     1     5    10    10     5     1
     0     0     0     0     0     1     6    15    20    15     6     1
     0     0     0     0     1     7    21    35    35    21     7     1
     0     0     0     1     8    28    56    70    56    28     8     1
     0     0     1     9    36    84   126   126    84    36     9     1
     0     1    10    45   120   210   252   210   120    45    10     1
     1    11    55   165   330   462   462   330   165    55    11     1
F =
     1     1     2     3     5     8    13    21    34    55    89   144

pi

The elements in the third column of lower triangular Pascal matrix are the triangle numbers. The n-th triangle number is the number of bowling pins in the n-th row of an array of bowling pins. $$t_n = {{n+1} \choose {2}}$$
   L = pascal(12,1);
   t = L(3:end,3)'
t =
     1     3     6    10    15    21    28    36    45    55
Here's an unusual series relating the triangle numbers to $\pi$. The signs go + + - - + + - - .
pi - 2 = 1 + 1/3 - 1/6 - 1/10 + 1/15 + 1/21 - 1/28 - 1/36 + 1/45 + 1/55 - ...
    type pi_pascal
function pie = pi_pascal(n)
tk = 1;
s = 1;
for k = 2:n
    tk = tk + k;
    if mod(k+1,4) > 1
        s = s + 1/tk;
    else
        s = s - 1/tk;
    end
end
pie = 2 + s;
Ten million terms gives $\pi$ to 14 decimal places.
    format long
    pie = pi_pascal(10000000)
    err = pi - pie
pie =
   3.141592653589817
err =
    -2.398081733190338e-14

Matrix Exponential

Finally, I love this one. The solution to the (potentially infinite) set of ordinary differential equations $\dot{x_1} = x_1$ $\dot{x_j} = x_j + (j-1) x_{j-1}$ is $x_j = e^t (t + 1)^{j-1}$ This means that the matrix exponential of the simple diagonal matrix
  D = diag(1:7,-1)
D =
     0     0     0     0     0     0     0     0
     1     0     0     0     0     0     0     0
     0     2     0     0     0     0     0     0
     0     0     3     0     0     0     0     0
     0     0     0     4     0     0     0     0
     0     0     0     0     5     0     0     0
     0     0     0     0     0     6     0     0
     0     0     0     0     0     0     7     0
is
  expm_D = round(expm(D))
expm_D =
     1     0     0     0     0     0     0     0
     1     1     0     0     0     0     0     0
     1     2     1     0     0     0     0     0
     1     3     3     1     0     0     0     0
     1     4     6     4     1     0     0     0
     1     5    10    10     5     1     0     0
     1     6    15    20    15     6     1     0
     1     7    21    35    35    21     7     1

Thanks

Thanks to Nick Higham for pascal.m, gallery.m and section 28.4 of N. J. Higham, Accuracy and Stability of Numerical Algorithms, Second edition, SIAM, 2002. http://epubs.siam.org/doi/book/10.1137/1.9780898718027

Get the MATLAB code

Published with MATLAB® R2018a

Note

Comments are closed.

3 CommentsOldest to Newest

Niall Madden replied on : 1 of 3
Some lovely examples here: I look forward to using them in class. However, I'm a little surprised by the assertion that MATLAB's nchoosek uses Pascal's identity. I have R2014a. For that version, it appears to use factorials (though not explicitly computed). For example, when the data are of type int64, it calls binCoef: which is described as % For integers, compute N!/((N-K)! K!) using prime factor cancellations It also has the curious behaviour of not allowing the first argument to be zero, even though this is well-defined, mathematically. Still: a fun post. Thanks for sharing. Niall.
José-Javier Martínez replied on : 3 of 3
Dear Professor Moler: After reading your beautiful post on Pascal matrices, I would like to add a new interesting aspect of them, already commented by Professor Higham in Section 28.4 of his book: “A more subtle property of the Pascal matrix is that it is totally positive”. In connection wit MATLAB, an “invisible hand” mysteriously invented the function ones(n,n). It happens that the matrix B = ones(n,n) is precisely the bidiagonal decomposition of the Pascal matrix of order n, and the work of Demmel and Koev has shown us the usefulness of that parameterization of totally nonnegative matrices. So, if we visit the software package of Plamen Koev “Accurate Computations with Totally Nonnegative Matrices” (http://www.math.sjsu.edu/~koev/software/TNTool.html) we learn that the command P=TNExpand(B) computes the Pascal matrix of order n. More importantly, without computing the Pascal matrix the command TNEigenvalues(B) computes (with high relative accuracy) the corresponding eigenvalues, even for not small values of n. Several additional functions are included in the package. Thank you, always, for your work. José-Javier Martínez