/*
Functions for generating and manipulating Pascal matrices
 in Maxima (descendant of DOE Macsyma).

Copyright (C) 2005 Exstrom Laboratories LLC
This work is freely redistributable under the terms of the
 GNU General Public License as published by the Free Software
 Foundation. This work comes with ABSOLUTELY NO WARRANTY.
*/

/*
This function generates the Pascal matrix of the specified
 size, and raised to the specified power: P(n)^t.
Usage: pascal(n,t)
  n = size of the (square) Pascal matrix, n > 0.
  t = power of the Pascal matrix = any real number.

Note: The Pascal matrix is a lower triangular matrix.
 To make an upper triangular Pascal matrix, just transpose it.
 To make a symmetric Pascal matrix, multiply the Pascal matrix
 by its transpose on the right: P.transpose(P), or just use
 the function pascalsym(n).

Date: Mar 20, 2005
*/
pascal(n,t):=
block([P,i,j],
P:zeromatrix(n,n),
for i:0 step 1 while i<n do(
  for j:0 step 1 while j<=i do(
    P[i+1,j+1]:t^(i-j)*binomial(i,j)
  )
),
return(P)
)$

/*
This function generates the symmetric Pascal matrix
 of the specified size. It is equivalent to doing
 the following with the regular Pascal matrix:
 P.transpose(P).

Usage: pascalsym(n)
  n = size of the symmetric Pascal matrix, n > 0.

Date: Mar 24, 2005
*/
pascalsym(n):=
block([P,i,j],
P:zeromatrix(n,n),
for i:0 step 1 while i<n do(
  for j:0 step 1 while j<=i do(
    P[i+1,j+1]:binomial(i+j,j),
    P[j+1,i+1]:P[i+1,j+1]
  )
),
return(P)
)$

/*
This function will generate a shifted symmetric
 Pascal matrix of the specified size.

Usage: pascalsymshift(n,s)
  n = size of the symmetric Pascal matrix, n > 0.
  s = an integer, which is the shifting parameter.

Date: Mar 24, 2005
*/
pascalsymshift(n,s):=
block([P,i,j],
P:zeromatrix(n,n),
for i:0 step 1 while i<(n-s) do(
  for j:0 step 1 while j<=i do(
    P[i+1,j+1]:binomial(i+j+2*s,j+s),
    P[j+1,i+1]:P[i+1,j+1]
  )
),
return(P)
)$

/*
This function generates the shifted Pascal matrix of the specified
 size, and raised to the specified power: Ps(n)^t.
Usage: pascals(n,s,t)
  n = size of the shifted Pascal matrix, n > 0.
  s = an integer, which is the shifting parameter.
  t = power of the Pascal matrix = any real number.

Note: The shifted Pascal matrix is a lower triangular matrix.
 To make an upper triangular shifted Pascal matrix, just transpose it.
 To make a symmetric shifted Pascal matrix, multiply the matrix
 by its transpose on the right: Ps.transpose(Ps)

Date: Mar 23, 2005
*/
pascals(n,s,t):=
block([Ps,i,j],
Ps:zeromatrix(n,n),
for i:0 step 1 while i<n do(
  for j:0 step 1 while j<=i do(
    Ps[i+1,j+1]:t^(i-j)*binomial(i+s,j+s)
  )
),
return(Ps)
)$

/*
This function generates a shift matrix, composed of zeros
 everywhere except for a line of 1's just BELOW the main
 diagonal, raised to a specified power.
 Ref: Lakshmikantham and Trigiante, "Theory of  Difference
 Equations", 2002, p17.

Usage: shiftl(n,t)
  n = size of the shift matrix, n > 0.
  t = power matrix is raised to, t <= n-1.

Note: If the matrix to shift is P and the shift matrix is K,
then the result of multiplication is:
 P.K = P shifted to left, with column of zeros on right.
 K.P = P shifted down, with row of zeros on top.
To shift right or up, use shiftu.

Date: Mar 24, 2005
*/
shiftl(n,t):=
block([K,j],
K:zeromatrix(n,n),
for j:1 thru n-t do K[j+t,j]:1,
return(K)
)$

/*
This function generates a shift matrix, composed of zeros
 everywhere except for a line of 1's just ABOVE the main
 diagonal, raised to a specified power.
 Ref: Lakshmikantham and Trigiante, "Theory of Difference
 Equations", 2002, p17.

Usage: shiftu(n,t)
  n = size of the shift matrix, n > 0
  t = power matrix is raised to, t <= n-1.

Note: If the matrix to shift is P and the shift matrix is K,
then the result of multiplication is:
 P.K = P shifted to right, with column of zeros on left.
 K.P = P shifted up, with row of zeros on bottom.
To shift left or down, use shiftl.

Date: Mar 24, 2005
*/
shiftu(n,t):=
block([K,i],
K:zeromatrix(n,n),
for i:1 thru n-t do K[i,i+t]:1,
return(K)
)$