\\% Circulant matrices, shift matrix.
\\ Author: Joerg Arndt
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2011-December-07 (12:05)


matcirc(v)=
\\ Example:
\\ matcirc([a,b,c,d]) ==
\\  [a b c d]
\\  [d a b c]
\\  [c d a b]
\\  [b c d a]
\\ Let a=[a0,a1,a2,a3], b=[b0,b1,b2,b3],
\\ R=matcirc(a), and C=R~, then
\\ R*b~ == (b*C)~ == correlation(a,b) ==
\\  [b0*a0 + (b1*a1 + (b2*a2 + b3*a3)),
\\   b1*a0 + (b2*a1 + (b3*a2 + b0*a3)),
\\   b2*a0 + (b3*a1 + (b0*a2 + b1*a3)),
\\   b3*a0 + (b0*a1 + (b1*a2 + b2*a3))]~
\\ and C*b~ == (b*R)~ == convolution(a,b) ==
\\  [b0*a0 + (b3*a1 + (b2*a2 + b1*a3)),
\\   b1*a0 + (b0*a1 + (b3*a2 + b2*a3)),
\\   b2*a0 + (b1*a1 + (b0*a2 + b3*a3)),
\\   b3*a0 + (b2*a1 + (b1*a2 + b0*a3))]~
{
    local(M, n);
    n = #v;
    M = matrix(n,n);
    for (r=1,n,
        for (c=r,n, M[r, c]=v[c-r+1] );  \\ upper right triangle and diagonal
        for (c=1,r-1, M[r, c]=v[n+c-r+1] );  \\ lower left triangle
    );
    return( M );
} /* ----- */

matcirc2(v)=
\\ Example:
\\ matcirc2([a,b,c,d]) ==
\\  [a b c d]
\\  [b c d c]
\\  [c d c b]
\\  [d c b a]
{
    local(M, n);
    n = #v;
    M = matrix(n,n);
    for (r=1,n,
        for (c=1,n+1-r, M[r, c]=v[r+c-1] );   \\ upper left triangle and diagonal
        for (c=n+2-r,n, M[r, c]=v[2*n+1-r-c] );  \\ lower right triangle
    );
    return( M );
} /* ----- */

mathankel(n, v)=
\\ Example:
\\ mathankel(4, [a,b,c,d,e,f,g,...]) ==
\\  [a b c d]
\\  [b c d e]
\\  [c d e f]
\\  [d e f g]
{
    local(M);
    M = matrix(n,n);
    for (r=1,n,
        for (c=1,n, M[r, c]=v[r+c-1] );
    );
    return( M );
} /* ----- */


matshift(n, s=1)=
\\ Example:
\\ matshift(4, 1) ==
\\  [0 1 0 0]
\\  [0 0 1 0]
\\  [0 0 0 1]
\\  [1 0 0 0]
\\ A * S == matrix A with rows shifted to the right
\\ S * A == matrix A with columns shifted up
\\ Let T = S~ ( = S^-1 ), then
\\ T * A == matrix A with columns shifted down
\\ A * T == matrix A with rows shifted to the left
{
    local(M);
    s = s % n;
    M = matrix(n,n);
    for (r=1, n-s,  M[r,r+s] = 1 );
    for (r=n-s+1, n,  M[r,r+s-n] = 1 );
    return( M );
} /* ----- */


matgetdiag(M, d=0)=
\\ Get shifted diagonal d of M (usual diagonal with d=0):
\\ return vector with elements of M from positions of
\\  nonzero elements of matshift( , d).
\\ Same as: diagonal of M*S^-d where S=matshift(n, 1)
{
    local(n, v);
    n = matsize(M)[1];
    d = d % n;
    v = vector(n);
    for (r=1, n-d,  v[r] = M[r,r+d] );
    for (r=n-d+1, n,  v[r] = M[r,r+d-n] );
    return( v );
} /* ----- */


matsetdiag(M, v, d=0)=
\\ Set shifted diagonal d of M (usual diagonal with d=0).
{
    local(n);
    n = matsize(M)[1];
    d = d % n;
    for (r=1, n-d,  M[r,r+d] = v[r] );
    for (r=n-d+1, n,  M[r,r+d-n] = v[r] );
    return( M );
} /* ----- */


\\ ==== end of file ====