![[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$