# Create a Pascal Matrix

I am trying to write a function that can create a Pascal matrix. So far I have

P[m_, n_] := Table[(#1!/(#2! (#1 - #2)!)) & /@ Thread[i + j - 1, i - 1],
{i, 1, m}, {j, 1, n}] // MatrixForm


but this doesn't give me the right matrix and I am not sure what is wrong.

Can someone tell me where the error is? Thanks.

A much simpler solution:

P[m_, n_] := MatrixForm@Table[Binomial[i + j-2, i-1], {i, 1, m}, {j, 1, n}];


or something closer to your attempt:

P[m_, n_] :=  Table[Apply[(#1!/(#2! (#1 - #2)!)) &, {i + j - 2, i - 1}],
{i, 1, m}, {j, 1, n}] // MatrixForm

• Can you tell what is wrong with my code? Nov 24, 2015 at 3:59
• @Simple You used Thread wrong. You might have confused Threadwith Apply? You also need to use i+j-2. I'll add a solution that's closer to your style.
– paw
Nov 24, 2015 at 4:05

Oh poor Stephen Wolfram.. no one honors his life's work by giving an answer that uses a cellular automaton. Let me fix this. Here a version for square matrices

P[n_] := #.Transpose[#] &@
CellularAutomaton[{{i_, j_} :> i + j}, {{1}, 0}, n - 1] And for rectangular arrays, just truncate the above solution

P[n_, m_] := P[Max[n, m]][[;; n, ;; m]]

• I have a solution using Feynman diagrams too, but seems off topic :) Nov 24, 2015 at 5:16

From here

l[n_] := SparseArray[{i_, j_} -> Binomial[i - 1, j - 1], n]
u[n_] := SparseArray[{i_, j_} -> Binomial[j - 1, i - 1], n]
s[n_] := l[n].u[n]

l // MatrixForm
u // MatrixForm
s // MatrixForm Or equivalently:

l[n_] := MatrixExp@SparseArray[Band[{2, 1}] -> Range[n - 1], n]
u[n_] := MatrixExp@SparseArray[Band[{1, -n + 1}] -> Range[n - 1], n]
s[n_] := l[n].u[n]

l // MatrixForm
u // MatrixForm
s // MatrixForm

• Alternatively: Array[1/((#1 + #2 - 1) Beta[#1, #2]) &, {7, 7}]. Then the triangles can be extracted via CholeskyDecomposition[]. Nov 24, 2015 at 4:22

Here is a nicely compact solution, whose proof is left to the interested reader:

pascal[n_Integer?Positive] := NestList[Accumulate, ConstantArray[1, n], n - 1]


It's surprisingly quick:

With[{n = 50},
AbsoluteTiming[Array[Binomial[#1 + #2, #1] &, {n, n}, {0, 0}];]]
{0.054873, Null}

AbsoluteTiming[pascal;]
{0.001829, Null}


You can use LinearAlgebraPascalMatrix.

Example:

In:= LinearAlgebraPascalMatrix

Out= {{1, 1, 1}, {1, 2, 3}, {1, 3, 6}}


It is very quick too:

In:= AbsoluteTiming[LinearAlgebraPascalMatrix;]

Out= {0.000190, Null}


It has one option:

In:= Options[LinearAlgebraPascalMatrix]

Out= {WorkingPrecision -> ∞}