Fun With The Pascal Triangle

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

Published with MATLAB® R2018a
|
  • print

评论

要发表评论,请点击 此处 登录到您的 MathWorks 帐户或创建一个新帐户。