Return to ok_matrix.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

poly_lcm, poly_denominator(List), matrix_poly_to_matrix(Poly,Rule) are added.
Example. matrix_poly_to_matrix(x^2-1,[[x,newmat(2,2,[[a,0],[0,b]])]]);

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/ok_matrix.rr,v 1.8 2020/08/31 23:46:46 takayama Exp $ */
/* Old: Matrix */

/*-------------------------------*/
/* Package for matrix operations */
/*-------------------------------*/

#define _DEBUG_Matrix_ 0

/*---------*/
/* Utility */
/*---------*/
def omatrix_clone(M) {
  if (type(M) == 5) {
    N = size(M)[0];
    A = newvect(N);
    for (I=0; I<N; I++) {
      A[I] = omatrix_clone(M[I]);
    }
    return A;
  }else if (type(M) == 6) {
    P = size(M)[0];
    N = size(M)[1];
    A = newmat(P,N);
    for (I=0; I<P; I++) {
      for (J=0; J<N; J++) {
        A[I][J] = omatrix_clone(M[I][J]);
      }
    }
    return A;
  }else{
    return(M);
  }
}
def omatrix_ltom(L) {
  if (type(L) != 4 && type(L) != 5 && type(L) !=6) {
    error("parameter type error.");
  }
  if (type(L) == 5 || type(L) == 6) return L;
  if (type(L[0]) != 4) return omatrix_ltov(L);
  return newmat(length(L), length(L[0]), L);
}

def omatrix_mtol(M) {
  if (type(M) == 5) { /* vector case */
    return vtol(M);
  }
  if (type(M) != 6) {
    error("parameter type error.");
  }
  Size = size(M);
  L = [];
  for (I = 0; I < Size[0]; I++) {
    Row = [];
    for (J = 0; J < Size[1]; J++) {
      Row = append(Row, [M[I][J]]);
    }
    L = append(L, [Row]);
  }
  return L;
}

def omatrix_matrix_to_list(Y) {
  if ((type(Y) < 4) || (type(Y)>6)) return(Y);
  if (type(Y) == 4) return(map(omatrix_matrix_to_list,Y));
  if (type(Y) == 5) {
    Ans = newvect(length(Y));
    for (I=0; I<length(Y); I++) Ans[I]=omatrix_matrix_to_list(Y[I]);
    return(vtol(Ans));
  }
  if (type(Y) == 6) {
    V = newvect(size(Y)[1]);
    Ans=[];
    for (I=0; I<size(Y)[0]; I++) {
      for (J=0; J<size(Y)[1]; J++) {
        V[J] = omatrix_matrix_to_list(Y[I][J]);
      }
      V2 = vtol(V);
      Ans=cons(V2,Ans);
    }
    return(reverse(Ans));
  }
  error("Unknown type for omatrix_matrix_to_list");
}

def omatrix_ltov(L) {
  if (type(L) != 4) {
    error("parameter type error.");
  }
  return newvect(length(L), L);
}

/*-------------------*/
/* matrix operations */
/*-------------------*/
def omatrix_is_0(L) {
  if (type(L) == 5) {
    L = vtol(L);
  }
  if (type(L) != 4) {
    error("parameter type error.");
  }
  for (I = 0; I < length(L); I++) {
    if (L[I] != 0) {
      return 0;
    }
  }
  return 1;
}

def omatrix_1(N) {
  R = newmat(N, N);
  for (I = 0; I < N; I++) {
    R[I][I] = 1;
  }
  return R;
}

def omatrix_diag(L) {
  if (type(L) == 5) {
    L = vtol(L);
  }
  Size = length(L);

  R = newmat(Size, Size);
  for (I = 0; I < Size; I++) {
    R[I][I] = L[I];
  }
  return R;
}

def omatrix_mini(A, Row, Col) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }
  if (type(A) != 6) {
    error("parameter type is wrong.");
  }
  Size_A = size(A);
  if (Size_A[0] <= 1 || Size_A[1] <= 1) {
    error("a size of matrix A is too small.");
  }
  R = newmat(Size_A[0]-1, Size_A[1]-1);
  for (II = 0, I = 0; I < Size_A[0]; I++) {
    if (I != Row) {
      for (JJ = 0, J = 0; J < Size_A[1]; J++) {
	if (J != Col) {
	  R[II][JJ] = A[I][J];
	  JJ++;
	}
      }
      II++;
    }
  }
  return R;
}

/*&usage begin: omatrix_det2(A)
 It computes the determinant of  {A} by using the method of smaller determinant
 expansion.
 example: omatrix_det2([[1,x],[1,x^2]]);
end: */
def omatrix_det2(A) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }
  if (type(A) != 6) {
    error("parameter type is wrong.");
  }
  Size_A = size(A);
  if (Size_A[0] != Size_A[1]) {
    error("A is not NxN.");
  }
  Size = Size_A[0];

  if (Size == 1) {
    return A[0][0];
  }
  else {
    Det = 0;
    for (I = 0; I < Size; I++) {
      Det += (-1)^I*A[0][I]*omatrix_det2(omatrix_mini(A, 0, I));
    }
    return red(Det);
  }
}

/*&usage begin: omatrix_det(A)
 It computes the determinant of {A} by calling asir native det function.
 example: omatrix_det([[1,x],[1,x^2]]);
end: */
def omatrix_det(A) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }
  if (type(A) != 6) {
    error("parameter type is wrong.");
  }
  Size_A = size(A);
  if (Size_A[0] != Size_A[1]) {
    error("A is not NxN.");
  }
  Size = Size_A[0];

  B = newmat(Size,Size);
  D = 1;
  for (I=0; I<Size; I++) {
    for (J=0; J<Size; J++) {
      D = red(D*dn(A[I][J])/(gcd(D,dn(A[I][J]))));
    }
  }
  for (I=0; I<Size; I++) {
    for (J=0; J<Size; J++) {
      B[I][J] = red(A[I][J]*D);
    }
  }
  return red(det(B)/D^Size);
}

def omatrix_cofactor(A, Row, Col) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }
  Size_A = size(A);
  if (Size_A[0] != Size_A[1]) {
    error("A is not NxN.");
  }
  Size = Size_A[0];

  if (Size == 1) {
    return 1;
  }
  else {
    return (-1)^(Row + Col)*omatrix_det(omatrix_mini(A, Row, Col));
  }
}

def omatrix_adjugate(A) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }
  Size_A = size(A);
  if (Size_A[0] != Size_A[1]) {
    error("A is not NxN.");
  }
  Size = Size_A[0];

  R = newmat(Size, Size);
  for (I = 0; I < Size; I++) {
    for (J = 0; J < Size; J++) {
      R[I][J] = omatrix_cofactor(A, J, I);
    }
  }
  return R;
}

/*&usage begin: omatrix_inverse(A)
 It returns the inverse matrix of {A} by calling the asir native function
  invmat.  bug -- matrix of polynomial with complex coef, matrix of complex
  floating coef.
 example: omatrix_inverse([[1,1/(x+1)],[1,1/(x+1)^2]]);
end: */
def omatrix_inverse(A) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }
  if (type(A) != 6) {
    error("parameter type is wrong.");
  }
  Size_A = size(A);
  if (Size_A[0] != Size_A[1]) {
    error("A is not NxN.");
  }
  Size = Size_A[0];

  if (omatrix_float(A)) {
    C=invmat(A);
    CM = C[0]; CD=C[1];
    return CM*(1/CD);
  }

  B = newmat(Size,Size);
  D = 1;
  for (I=0; I<Size; I++) {
    for (J=0; J<Size; J++) {
      D = red(D*dn(A[I][J])/(gcd(D,dn(A[I][J]))));
    }
  }
  for (I=0; I<Size; I++) {
    for (J=0; J<Size; J++) {
      B[I][J] = red(A[I][J]*D);
    }
  }
  C = invmat(B);
  CM = C[0]; CD=C[1];
  DD = red(D/CD);
  return base_cancel(DD*CM);
}

/*&usage begin: omatrix_inverse2(A)
 It returns the inverse matrix of {A} by using adjunction matrix.
 example: omatrix_inverse2([[1,1/(x+1)],[1,1/(x+1)^2]]);
end: */
def omatrix_inverse2(A) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }
  Size_A = size(A);
  if (Size_A[0] != Size_A[1]) {
    error("A is not NxN.");
  }
  Size = Size_A[0];

  Det_A = omatrix_det(A);
  if (Det_A == 0) {
    error("Det(A) is zero.");
  }

  return omatrix_adjugate(A)/Det_A;
}

/*&usage begin: omatrix_inverse3(A)
 It computes the inverse matrix of {A} by using the Gaussian elimination.
 example: omatrix_inverse3([[1,1/(x+1)],[1,1/(x+1)^2]]);
end: */
/* Written by Nakayama, 2002 */
def omatrix_inverse3(A)
{
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }else {
    A = omatrix_clone(A);
  }
  Size_A = size(A);
  if (Size_A[0] != Size_A[1]) {
    error("A is not NxN.");
  }
  Size = Size_A[0];

  Det_A = omatrix_det(A);  /* BUG: should be removed. */
  if (Det_A == 0) {
    error("Det(A) is zero.");
  }
 B=A;
 N = size(B)[0];
 C = newmat(N, N);

 for (I = 0; I < N; I++)
  C[I][I] = 1;

 for (I = 0; I < N; I++) {
  J = I;
  while (J < N && B[J][I] == 0)
   J++;
   for (K = I; K < N; K++) {
     T = B[J][K];
     B[J][K] = B[I][K];
     B[I][K] = T;
   }
   for (K = 0; K < N; K++) {
     T = C[J][K];
     C[J][K] = C[I][K];
     C[I][K] = T;
   }
   V = B[I][I];
   B[I][I] = 1;
   for (K = I + 1; K < N; K++)
     B[I][K] = B[I][K] / V;
   for (K = 0; K < N; K++)
     C[I][K] = C[I][K] / V;
   for (J = 0; J < N; J++) {
     if (J == I)
       continue;
     V = B[J][I];
     B[J][I] = 0;
     for (K = I + 1; K < N; K++)
       B[J][K] = B[J][K] - B[I][K] * V;
     for (K = 0; K < N; K++)
       C[J][K] = C[J][K] - C[I][K] * V;
   }
 }
 return C;
}


def omatrix_trans(A) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }
  Size_A = size(A);
  Size_row = Size_A[1];
  Size_col = Size_A[0];

  R = newmat(Size_row, Size_col);
  for (I = 0; I < Size_row; I++) {
    for (J = 0; J < Size_col; J++) {
      R[I][J] = A[J][I];
    }
  }
  return R;
}

/* Solve A*X=Y */
def omatrix_solve(A, X, Y) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }
  if (type(X) == 4) {
    X = omatrix_ltov(X);
  }
  if (type(Y) == 4) {
    Y = omatrix_ltov(Y);
  }

  Size_A = size(A);
  Size_X = size(X);
  Size_Y = size(Y);

  /* ord(reverse(vtol(X))); */

  Sol = newvect(Size_X[0]);
  for (I = 0; I < Size_X[0]; I++) {
    Sol[I] = [X[I], X[I]];
  }

  Z = A*X - Y;
  Size_Z = size(Z);
  F = [];
  for (I = 0; I < Size_Z[0]; I++) {
    if (Z[I] != 0) {
      if (type(Z[I]) == 3) {
        F = append(F,[ nm(red(Z[I])) ]);
      }else{
        F = append(F, [Z[I]]);
      }
    }
  }

  if (F == []) {
    return vtol(Sol);
  }

  G = nd_gr(F, reverse(vtol(X)),0, 2);
#if _DEBUG_Matrix_
  print("G = ", 0);print(G);
#endif

  if (G[0] == 1 || G[0] == -1) {
    return [];
  }

  for (I = 0; I < length(G); I++) {
    for (J = Size_X[0] - 1; J >= 0; J--) {
      Nm = nm(G[I]); Dn = dn(G[I]);
      if (coef(Nm, 1, X[J]) != 0) {
	Sol[J] = [X[J], X[J] - Nm/coef(Nm, 1, X[J])];
	break;
      }
    }
  }

  return vtol(Sol);
}

/* A basis of Ker(A) */
def omatrix_kernel(A) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }

  Size_A = size(A);

  X = newvect(Size_A[1]);
  for (I = 0; I < Size_A[1]; I++) {
    X[I] = uc();
  }
  /* ord(reverse(vtol(X))); */

  Z = A*X;
  Size_Z = size(Z);
  F = [];
  for (I = 0; I < Size_Z[0]; I++) {
    if (Z[I] != 0) {
      F = append(F, [Z[I]]);
    }
  }

  if (F == []) {
    E = omatrix_1(Size_A[1]);
    Basis = newvect(Size_A[1]);
    for (I = 0; I < Size_A[1]; I++) {
      Basis[I] = E[I];
    }
    return append(size(Basis), vtol(Basis));
  }

  G = nd_gr(F, reverse(vtol(X)),0, 2);
#if _DEBUG_Matrix_
  print("G = ", 0);print(G);
#endif

  if (G[0] == 1 || G[0] == -1) {
    return [];
  }

  Sol = newvect(Size_A[1]);
  for (I = 0; I < Size_A[1]; I++) {
    Sol[I] = [X[I], X[I], 1];
  }

  for (I = 0; I < length(G); I++) {
    for (J = Size_A[1] - 1; J >= 0; J--) {
      Nm = nm(G[I]); Dn = dn(G[I]);
      if (coef(Nm, 1, X[J]) != 0) {
	Sol[J] = [X[J], coef(Nm, 1, X[J])/Dn*X[J] - G[I], coef(Nm, 1, X[J])/Dn];
	break;
      }
    }
  }
#if _DEBUG_Matrix_
  print("Sol = ",0);print(Sol);
#endif

  if (Size_A[1] - length(G) == 0) {
    return [0];
  }
  else {
    Dim = Size_A[1] - length(G);
    Basis = [];
    for (I = 0; I < Size_A[1]; I++) {
      E = [];
      for (J = 0; J < Size_A[1]; J++) {
	Sol_Nm = nm(Sol[J][1]); Sol_Dn = dn(Sol[J][1]);
	E = append(E, [red(coef(Sol_Nm,1,X[I])/Sol_Dn)/Sol[J][2]]);
      }
      if (!omatrix_is_0(E)) {
	Basis = append(Basis, [E]);
      }
    }
    if (Dim != length(Basis)) {
      error("algorithm error in omatrix_kernel().");
    }
    return append([Dim], [Basis]);
  }
}

def omatrix_eigenvalues(M) {
  if (type(M) == 4) {
    M = omatrix_ltom(M);
  }
  N = size(M)[0];
  L = newvect(N);
  for (I=0; I<N; I++) {
    L[I] = omatrix_v;
  }
  M2 = omatrix_diag(L);
  M2 = M-M2;
  L = omatrix_det(M2);
  return taka_poly_solve_poly1(L,omatrix_v|option_list=getopt());
}

def omatrix_rank(M) {
  if (type(M) == 4) {
    M = omatrix_ltom(M);
  }
  N = size(M)[1];
  A = omatrix_kernel(M);
  return N-A[0];
}

def omatrix_inner_product(A,B) {
  if (type(A) == 4) {
    N = length(A);
  }else{
    N = size(A)[0];
  }
  P = 0;
  for (I=0; I<N; I++) {
    P += A[I]*B[I];
  }
  return P;
}

def omatrix_submatrix(M,Ind) {
  if (type(M) == 6 || type(M) == 5) {
    Flag = 1;
    M = matrix_matrix_to_list(M);
  }
  R = [];
  N = length(Ind);
  for (I=N-1; I>=0; I--) {
    R = cons(M[Ind[I]],R);
  }
  /* Return the value in the same data type with the argument. */
  if (Flag) {
    if (R == []) return 0;
    return matrix_list_to_matrix(R);
  }else{
    return R;
  }
}

def omatrix_float(A) {
  if (type(A)<=1) {
    if ((ntype(A) == 1) || (ntype(A) == 3)) return 1;
    else return 0;
    /* BUG: case of complex numbers */
  }
  if (type(A) == 4) {
    for (I=0; I<length(A);  I++) {
      if (omatrix_float(A[I])) return 1;
    }
    return 0;
  }
  if (type(A) == 6 || type(A) == 5) {
    S = size(A);
    if (length(S) == 0) return 0;
    for (I=0; I<S[0]; I++) {
      if (omatrix_float(A[I])) return 1;
    }
    return 0;
  }
  return 0;
}

def omatrix_image(A) {
  if (type(A) == 4) {
    A = omatrix_ltom(A);
  }

  Size_A = size(A);

  X = newvect(Size_A[1]);
  for (I = 0; I < Size_A[1]; I++) {
    X[I] = uc();
  }
  /* ord(reverse(vtol(X))); */

  Z = A*X;
  Size_Z = size(Z);
  F = [];
  for (I = 0; I < Size_Z[0]; I++) {
    if (Z[I] != 0) {
      F = append(F, [Z[I]]);
    }
  }

  if (F == []) {
    return [[]];
  }

  G = nd_gr(F, reverse(vtol(X)),0, 2);
#if _DEBUG_Matrix_
  print("G = ", 0);print(G);
#endif

  if (G[0] == 1 || G[0] == -1) {
    error("omatrix_image: internal error."); 
    return [[]];
  }

  Sol = newvect(length(G));
  for (I = 0; I < length(Sol); I++) {
    T = newvect(length(X));
    for (J=0; J<length(X); J++) {
      T[J] = coef(G[I],1,X[J]);
    }
    Sol[I] = T;
  }
  return omatrix_mtol(map(omatrix_mtol,Sol));
}

def omatrix_gauge_transformation(M,T,V) {
  M=matrix_list_to_matrix(M);
  T=matrix_list_to_matrix(T);
  IT = omatrix_inverse(T);
  Ans=IT*M*T-IT*map(diff,T,V);
  return(map(red,Ans));
  
}

def ok_matrix_ij(N,II,JJ) {
  P=newmat(N,N);
  for (I=0; I<N; I++) {
    if (I == II) {
      P[I][JJ] = 1;
    }else if (I == JJ) {
      P[I][II] = 1;
    }else{
      P[I][I]=1;
    }
  }
  return(P);
}
def ok_matrix_zero_row(Mat) {
  D=size(Mat)[0];
  N=size(Mat)[1];
  for (I=0; I<D; I++) {
    Idx=I;
    for (J=0; J<N; J++) {
      if (Mat[I][J] != 0) {Idx=-1; break;}
    }
    if (Idx >= 0) return(Idx);
  }
  return(-1);
}
def ok_matrix_inverse_singular(Mat) {
  if (det(Mat) != 0) return(matrix_inverse(Mat));
  N=size(Mat)[0];
  Idx = ok_matrix_zero_row(Mat);
  if (Idx < 0) debug();
  P=ok_matrix_ij(N,Idx,N-1);
  Idx1=Idx; P1=P;
  Mat=matrix_transpose(P*Mat);
  Idx = ok_matrix_zero_row(Mat);
  if (Idx < 0) debug();
  P=ok_matrix_ij(N,Idx,N-1);
  Idx2=Idx; P2=P;
  Mat=matrix_transpose(P*Mat);

  Mat0=newmat(N-1,N-1);
  for (I=0; I<N-1; I++) for (J=0; J<N-1; J++) Mat0[I][J]=Mat[I][J];
  if (det(Mat0)==0) Inv0=ok_matrix_inverse_singular(Mat0);
  else Inv0=matrix_inverse(Mat0);
  Inv=newmat(N,N);
  for (I=0; I<N-1; I++) for (J=0; J<N-1; J++) Inv[I][J]=Inv0[I][J];

  Inv=P1*Inv*P2;
  return(Inv);
}

/*
  gauge.rr, matrix_replace() --> poly_to_matrix()
  E=matrix_identity_matrix(2)$ poly_to_matrix(5*a^2*b^3-3*b-5,[[a,A],[b,B]])-(5*A^2*B^3-3*B-5*E);
  Rule =[[scalar variable,matrix],...]
*/
def ok_poly_to_matrix(F,Rule) 
{
  V=[];
  Rule2=[];
  for (I=0; I<length(Rule); I++) {
    V=cons(Rule[I][0],V);
    Rule2=cons([Rule[I][0],T=matrix_list_to_matrix(Rule[I][1])],Rule2);
    N = size(T)[0];
  }
  V = reverse(V); 
  Rule = reverse(Rule2);
  F = dp_ptod(F,V);
  Ans=newmat(N,N);
  E = matrix_identity_matrix(N);
  while (F != 0) {
    T = dp_hc(F)*E;
    Exp = dp_etov(F);
    for (I=0; I<length(Exp); I++) {
      T = T*Rule[I][1]^Exp[I];
    }
    Ans = Ans + T;
    F = dp_rest(F);
  }
  return Ans;
}


Loaded_ok_matrix = 1 $
end$