/* -*- gp-script -*- */ \\% Circulant matrix, Toeplitz matrix, Hankel matrix, shift matrix. \\ Author: Joerg Arndt \\ License: GPL version 3 or later \\ online at http://www.jjj.de/pari/ \\ version: 2014-October-16 (18:31) matcirc(v)= \\ Return circulant matric corresponding to cyclic correlation. \\ The transpose is the circulant matric corresponding to cyclic convolution. \\ 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))]~ { my(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] { my(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 ); } /* ----- */ mattoeplitz(v)= \\ Return n x n Toeplitz matrix, v[] must have 2*n+1 entries. \\ Example: \\ mattoeplitz([A,B,C,D,e,f,g]) == \\ [A B C D] \\ [e A B C] \\ [f e A B] \\ [g f e A] { my(n,M); if ( (#v+1)%2 != 0, error("Need an odd number of elements")); n = (#v + 1) / 2; M = matrix(n,n); for (d=1, n, M[1,d] = v[d] ); for (r=2, n, for (c=r,n, M[r,c] = M[r-1,c-1] ); ); for (d=2, n, M[d,1] = v[n-1+d] ); for (r=3, n, for (c=2,r-1, M[r,c] = M[r-1,c-1] ); ); return( M ); } /* ----- */ mattoeplitz2(v)= \\ Return n x n Toeplitz matrix, v[] must have 2*n+1 entries. \\ Example: \\ mattoeplitz([A,B,C,D,e,f,g]) == \\ [A B C D] \\ [g A B C] \\ [f g A B] \\ [e f g A] { if ( (#v+1)%2 != 0, error("Need an odd number of elements")); my(n, nh, t, a, b); n = (#v + 1) / 2; nh = (n-1) \ 2; for (j=1, nh, a=n+j; b=#v+1-j; t=v[a]; v[a]=v[b]; v[b]=t ); return( mattoeplitz(v) ); } /* ----- */ 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] { my(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 { my(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) { my(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 cyclically shifted diagonal d of M (usual diagonal with d=0). { my(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 ====