![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mat CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / lib|
Return to mat CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/lib/mat,v 1.1.1.1 1999/11/10 08:12:30 noro Exp $ */
2: /* fraction free gaussian elimination; detructive */
3:
4: def deter(MAT)
5: {
6: S = size(MAT);
7: if ( car(S) != car(cdr(S)) )
8: return 0;
9: N = car(S);
10: for ( J = 0, D = 1; J < N; J++ ) {
11: for ( I = J; (I<N)&&(MAT[I][J] == 0); I++ );
12: if ( I != N ) {
13: for ( L = I; L < N; L++ )
14: if ( MAT[L][J] && (nmono(MAT[L][J]) < nmono(MAT[I][J])) )
15: I = L;
16: if ( J != I )
17: for ( K = 0; K < N; K++ ) {
18: B = MAT[J][K]; MAT[J][K] = MAT[I][K]; MAT[I][K] = B;
19: }
20:
21: for ( I = J + 1, V = MAT[J][J]; I < N; I++ )
22: for ( K = J + 1, U = MAT[I][J]; K < N; K++ )
23: MAT[I][K] = sdiv(MAT[I][K]*V-MAT[J][K]*U,D);
24: D = V;
25: } else
26: return ( 0 );
27: }
28: return (D);
29: }
30:
31: /* characteristic polynomial */
32:
33: def cp(M)
34: {
35: tstart;
36: S = size(M);
37: if ( car(S) != car(cdr(S)) )
38: return 0;
39: N = car(S);
40: MAT = newmat(N,N);
41: for ( I = 0; I < N; I++ )
42: for ( J = 0; J < N; J++ )
43: if ( I == J )
44: MAT[I][J] = red(M[I][J]-x);
45: else
46: MAT[I][J] = red(M[I][J]);
47: D = deter(MAT);
48: tstop;
49: return (D);
50: }
51:
52: /* calculation of charpoly by danilevskii's method */
53:
54: def da(MAT)
55: {
56: tstart;
57: S = size(MAT);
58: if ( car(S) != car(cdr(S)) )
59: return 0;
60: N = car(S);
61: M = newmat(N,N);
62: for ( I = 0; I < N; I++ )
63: for ( J = 0; J < N; J++ )
64: M[I][J] = red(MAT[I][J]);
65:
66: for ( W = newvect(N), J = 0, K = 0, D = 1; J < N; J++ ) {
67: for ( I = J + 1; (I<N) && (M[I][J] == 0); I++ );
68: if ( I == N ) {
69: for ( L = J, S = 1; L >= K; L-- )
70: S = S*x-M[L][J];
71: D *= S;
72: K = J + 1;
73: } else {
74: B = J + 1;
75: for ( L = 0; L < N; L++ ) {
76: T = M[I][L]; M[I][L] = M[B][L]; M[B][L] = T;
77: }
78: for ( L = 0; L < N; L++ ) {
79: T = M[L][B]; M[L][B] = M[L][I]; M[L][I] = T;
80: W[L] = M[L][J];
81: }
82: for ( L = K, T = red(1/M[B][J]); L < N; L++ )
83: M[B][L] *= T;
84: for ( L = K; L < N; L++ )
85: if ( W[L] && (L != J + 1) )
86: for ( B = K, T = W[L]; B < N; B++ )
87: M[L][B] -= M[J+1][B]*T;
88: for ( L = K; L < N; L++ ) {
89: for ( B = 0, T = 0; B < N ; B++ )
90: T += M[L][B] * W[B];
91: M[L][J + 1] = T;
92: }
93: }
94: }
95: tstop;
96: return ( D );
97: }
98:
99: def newvmat(N) {
100: M = newmat(N,N);
101: for ( I = 0; I < N; I++ )
102: for ( J = 0; J < N; J++ )
103: M[I][J] = strtov(rtostr(x)+rtostr(I))^J;
104: return M;
105: }
106:
107: def newssmat(N) {
108: M = newmat(N,N);
109: for ( I = 0; I < N; I++ )
110: for ( J = 0; J < N; J++ )
111: M[I][J] = strtov(rtostr(x)+rtostr(I)+"_"+rtostr(J));
112: return M;
113: }
114:
115: def newasssmat(N) {
116: N *= 2;
117: M = newmat(N,N);
118: for ( I = 0; I < N; I++ )
119: for ( J = 0; J < I; J++ )
120: M[I][J] = strtov(rtostr(x)+rtostr(I)+"_"+rtostr(J));
121: for ( I = 0; I < N; I++ )
122: for ( J = I + 1; J < N; J++ )
123: M[I][J] = -M[J][I];
124: return M;
125: }
126:
127: /* calculation of determinant by minor expansion */
128:
129: def edet(M) {
130: S = size(M);
131: if ( S[0] == 1 )
132: return M[0][0];
133: else {
134: N = S[0];
135: L = newmat(N-1,N-1);
136: for ( I = 0, R = 0; I < N; I++ ) {
137: for ( J = 1; J < N; J++ ) {
138: for ( K = 0; K < I; K++ )
139: L[J-1][K] = M[J][K];
140: for ( K = I+1; K < N; K++ )
141: L[J-1][K-1] = M[J][K];
142: }
143: R += (-1)^I*edet(L)*M[0][I];
144: }
145: return R;
146: }
147: }
148:
149: /* sylvester's matrix */
150:
151: def syl(V,P1,P2) {
152: D1 = deg(P1,V); D2 = deg(P2,V);
153: M = newmat(D1+D2,D1+D2);
154: for ( J = 0; J <= D2; J++ )
155: M[0][J] = coef(P2,D2-J,V);
156: for ( I = 1; I < D1; I++ )
157: for ( J = 0; J <= D2; J++ )
158: M[I][I+J] = M[0][J];
159: for ( J = 0; J <= D1; J++ )
160: M[D1][J] = coef(P1,D1-J,V);
161: for ( I = 1; I < D2; I++ )
162: for ( J = 0; J <= D1; J++ )
163: M[D1+I][I+J] = M[D1][J];
164: return M;
165: }
166:
167: /* calculation of resultant by edet() */
168:
169: def res_minor(V,P1,P2)
170: {
171: D1 = deg(P1,V); D2 = deg(P2,V);
172: M = newmat(D1+D2,D1+D2);
173: for ( J = 0; J <= D2; J++ )
174: M[0][J] = coef(P2,D2-J,V);
175: for ( I = 1; I < D1; I++ )
176: for ( J = 0; J <= D2; J++ )
177: M[I][I+J] = M[0][J];
178: for ( J = 0; J <= D1; J++ )
179: M[D1][J] = coef(P1,D1-J,V);
180: for ( I = 1; I < D2; I++ )
181: for ( J = 0; J <= D1; J++ )
182: M[D1+I][I+J] = M[D1][J];
183: return edet(M);
184: }
185: end$