Return to noro_matrix.rr CVS log

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

Fixed a bug in linalg.minipoly_check_mat.

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/noro_matrix.rr,v 1.16 2017/04/27 22:41:33 noro Exp $ */

/* functions for linear algebra */

load("sp")$

module linalg;

localf random_mat,random_smat,random_unimodular_mat,random_rmat,random_vect;
localf unit_mat,subst_mat,mulsub_mat,dup_mat,ptozp_mat,zero_mat;
localf mulsub_vect,muladd_vect,dup_vect,mul_smat_vect,mul_svect_vect;
localf mat_to_smat,smat_to_mat,vect_to_svect,zero_vect,ptozp_vect;
localf etrans_mat,swap_mat,sample_mat;
localf reduced_form_mat,mat_to_col;
localf companion_mat;

localf minipoly_mat,minipoly_mat_direct,reduce_mat;
localf minipoly_cand_mat,minipoly_cand_smat;
localf minipoly_check_mat,minipoly_check_smat;
localf minipoly_mat_direct;

localf reduce_vect,reduce_vect_hist,insert_to_base_hist,insert_to_base;
localf compute_kernel,compute_image,create_basis;
localf const_term,get_alg,is_integer;

localf jordan_canonical_form,compute_jordan_block_basis,jordan_decomp;
localf partial_fraction;
localf print_kernel,print_vect;

localf next_combination,reduce_symmetric_form,generate_elementary_symmetric_form;

localf extended_euclid2,extended_euclid;
localf binomial_coef;
localf power_mat,exp_mat;

def companion_mat(F,V)
{
	D = deg(F,V);
	A = matrix(D,D);
	for ( I = 0; I < D-1;I++ )
		A[I][I+1] = 1;
	for ( J = 0; J < D; J++ )
		A[D-1][J] = -coef(F,J,V);
	return A;
}

def random_mat(N,Bound)
{
	A = matrix(N,N);
	for ( I = 0; I < N; I++ )
		for ( J = 0; J < N; J++ )
			A[I][J] = (random()%Bound)*(random()%2?-1:1);
	return A;
}

def random_rmat(M,N,Bound)
{
	A = matrix(M,N);
	for ( I = 0; I < M; I++ )
		for ( J = 0; J < N; J++ )
			A[I][J] = (random()%Bound)*(random()%2?-1:1);
	return A;
}

def random_smat(N,Bound,Percent)
{
	R = [];
	for ( I = N-1; I >= 0; I-- ) {
		V = [];
		for ( J = N-1; J >= 0; J-- ) {
			if ( (random()%100) < Percent ) {
				C = (random()%Bound)*(random()%2?-1:1);
				V = cons([J,C],V);
			}
		}
		R = cons([I,V],R);
	}
	return cons([N,N],R);
}

def etrans_mat(A,I,J,C,Left)
{
	M = size(A)[0];
	N = size(A)[1];
	if ( Left ) {
		for ( K = 0; K < N; K++ )
			A[I][K] += C*A[J][K];
	} else {
		for ( K = 0; K < M; K++ )
			A[K][J] += C*A[K][I];
	}
}

def swap_mat(A,I,J,Row)
{
	M = size(A)[0];
	N = size(A)[1];
	if ( Row ) {
		for ( K = 0; K < N; K++ ) {
			T = A[I][K]; A[I][K] = A[J][K]; A[J][K] = T;
		}
	} else {
		for ( K = 0; K < M; K++ ) {
			T = A[K][I]; A[K][I] = A[K][J]; A[K][J] = T;
		}
	}
}

def random_unimodular_mat(N)
{
	P = unit_mat(N);
	PI = unit_mat(N);
	for ( I = 0; I < N; I++ ) {
		K = random()%N;
		while ( 1 ) {
			L = random()%N;
			if ( K != L ) break;
		}
		C = random()%10;
		Left = random()%2;
		etrans_mat(P,K,L,C,Left);
		etrans_mat(PI,K,L,-C,!Left);
		K = random()%N;
		while ( 1 ) {
			L = random()%N;
			if ( K != L ) break;
		}
		swap_mat(P,K,L,1);
		swap_mat(PI,K,L,0);
	}
	return [P,PI];
}

def subst_mat(P,A)
{
	N = size(A)[0];
	if ( type(P) <= 1 )
		return P*unit_mat(N);
	else {
		V = var(P);
		QR = sqr(P,V);
		/* P = Q*V+R */
		return subst_mat(QR[0],A)*A+QR[1]*unit_mat(N);
	}
}

/* A = A*C1-B*C2 */

def mulsub_mat(A,C1,B,C2)
{
	S = size(A); R = S[0]; C = S[1];
	for( I = 0; I < R; I++ )
		for ( J = 0; J < C; J++ )
			A[I][J] = A[I][J]*C1-B[I][J]*C2;
}

def dup_mat(A)
{
	S = size(A); R = S[0]; C = S[1];
	M = matrix(R,C);
	for( I = 0; I < R; I++ )
		for ( J = 0; J < C; J++ )
			M[I][J] = A[I][J];
	return M;
}

def unit_mat(N)
{
	M = matrix(N,N);
	for ( I = 0; I < N; I++ ) M[I][I] = 1;
	return M;
}

/* Base = [[Mat,I,J,Poly],...] */

def minipoly_mat_direct(A)
{
	L = ptozp_mat(A);
	/* A = B/DN */
	B = L[0]; DN = L[1];
	N = size(B)[0];
	E = unit_mat(N);
	Base = [[E,0,0,1]];
	Bk = E;
	for ( K = 1; ; K++ ) {
		Bk *= B;	
		R = reduce_mat(Bk,x^K,Base);
		if ( !R[0] ) {
			P = ptozp(R[3]);
			if ( coef(P,deg(P,x)) < 0 )
				P = -P;
			break;
		} else
			Base = insert_to_base_hist(R,Base);
	}
	/* m_B(x)=P  => m_A(x)=m_B(DN*x) */
	return ptozp(subst(P,x,DN*x));
}

def minipoly_mat(A)
{
	L = ptozp_mat(A);
	/* A = B/DN */
	B = L[0]; DN = L[1];
	while ( 1 ) {
		P = minipoly_cand_mat(B);
		if ( P && minipoly_check_mat(B,P) ) break;
	}
	/* m_B(x)=P  => m_A(x)=m_B(DN*x) */
	return ptozp(subst(P,x,DN*x));
}

def reduce_mat(M,P,Base)
{
	N = size(M)[0];
	M = dup_mat(M);
	for ( T = Base; T != []; T = cdr(T) ) {	
		B = car(T);
		BM = B[0]; I = B[1]; J = B[2]; BP = B[3];
		if ( M[I][J] ) {
			/* M <- M*BM[I][J]-M[I][J]*BM */
			/* P <- P*BM[I][J]-M[I][J]*BP */
			C1 = BM[I][J]; C2 = M[I][J];
			G = igcd(C1,C2);
			C1 = idiv(C1,G);
			C2 = idiv(C2,G);
			mulsub_mat(M,C1,BM,C2);
			P = P*C1-BP*C2;
		}
	}
	for ( I = 0; I < N; I++ )
		for ( J = 0; J < N; J++ )
			if ( M[I][J] ) return [M,I,J,P];
	return [0,-1,-1,P];
}

def insert_to_base_hist(R,Base)
{
	if ( Base == [] )
		return [R];
	Top = car(Base);
	if ( R[1] < Top[1] || (R[1] == Top[1] && R[2] < Top[2]) )
		return cons(R,Base);
	else
		return cons(Top,insert_to_base_hist(R,cdr(Base)));
}

/* sparse matrix */

def mat_to_smat(A)
{
	N = size(A)[0];
	R = [];
	for ( I = N-1; I >= 0; I-- ) {
		V = vect_to_svect(A[I]);
		R = cons([I,V],R);
	}
	return cons(size(A),R);
}

def smat_to_mat(A)
{
	Size = A[0];
	R = matrix(Size[0],Size[1]);
	for ( T = cdr(A); T != []; T = cdr(T) ) {
		P = car(T);
		I = P[0]; V = P[1];
		for ( S = V; S != []; S = cdr(S) ) {
			P = car(S);
			R[I][P[0]] = P[1];
		}
	}
	return R;
}

def vect_to_svect(V)
{
	N = size(V)[0];
	R = [];
	for ( I = N-1; I >= 0; I-- )
		if ( V[I] )
			R = cons([I,V[I]],R);
	return R;
}

def random_vect(N,Bound)
{
	A = vector(N);
	for ( I = 0; I < N; I++ )
		A[I] = (random()%Bound)*(random()%2?-1:1);
	return A;
}

/* A = A*C1-B*C2 */

def mulsub_vect(A,C1,B,C2)
{
	S = size(A); R = S[0];
	for( I = 0; I < R; I++ )
		A[I] = A[I]*C1-B[I]*C2;
}

/* A += B*C */

def muladd_vect(A,B,C)
{
	S = size(A); R = S[0];
	for( I = 0; I < R; I++ )
		A[I] += B[I]*C;
}

def dup_vect(A)
{
	S = size(A); R = S[0];
	M = vector(R);
	for( I = 0; I < R; I++ )
		M[I] = A[I];
	return M;
}

/* A = [[R,C],[I,[[J,AIJ],...]],...] */

def mul_smat_vect(A,V)
{
	N = size(V)[0];
	Row = A[0][0];
	R = vector(Row);
	for ( T = cdr(A), I = 0; T != []; T = cdr(T), I++ )
		R[I] = mul_svect_vect(car(T)[1],V);
	return R;
}

/* A = [[J,AIJ],...] */

def mul_svect_vect(A,V)
{
	R = 0;
	for ( T = A; T != []; T = cdr(T) ) {
		S = car(T);
		R += S[1]*V[S[0]];
	}
	return R;
}

/* Base = [[Vect,I,Poly],...] */
def minipoly_cand_mat(A)
{
	N = size(A)[0];
	V = random_vect(N,50);
	Base = [[V,0,1]];
	AkV = V;
	for ( K = 1; ; K++ ) {
		AkV = A*AkV;
			R = reduce_vect_hist(AkV,x^K,Base);
			if ( !R[0] ) {
				P = ptozp(R[2]);
				if ( coef(P,deg(P,x)) < 0 )
					P = -P;
				return P;
			} else
			Base = insert_to_base(R,Base);
	}
}

/* Base = [[Vect,I,Poly],...] */
def minipoly_cand_smat(A)
{
	N = A[0][0];
	V = random_vect(N,50);
	Base = [[V,0,1]];
	AkV = V;
	for ( K = 1; ; K++ ) {
		AkV = mul_smat_vect(A,AkV);
		R = reduce_vect_hist(AkV,x^K,Base);
		if ( !R[0] ) {
			P = ptozp(R[2]);
			if ( coef(P,deg(P,x)) < 0 )
				P = -P;
			return P;
		} else
			Base = insert_to_base(R,Base);
	}
}

/* reduce a vector M by a sorted base Base.
   P is a symbol representing M.
   returns [M',Head,P'], where Head is the
   position of the leftmost nonzero element in M',
   M' is the reduced vector, P' represents M'.
*/

def reduce_vect_hist(M,P,Base)
{
	N = size(M)[0];
	M = dup_vect(M);
	for ( T = Base; T != []; T = cdr(T) ) {	
		B = car(T);
		BM = B[0]; I = B[1]; BP = B[2];
		if ( M[I] ) {
			/* M <- M*BM[I]-M[I]*BM */
			/* P <- P*BM[I]-M[I]*BP */
			C1 = BM[I]; C2 = M[I];
			G = igcd(C1,C2);
			C1 = idiv(C1,G);
			C2 = idiv(C2,G);
			mulsub_vect(M,C1,BM,C2);
			P = P*C1-BP*C2;
		}
	}
	for ( I = 0; I < N; I++ )
		if ( M[I] ) return [M,I,P];
	return [0,-1,P];
}

def insert_to_base(R,Base)
{
	if ( Base == [] )
		return [R];
	Top = car(Base);
	if ( R[1] < Top[1] )
		return cons(R,Base);
	else
		return cons(Top,insert_to_base(R,cdr(Base)));
}

def minipoly_check_mat(S,P)
{
  N = size(S)[0];
	B = subst_mat(P,S);
	for ( I = 0; I < N; I++ )
		for ( J = 0; J < N; J++ )
			if ( B[I][J] ) return 0;
	return 1;
}

def minipoly_check_smat(S,P)
{
	Size = S[0];
	N = Size[0];
	V = vector(N);
	X = var(P);
	D = deg(P,X);
	C = vector(D+1);
	for ( I = 0; I <= D; I++ )
		C[I] = coef(P,I,X);
	for ( I = 0; I < N; I++ ) {
		for ( J = 0; J < N; J++ )
			V[J] = J==I?1:0;
		AkV = dup_vect(V);
		V *= C[0];
		for ( K = 1; K <= D; K++ ) {
			AkV = mul_smat_vect(S,AkV);
			muladd_vect(V,AkV,C[K]);
		}
		for ( J = 0; J < N; J++ )
			if ( V[J] ) break;
		if ( J < N ) return 0;
	}
	return 1;
}

/* A = 1/D*B => return [B,D] */

def ptozp_mat(A)
{
	Size = size(A);
	M = Size[0]; N = Size[1];
	D = 0; K = 0;
	for ( I = 0; I < M; I++ )
		for ( J = 0; J < N; J++ ) {
			D += A[I][J]*<<K>>; K++;
		}
	if ( !D ) {
		B = matrix(M,N);
		return [B,1];
	} else {
		E = dp_ptozp(D);
		C = dp_hc(E)/dp_hc(D);
		return [C*A,C];
	}
}

def ptozp_vect(A)
{
	M = size(A)[0];
	D = 0; K = 0;
	for ( I = 0; I < M; I++ ) {
		D += A[I]*<<K>>; K++;
	}
	if ( !D ) {
		B = vector(M);
		return [B,1];
	} else {
		E = dp_ptozp(D);
		C = dp_hc(E)/dp_hc(D);
		return [C*A,C];
	}
}

def create_basis(B)
{
	N = length(B);
	Base = [];
	for ( K = 0; K < N; K++ ) {
		AkV = ptozp_vect(B[K])[0];
		R = reduce_vect_hist(AkV,<<K>>,Base);
		if ( R[0] )
			Base = insert_to_base(R,Base);
	}
	return Base;
}

def sample_mat(Data)
{
	Asis = getopt(asis);
	/* Data = [[Alpha,M],...] */
	N = 0;
	for ( T = Data; T != []; T = cdr(T) )
		N += car(T)[1];
	A = matrix(N,N);
	for ( T = Data, I = 0; T != []; T = cdr(T) ) {
		B = car(T);
		Alpha = B[0];
		M = B[1];
		for ( J = 0; J < M; J++ ) A[I+J][I+J] = Alpha;
		for ( J = 1; J < M; J++ ) A[I+J-1][I+J] = 1;
		I += M;
	}
	if ( type(Asis) != -1 )
		return A;

	/* L = [P,P^-1] */
	L = random_unimodular_mat(N);
	return L[1]*A*L[0];
}

def jordan_decomp(A)
{
	N = size(A)[0];
	M = minipoly_mat(A);
	FM = fctr(M);
	/* skip const */
	FM = cdr(FM);
	L = length(FM);

	/* extracting roots and multiplicities */
	Alpha = vector(L);
	E = vector(L);
	for ( I = 0; I < L; I++ )
		if ( deg(F=FM[I][0],x) != 1 ) 
			error("irrational eigenvalue : not implemented yet");
		else {
			Alpha[I] = -coef(F,0,x)/coef(F,1,x);
			E[I] = FM[I][1];
		}
	F = vector(L);
	M = sdiv(M,coef(M,deg(M,x)));
	for ( I = 0; I < L; I++ )
		F[I] = sdiv(M,(x-Alpha[I])^E[I]);

	/* computation of coefficient polys */
	T = extended_euclid(F);
	G = T[0]; GCD = T[1];
	G /= GCD;
	for ( I = 0; I < L; I++ )
		G[I] = srem(G[I],(x-Alpha[I])^E[I]);

	/* projection operators */
	P = vector(L);
	for ( I = 0; I < L; I++ )
		P[I] = subst_mat(G[I]*F[I],A);

	/* computation of S, N */
	S = matrix(N,N);
	for ( I = 0; I < N; I++ )
		for ( J = 0; J < N; J++ ) {
			for ( T = 0, K = 0; K < L; K++ )
				T += Alpha[K]*P[K][I][J];
			S[I][J] = T;
		 }
	return [S,A-S];
}

def binomial_coef(M,I)
{
	R = 1;
	for ( J = 0; J < I; J++ ) R *= M-J;
	return R/fac(I);
}

def power_mat(A)
{
	Solve = getopt(solve);
	if ( type(Solve)== -1 ) Solve = 0;

	/* JCF = [P,Block,DefIdeal] */
	JCF = jordan_canonical_form(A);
	P = JCF[0];
	Block = JCF[1];
	DefIdeal = JCF[2];
	if ( !Solve ) {
		Inv = invmat(P); PI = map(red,Inv[0]/Inv[1]);
	}
	N = size(A)[0];
	PAPM = matrix(N,N);
	for ( T = Block, I = 0; T != []; T = cdr(T) ) {
		BI = car(T); Alpha = BI[0]; K = BI[1]; R = BI[2];
		JM = matrix(K,K);
		for ( Y = 0; Y < K; Y++ ) {
			Ent = (Alpha==1?1:Alpha^(m-Y))*binomial_coef(m,Y);
			for ( X = 0; X+Y < K; X++ )
				JM[X][X+Y] = Ent;
		}
		for ( J = 0; J < R; J++, I += K )
			for ( X = 0; X < K; X++ )
				for ( Y = X; Y < K; Y++ )
					PAPM[I+X][I+Y] = JM[X][Y];
	}
	if ( Solve ) {
		S = map(red,P*PAPM);
		for ( J = N-1, Basis = []; J >= 0; J-- ) {
			SI = vector(N);
			for ( I = 0; I < N; I++ )
				SI[I] = S[I][J];
			Basis = cons(SI,Basis);
		}
		return Basis;
	} else
		return map(red,P*map(red,PAPM*PI));
}

def exp_mat(A)
{
	Solve = getopt(solve);
	if ( type(Solve)== -1 ) Solve = 0;

	/* JCF = [P,Block,DefIdeal] */
	JCF = jordan_canonical_form(A);
	P = JCF[0];
	Block = JCF[1];
	DefIdeal = JCF[2];
	if ( !Solve ) {
		Inv = invmat(P); PI = map(red,Inv[0]/Inv[1]);
	}
	N = size(A)[0];
	ExpJ = matrix(N,N);
	for ( T = Block, I = 0; T != []; T = cdr(T) ) {
		BI = car(T); Alpha = BI[0]; K = BI[1]; R = BI[2];
		JM = matrix(K,K);
		Exp = Alpha==0?1:exp(Alpha*t);
		for ( Y = 0; Y < K; Y++ ) {
			Ent = Exp*t^Y/fac(Y);
			for ( X = 0; X+Y < K; X++ )
				JM[X][X+Y] = Ent;
		}
		for ( J = 0; J < R; J++, I += K )
			for ( X = 0; X < K; X++ )
				for ( Y = X; Y < K; Y++ )
					ExpJ[I+X][I+Y] = JM[X][Y];
	}
	if ( Solve ) {
		S = map(red,P*ExpJ);
		for ( J = N-1, Basis = []; J >= 0; J-- ) {
			SI = vector(N);
			for ( I = 0; I < N; I++ )
				SI[I] = S[I][J];
			Basis = cons(SI,Basis);
		}
		return Basis;
	} else
	return map(red,P*map(red,ExpJ*PI));
}

def compute_image(A)
{
	B = mat_to_col(A);
	return create_basis(B);
}

/* 
 *	Mat : matrix of size MxN
 *	Index : vector of length M
 *  Index[I] = head position of Mat[I]
 */

def reduced_form_mat(Mat,Index,M,N) 
{
	for ( J = 0, L = 0, D = 1; J < N; J++ ) {
		for ( I = L; I < M && !Mat[I][J]; I++ );
		if ( I == M )
			continue;
		Index[L] = J;
		for ( K = 0; K < N; K++ ) {
			T = Mat[I][K]; Mat[I][K] = Mat[L][K]; Mat[L][K] = T;
		}
		for ( I = L + 1, V = Mat[L][J]; I < M; I++ )
			for ( K = J, U = Mat[I][J]; K < N; K++ )
				Mat[I][K] = sdiv(Mat[I][K]*V-Mat[L][K]*U,D);
		D = V; L++;
	}
	for ( I = L; I < M; I++ )
		for ( J = 0; J < N; J++ )
			if ( Mat[I][J] )
				return [-1,0];
	for ( I = L - 2, W = vector(N); I >= 0; I-- ) {
		for ( J = 0; J < N; J++ ) W[J] = 0;
		for ( G = I + 1; G < L; G++ )
			for ( H = Index[G], U = Mat[I][H]; H < N; H++ )
				W[H] += Mat[G][H]*U;
		for ( J = Index[I], U = Mat[I][J]; J < N; J++ )
			Mat[I][J] = sdiv(Mat[I][J]*D-W[J],U);
	}
	return [L,D];
}

def mat_to_col(A)
{
	R = [];
	N = size(A)[0];
	for ( J = N-1; J >= 0; J-- ) {
		T = vector(N);
		for ( I = 0; I < N; I++ )
			T[I] = A[I][J];
		R = cons(T,R);
	}
	return R;
}

def zero_vect(A)
{
	S = size(A);
	M = S[0];
	for ( I = 0; I < M; I++ )
		if ( A[I] ) return 0;
	return 1;
}

def zero_mat(A)
{
	S = size(A);
	M = S[0]; N = S[1];
	for ( I = 0; I < M; I++ )
		for ( J = 0; J < N; J++ )
			if ( A[I][J] ) return 0;
	return 1;
}

def extended_euclid2(F,G)
{
	F1 = F; F2 = G;
	A1 = 1; A2 = 0;
	B1 = 0; B2 = 1;
	while ( 1 ) {
		Q = sdiv(F1,F2);
		F3 = F1-Q*F2;
		if ( !F3 ) {
			return [A2,B2,F2];
		}
		A3 = A1-Q*A2;
		B3 = B1-Q*B2;
		F1 = F2; F2 = F3;
		A1 = A2; A2 = A3;
		B1 = B2; B2 = B3;
	}
}

def extended_euclid(F)
{
	N = size(F)[0];
	G = vector(N);
	GCD = F[N-1];
	G[N-1]=1;
	for ( I = N-2; I >= 0; I-- ) {
		/* G[I+1]*F[I+1]+...+G[N-1]*F[N-1]=GCD */
		/* C[0]*F[I]+C[1]*GCD=C[2] */
		C = extended_euclid2(F[I],GCD);
		if ( C[0]*F[I]+C[1]*GCD-C[2] )
			error("afo");
		/* C[0]*F[I]+C[1]*G[I+1]*F[I+1]+...+C[1]*G[N-1]*F[N-1]=C[2] */
		G[I] = C[0];
		for ( J = I+1; J < N; J++ )
			G[J] *= C[1];
		GCD = C[2];

		for ( K = I, T = 0; K < N; K++ )
			T += G[K]*F[K];
		if ( T != GCD )
			error("bfo");
	}
	return [G,GCD];
}


def const_term(F)
{
	V = var(F);
	if ( V )
		return const_term(coef(F,0,V));
	else
		return F;
}

def get_alg(L)
{
	for ( T = L, R = []; T != []; T = cdr(T) )
		R = union_sort(R,getalgtreep(car(T)));
	return R;
}

def compute_kernel(A)
{
	B = getopt(rhs);
	if ( type(B) == -1 || (type(B)==5 && zero_vect(B)) ) B = 0;
	S = size(A); M = S[0]; N = S[1];
	for ( I = 0, V = vector(N); I < N; I++ ) V[I] = strtov("x"+rtostr(I));
	VL = reverse(vtol(V));
	E = B ? vtol(A*V-B) : vtol(A*V);
	Alg = get_alg(E);
	if ( Alg != [] ) {
		AlgV = map(algptorat,Alg);
		Defpoly = map(defpoly,Alg);
		Eext = append(reverse(Defpoly),map(algptorat,E));
		VLext = append(VL,AlgV);
		Ord = [[0,N],[0,length(Alg)]];
		G = nd_gr(Eext,VLext,0,Ord);
		G = map(rattoalgp,G,Alg);
		for ( T = G, G0 = []; T != []; T = cdr(T) )
			if ( vars(car(T)) != [] ) G0 = cons(car(T),G0);
		G = reverse(G0);
	} else {
		dp_ord(0);
		G = nd_gr(E,VL,0,0);
	}
	D = map(dp_ptod,G,VL);
	for ( LCM = 1, T = D; T != []; T = cdr(T) )
		LCM = ilcm(dp_hc(car(T)),LCM);
	for ( T = D, Sol = LCM*V; T != []; T = cdr(T) ) {
		P = car(T);
		EV = dp_etov(dp_ht(P));
		for ( I = 0; !EV[I]; I++ );
		Sol[N-I-1] = dp_dtop((-LCM/dp_hc(P))*dp_rest(P),VL);
	}
	VSol = vars(Sol);
	for ( T = [], I = 0; I < N; I++ )
		if ( member(V[I],VSol) ) T = cons(V[I],T);
	for ( Ker = []; T != []; T = cdr(T) ) {
		Sol0 = map(coef,Sol,1,car(T));
		if ( Alg == [] )
			Sol0 = ptozp_vect(Sol0)[0];
		for ( I = 0; I < N && !Sol0[I]; I++ );
		Ker = cons([Sol0,I],Ker);
	}
	return B ? [map(const_term,Sol)/LCM,Ker] : Ker;
		
}

def jordan_canonical_form(A)
{
	Verbose = getopt(verbose);
	Verbose = type(getopt(verbose))==-1?0:1;

	N = size(A)[0];
	M = minipoly_mat(A);
	if ( type(Ext=getopt(ext)) != -1 ) {
		FM = af(M,Ext);
	} else {
		FM = fctr(M);
	}
	/* skip const */
	FM = cdr(FM);
	L = length(FM);
	if ( Verbose ) {
		print("minimal polynomial = ",0);
		for ( I = 0; I < L; I++ ) {
			F = FM[I][0];
			E = FM[I][1];
			print("(",0); print(F,0); print(")",0);
			if ( E > 1 )
				print("^",0); print(E,0);
		}
		print("");
	}
	/* extracting roots and multiplicities */
	Alpha = vector(L);
	E = vector(L);
	D = vector(L);
	for ( I = 0; I < L; I++ ) {
		F = FM[I][0];
		D[I] = deg(F,x);
		Alpha[I] = D[I]==1 ? -coef(F,0,x)/coef(F,1,x) : newalg(F);
		E[I] = FM[I][1];
	}

	Q = matrix(N,N);
	Prefix = strtoascii("a")[0];
	DefIdeal = [];
	JCF = [];
	for ( K = 0, J = 0; K < L; K++ ) {
		if ( Verbose ) {
			print("computing jordan blocks for Alpha = ",0); print(Alpha[K]);
		}
		/* JBK = [P,BlockInfo] */
		JBK = compute_jordan_block_basis(A,Alpha[K],E[K],Verbose);
		if ( (DK = D[K]) == 1 ) {
			for ( T = JBK[0]; T != []; T = cdr(T), J++ ) {
				V = car(T);
				for ( I = 0; I < N; I++ )
					Q[I][J] = V[I];
			}
			JCF = append(reverse(JBK[1]),JCF);
		} else {
			AlgV = algptorat(Alpha[K]);
			AlphaV = vector(DK);
			Sym = 1;
			Block = JBK[1];
			for ( U = 0; U < DK; U++ ) {
				AlphaV[U] = strtov(asciitostr([Prefix])+rtostr(U));
				Sym *= (x-AlphaV[U]);
				for ( T = JBK[0]; T != []; T = cdr(T), J++ ) {
					V = car(T);
					for ( I = 0; I < N; I++ )
						Q[I][J] = subst(algptorat(V[I]),AlgV,AlphaV[U]);
				}
			}
			for ( U = 0; U < DK; U++ )
				for ( Tmp = Block; Tmp != []; Tmp = cdr(Tmp) )
					JCF = cons([AlphaV[U],car(Tmp)[1],car(Tmp)[2]],JCF);
			F = Sym-FM[K][0];
			for ( I = 0, T = []; I < DK; I++ )
				T = cons(coef(F,I,x),T);
			DefIdeal = cons(T,DefIdeal);
			Prefix++;
		}
	}
	return [Q,reverse(JCF),reverse(DefIdeal)];
}

def is_integer(A)
{
	if ( !A || (type(A)==1 && ntype(A)==0 && dn(A)==1 ) ) return 1;
	else return 0;
}

def reduce_vect(M,Base)
{
	N = size(M)[0];
	M = dup_vect(M);
	for ( T = Base; T != []; T = cdr(T) ) {	
		B = car(T);
		BM = B[0]; I = B[1];
		if ( M[I] ) {
			/* M <- M*BM[I]-M[I]*BM */
			C1 = BM[I]; C2 = M[I];
			if ( is_integer(C1) && is_integer(C2) ) {
				G = igcd(C1,C2);
				C1 = idiv(C1,G);
				C2 = idiv(C2,G);
			}
			mulsub_vect(M,C1,BM,C2);
		}
	}
	M = map(simpalg,M);
	for ( I = 0; I < N; I++ )
		if ( M[I] ) return [M,I];
	return [0,-1];
}

def print_kernel(Ker)
{
	print("{",0);
	for ( T = Ker; T != []; T = cdr(T) ) {
		print_vect(car(T));
		if ( cdr(T) != [] ) print(",",0);
	}
	print("}");
}

def print_vect(V)
{
	K = V[0]; S = V[1];
	N = size(K)[0];
	if ( S < 0 )
		for ( S = 0; S < N; S++ )
			if ( K[S] ) break;
	if ( S == N )
		print(0,0);
	else {
		print("("+rtostr(K[S])+")e_{"+rtostr(S+1)+"}",0);	
		for ( I = S+1; I < N; I++ )
			if ( K[I] ) print("+("+rtostr(K[I])+")e_{"+rtostr(I+1)+"}",0);	
	}
}

def compute_jordan_block_basis(A,Alpha,E,Verbose)
{
	NI = unit_mat(size(A)[0]);
	N = A-Alpha*NI;
	if ( Verbose ) Ns = rtostr("A-(")+rtostr(Alpha)+")E";
	Ker = vector(E+1);
	for ( I = 0; I <= E; I++ ) {
		if ( Verbose ) print("Ker("+Ns+")^"+rtostr(I)+"=",0);
		Ker[I] = compute_kernel(NI);
		if ( I < E ) NI *= N;
		if ( Verbose )
			print_kernel(Ker[I]);
	}
	/* Block = a list of block info : [Alpha,K,M] -> J(Alpha,K) x M */
	/* PW : previous W */
	for ( PW = AllBase = Block = [], I = E-1; I >= 0; I-- ) {
		/* NPW = N*PW */
		NPW = [];
		/* create a basis of Ker[I]+N*PW */
		if ( Verbose ) {
			print("Ker("+Ns+")^"+rtostr(I)+"+NW_"+rtostr(I+2)+"=",0);
		}
		for ( T = PW, Base = Ker[I]; T != []; T = cdr(T) ) {
			NPW = cons(N*car(T),NPW);
			Red = reduce_vect(car(NPW),Base);
			Base = insert_to_base([ptozp_vect(Red[0])[0],Red[1]],Base);
		}
		/* number of elements to be added */
		DW = length(Ker[I+1])-length(Ker[I])-length(NPW);
		if ( Verbose ) {
			print_kernel(Base);
			print("number of bases to be added=",0); print(DW);
		}
		if ( DW ) Block = cons([Alpha,I+1,DW],Block);
		for ( PW = NPW, T = Ker[I+1]; DW > 0; T = cdr(T) ) {
			/* take an element of Ker[I+1] and reduce it by Base */
			Red = reduce_vect(car(T)[0],Base);
			if ( Red[0] ) {
				if ( Verbose ) {
					print("a new basis element=",0); print_vect(car(T)); print("");
				}
				/* we found a new basis element of height I+1 */
				DW--;
				/* append the new element to PW */
				PW = cons(Red[0],PW);
				/* append the new element to Base */
				Base = insert_to_base(Red,Base);
				/* Kry = [N^I*v,N^(I-1)*v,...,v] */
				for ( Kry = [car(T)[0]], K = 0; K < I; K++ )
					Kry = cons(map(simpalg,N*car(Kry)),Kry);
				if ( Verbose ) {
					print("added basis={",0);
					for ( Tmp = Kry; Tmp != []; Tmp = cdr(Tmp) ) {
						print_vect([car(Tmp),-1]);
						if ( cdr(Tmp) != [] ) print(",",0);
					}
					print("}");
				}
				AllBase = append(reverse(Kry),AllBase);
			}
		}
	}
	return [reverse(AllBase),reverse(Block)];
}

def next_combination(N,K,C)
{
	if ( !K ) return 0;
	for ( I = K-1; I >= 0 && C[I] == N+I-K; I-- );
	if ( I < 0 ) return 0;
	T = ++C[I++];
	for ( T++ ; I < K; I++, T++ ) C[I] = T;
	return 1;
}

def reduce_symmetric_form(F,E,V)
{
	N = size(E)[0]-1;
	F = dp_ptod(F,V);
	W = vector(N+1);
	for ( I = 0; I <= N; I++ )
		W[I] = strtov("e"+rtostr(I));
	R = 0;
	while ( F ) {
		HT = dp_ht(F);
		S = E[0];
		SE = 1;
		while ( dp_td(HT) ) {
			for ( I = N; I >= 1; I-- )
				if ( dp_redble(HT,dp_ht(E[I])) ) {
					HT = dp_subd(HT,dp_ht(E[I]));
					S *= E[I];
					SE *= W[I];
					break;
				}
		}
		C = dp_hc(F);
		F -= S*C;
		R += SE*C;
	}
	return R;
}

def generate_elementary_symmetric_form(N)
{
	C = vector(N);
	V = getopt(vars);
	if ( type(V) == -1 ) {
		V = vector(N);
		for ( I = 0; I < N; I++ )
			V[I] = strtov("x"+rtostr(I));
	}
	R = [1];
	for ( I = 1; I <= N; I++ ) {
		for ( J = 0; J < I; J++ ) C[J] = J;
		EI = 0;
		do {
			for ( J = 0, E = 1; J < I; J++ )
				E *= V[C[J]];
			EI += E;
		} while ( next_combination(N,I,C) );
		R = cons(EI,R);
	}
	dp_ord(2);
	V = vtol(V);
	return [ltov(map(dp_ptod,reverse(R),V)),V];
}

def partial_fraction(F,E)
{
	V = var(F);
	N = deg(F,V);
	FE = F^E;
	R = -1;
	L = generate_elementary_symmetric_form(N);
	ES = L[0];
	ESV	=L[1];
	NE = N*E;
	for ( I = 0; I < N; I++ ) {
		for ( J = 0, H = 0; J < E; J++ ) {
			T = strtov("h"+rtostr(I)+"_"+rtostr(J));
			H += T*V^J;
		}
		Q = sdiv(FE,(V-ESV[I])^E);
		R += srem(Q*H,subst(F,V,ESV[I]));
	}
	HV = [];
	for ( J = E-1; J >= 0; J-- )
		for ( I = N-1; I >= 0; I-- )
			HV = cons(strtov("h"+rtostr(I)+"_"+rtostr(J)),HV);
	for ( I = 0, L2 = []; I < NE; I++ )
		L2 = cons(coef(R,I,V),L2);
	for ( I = 1, SGN = -1, L1 = []; I <= N; I++, SGN =-SGN )
		L1 = cons(dp_dtop(ES[I],ESV)-SGN*coef(F,I,V),L1);
	L = append(L1,L2);
#if 0
	G2 = dp_gr_main(L2,HV,0,0,2);
	G1 = dp_gr_main(L1,ESV,0,0,2);
	AllV = append(HV,ESV);
	G = nd_gr(append(G1,G2),AllV,0,2);
#else
	AllV = append(HV,ESV);
	G = nd_gr(L,AllV,0,2);
#endif
	R = 0;
	for ( I = 0; I < N; I++ ) {
		for ( J = 0, H = 0; J < E; J++ ) {
			T = strtov("h"+rtostr(I)+"_"+rtostr(J));
			H += T*V^J;
		}
		Q = sdiv(FE,(V-ESV[I])^E);
		T = p_true_nf(srem(ESV[I]*Q*H,subst(F,V,ESV[I])),G,AllV,2);
		R += T[0]/T[1];
	}
	return R;
}
endmodule$
end$