load("solve")$
load("gr")$

def nonposdegchk(Res){

	for(I=0;I<length(Res);I++)
		if(Res[I][1]<=0)
			return 0$

	return 1$
}

def resvars(Res,Vars){

	ResVars=newvect(length(Vars),Vars)$

	for(I=0;I<length(Res);I++){
	
		for(J=0;J<size(ResVars)[0];J++)
			if(Res[I][0]==ResVars[J])
				break$

		ResVars[J]=Res[I][1]$
	}

	return(ResVars)$
}

def makeret1(Res,Vars){

	VarsNum=length(Vars)$

	ResVec=newvect(VarsNum,Vars)$

	for(I=0,M=0;I<length(Res);I++){

		for(J=0;J<VarsNum;J++)
			if(Res[I][0]==Vars[J])
				break$

		if(J<VarsNum){
			ResVec[J]=Res[I][1]$

			if(type(ResVec[J])==1){
				if(M==0)
					M=ResVec[J]$
				else
					if(ResVec[J]<M)
						M=ResVec[J]$
			}
		}

	}

	for(F=0,I=0;I<VarsNum;I++)
		if(type(ResVec[I])!=1){
			F=1$
			break$
		}

	if(F==0)
		for(I=0;I<VarsNum;I++)
			ResVec[I]=ResVec[I]/M*1.0$

	for(I=0;I<VarsNum;I++)
		for(J=0;J<length(Vars);J++)
			ResVec[I]=subst(ResVec[I],Vars[J],
				strtov(rtostr(Vars[J])+"_deg"))$

	ResVec=cons(F,vtol(ResVec))$
	return ResVec$
}

def junban1(A,B){
	return (nmono(A)<nmono(B) ? -1:(nmono(A)>nmono(B) ? 1:0))$
}

def junban2(A,B){

	for(I=0;I<size(A)[0];I++){
		if(A[I]<B[I])
			return 1$
		
		if(A[I]>B[I])
			return -1$
	}

	return 0$
}

def roundret(V){

	VN=length(V)$
	RET0=newvect(VN,V)$

	for(I=1;I<1000;I++){
		RET1=I*RET0$
		for(J=0;J<VN;J++){
			X=drint(RET1[J])$
			if(dabs(X-RET1[J])<0.2)
				RET1[J]=X$
			else
				break$
		}
		if(J==VN)
			break$
	}
	
	if(I==1000)
		return []$
	else
		return RET1$
}

def chkou(L,ExpMat,CHAGORD){

	P=1$
	F=ExpMat[L]$

	for(I=0;I<L;I++){
		Q=ExpMat[L][CHAGORD[I]]$
		for(J=0;J<size(ExpMat[0])[0];J++){
			ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
				*ExpMat[L][CHAGORD[J]]-
					Q*ExpMat[I][CHAGORD[J]])/P)$
		}

		P=ExpMat[I][CHAGORD[I]]$
	}

	for(J=0;J<size(ExpMat[0])[0];J++)
		if(ExpMat[L][CHAGORD[J]]!=0)
			break$

	if(J==size(ExpMat[0])[0])
		return L$
	else{
		TMP=CHAGORD[L]$
		CHAGORD[L]=CHAGORD[J]$
		CHAGORD[J]=TMP$
		return (L+1)$
	}
}

def qcheck0(PolyList,Vars){

	RET=[]$
	PolyListNum=length(PolyList)$
	VarsNum=length(Vars)$

	ExpMat=newvect(VarsNum)$
	CHAGORD=newvect(VarsNum)$
	for(I=0;I<VarsNum;I++)
		CHAGORD[I]=I$

	L=0$
	for(I=0;I<PolyListNum;I++){
		Poly=dp_ptod(PolyList[I],Vars)$
		BASE0=dp_etov(dp_ht(Poly))$
		Poly=dp_rest(Poly)$
		for(;Poly!=0;Poly=dp_rest(Poly)){
			ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
			L=chkou(L,ExpMat,CHAGORD)$
			if(L==VarsNum-1){
				RET=cons(ExpMat,RET)$
				RET=cons(CHAGORD,RET)$
				RET=cons(L,RET)$
				return RET$
			}
		}	
	}
	
	RET=cons(ExpMat,RET)$
	RET=cons(CHAGORD,RET)$
	RET=cons(L,RET)$
	return RET$
}

def inner(A,B){

	SUM=0$
	for(I=0;I<size(A)[0];I++)
		SUM+=A[I]*B[I]$

	return SUM$
}

def checktd(PolyList,Vars,ResVars){

	PolyListNum=length(PolyList)$
	VarsNum=length(Vars)$

	L=0$
	for(I=0;I<PolyListNum;I++){
		Poly=dp_ptod(PolyList[I],Vars)$
		J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
		Poly=dp_rest(Poly)$
		for(;Poly!=0;Poly=dp_rest(Poly))
			if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
				return 0$
	}
	
	return 1$
}

def getgcd(A,B){

	VarsNumA=length(A)$
	VarsNumB=length(B)$

	C=newvect(VarsNumB,B)$

	for(I=0;I<VarsNumA;I++){

		for(J=0;J<VarsNumB;J++)
			if(C[J]==A[I][0])
				break$

		C[J]=A[I][1]$
	}

	D=0$
	for(I=0;I<VarsNumB;I++)
		D=gcd(D,C[I])$

	if(D!=0){

		for(I=0;I<VarsNumB;I++)
			C[I]=red(C[I]/D)$

	}

	for(L=1,D=0,I=0;I<VarsNumB;I++){

		if(type(C[I])==1){
			L=ilcm(L,dn(C[I]))$
			D=igcd(D,nm(C[I]))$
		}
		else
			break$

	}

	if(I==VarsNumB)
		for(I=0;I<VarsNumB;I++)
			C[I]=C[I]*L/D$

	RET=newvect(VarsNumB)$
	for(I=0;I<VarsNumB;I++){
		RET[I]=newvect(2)$
		RET[I][0]=B[I]$
		RET[I][1]=C[I]$
	}

	return vtol(map(vtol,RET))$
}

def qcheck(PolyList,Vars){

	RET=[]$
	Res=qcheck0(PolyList,Vars)$
	VarsNum=length(Vars)$

	IndNum=Res[0]$
	CHAGORD=Res[1]$
	ExpMat=Res[2]$

	SolveList=[]$
	for(I=0;I<IndNum;I++){
		TMP=0$
		for(J=0;J<VarsNum;J++)
			TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$

		SolveList=cons(TMP,SolveList)$
	}

	VarsList=[]$
	for(I=0;I<VarsNum;I++)
		VarsList=cons(Vars[CHAGORD[I]],VarsList)$

	Rea=vars(SolveList)$
	Res=solve(reverse(SolveList),reverse(VarsList))$

	if(nonposdegchk(Res)){

		Res=getgcd(Res,Rea)$
   	ResVars=resvars(Res,Vars)$

		if(checktd(PolyList,Vars,ResVars)==1){

			for(J=0;J<length(Vars);J++)
				ResVars=map(subst,ResVars,Vars[J],
					strtov(rtostr(Vars[J])+"_deg"))$

			RET=cons([vtol(ResVars),ResVars,[]],RET)$
			return cons(1,RET)$
		}
		else
			return cons(0,RET)$
	}
	else
		return cons(0,RET)$

}

def weight(PolyList,Vars){

	Vars0=vars(PolyList)$
	Vars1=[]$
	for(I=0;I<length(Vars);I++)
		if(member(Vars[I],Vars0))
			Vars1=cons(Vars[I],Vars1)$

	Vars=reverse(Vars1)$

	RET=[]$

	TMP=qcheck(PolyList,Vars)$

	if(car(TMP)==1){
		RET=cdr(TMP)$
		RET=cons(Vars,RET)$
		RET=cons(1,RET)$
		return RET$	
	}

	dp_ord(2)$

	PolyListNum=length(PolyList)$
	VPolyList=qsort(newvect(PolyListNum,PolyList),junban1)$
	VPolyList=vtol(VPolyList)$

	ExpMat=[]$
	for(I=0;I<PolyListNum;I++)
		for(Poly=dp_ptod(VPolyList[I],Vars);Poly!=0;Poly=dp_rest(Poly))
			ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$

	ExpMat=reverse(ExpMat)$
	ExpMat=newvect(length(ExpMat),ExpMat)$


/* first */

	ExpMatRowNum=size(ExpMat)[0]$
	ExpMatColNum=size(ExpMat[0])[0]$
	ExtMatColNum=ExpMatColNum+PolyListNum$

	OneMat=newvect(PolyListNum+1,[0])$
	for(I=0,SUM=0;I<PolyListNum;I++){
		SUM+=nmono(VPolyList[I])$
		OneMat[I+1]=SUM$
	}

	RevOneMat=newvect(ExpMatRowNum)$
	for(I=0;I<PolyListNum;I++)
		for(J=OneMat[I];J<OneMat[I+1];J++)
			RevOneMat[J]=I$

	NormMat=newmat(ExpMatColNum,ExtMatColNum)$

	for(I=0;I<ExpMatColNum;I++)
		for(J=0;J<ExpMatColNum;J++)
			for(K=0;K<ExpMatRowNum;K++)
				NormMat[I][J]+=ExpMat[K][I]*ExpMat[K][J]$

	for(I=0;I<ExpMatColNum;I++)
		for(J=0;J<PolyListNum-1;J++)
			for(K=OneMat[J];K<OneMat[J+1];K++)
				NormMat[I][J+ExpMatColNum]-=ExpMat[K][I]$

	for(I=0;I<ExpMatColNum;I++)
		for(J=OneMat[PolyListNum-1];J<OneMat[PolyListNum];J++)
			NormMat[I][ExtMatColNum-1]+=ExpMat[J][I]$

	NormMat2=newmat(PolyListNum-1,ExpMatColNum+1)$

	for(I=0;I<PolyListNum-1;I++)
		for(J=0;J<ExpMatColNum;J++)
			for(K=OneMat[I];K<OneMat[I+1];K++)
				NormMat2[I][J]-=ExpMat[K][J]$

	for(I=0;I<PolyListNum-1;I++)
		NormMat2[I][ExpMatColNum]=OneMat[I+1]-OneMat[I]$

	ExtVars=Vars$
	for(I=0;I<PolyListNum-1;I++)
		ExtVars=append(ExtVars,[uc()])$

	SolveList=[]$
	for(I=0;I<ExpMatColNum;I++){
		TMP=0$
		for(J=0;J<ExtMatColNum-1;J++)
			TMP+=NormMat[I][J]*ExtVars[J]$

		TMP-=NormMat[I][ExtMatColNum-1]$
		SolveList=cons(TMP,SolveList)$
	}

	for(I=0;I<PolyListNum-1;I++){
		TMP=0$
		for(J=0;J<ExpMatColNum;J++)
			TMP+=NormMat2[I][J]*ExtVars[J]$

		TMP+=NormMat2[I][ExpMatColNum]*ExtVars[I+ExpMatColNum]$

		SolveList=cons(TMP,SolveList)$
	}

	Rea=vars(SolveList)$
	Res=solve(SolveList,reverse(ExtVars))$

	if(nonposdegchk(Res)){
		Res=getgcd(Res,Rea)$
		TMP1=makeret1(Res,Vars);
		if(car(TMP1)==0){	
			TMP2=roundret(cdr(TMP1));
			TMP3=map(drint,cdr(TMP1))$
			RET=cons([cdr(TMP1),newvect(length(TMP3),TMP3),TMP2],RET)$
		}
		else
			RET=cons([cdr(TMP1),[],[]],RET)$
	}

/* second */

	NormMat=newmat(ExpMatColNum,ExpMatColNum+1)$

	for(I=0;I<ExpMatColNum;I++)
		for(J=0;J<ExpMatColNum;J++)
			for(K=0;K<ExpMatRowNum;K++)
				NormMat[I][J]+=ExpMat[K][I]*ExpMat[K][J]$

	for(I=0;I<ExpMatColNum;I++)
		for(J=0;J<ExpMatRowNum;J++)
			NormMat[I][ExpMatColNum]+=ExpMat[J][I]$

	SolveList=[]$
	for(I=0;I<ExpMatColNum;I++){
		TMP=0$
		for(J=0;J<ExpMatColNum;J++)
			TMP+=NormMat[I][J]*Vars[J]$

		TMP-=NormMat[I][ExpMatColNum]$
		SolveList=cons(TMP,SolveList)$
	}

	Rea=vars(SolveList)$
	Res=solve(SolveList,Vars)$

	if(nonposdegchk(Res)){
		Res=getgcd(Res,Rea)$
		TMP1=makeret1(Res,Vars);
		if(car(TMP1)==0){	
			TMP2=roundret(cdr(TMP1));
			TMP3=map(drint,cdr(TMP1))$
			RET=cons([cdr(TMP1),newvect(length(TMP3),TMP3),TMP2],RET)$
		}
		else
			RET=cons([cdr(TMP1),[],[]],RET)$
	}

/* third */

	ExpMat=qsort(ExpMat,junban2)$
	ExpMat2=[]$
	for(I=0;I<size(ExpMat)[0];I++)
		if(car(ExpMat2)!=ExpMat[I])
			ExpMat2=cons(ExpMat[I],ExpMat2)$

	ExpMat=newvect(length(ExpMat2),ExpMat2)$
	ExpMatRowNum=size(ExpMat)[0]$
	ExpMatColNum=size(ExpMat[0])[0]$

	NormMat=newmat(ExpMatColNum,ExpMatColNum+1)$

	for(I=0;I<ExpMatColNum;I++)
		for(J=0;J<ExpMatColNum;J++)
			for(K=0;K<ExpMatRowNum;K++)
				NormMat[I][J]+=ExpMat[K][I]*ExpMat[K][J]$

	for(I=0;I<ExpMatColNum;I++)
		for(J=0;J<ExpMatRowNum;J++)
			NormMat[I][ExpMatColNum]+=ExpMat[J][I]$

	SolveList=[]$
	for(I=0;I<ExpMatColNum;I++){
		TMP=0$
		for(J=0;J<ExpMatColNum;J++)
			TMP+=NormMat[I][J]*Vars[J]$

		TMP-=NormMat[I][ExpMatColNum]$
		SolveList=cons(TMP,SolveList)$
	}

	Rea=vars(SolveList)$
	Res=solve(SolveList,Vars)$

	if(nonposdegchk(Res)){
		Res=getgcd(Res,Rea)$
		TMP1=makeret1(Res,Vars);
		if(car(TMP1)==0){	
			TMP2=roundret(cdr(TMP1));
			TMP3=map(drint,cdr(TMP1))$
			RET=cons([cdr(TMP1),newvect(length(TMP3),TMP3),TMP2],RET)$
		}
		else
			RET=cons([cdr(TMP1),[],[]],RET)$
	}

	RET=cons(Vars,reverse(RET))$
	RET=cons(0,RET)$
	return RET$
}

def average(PolyList,Vars){

	RET=[]$
	dp_ord(2)$

	PolyListNum=length(PolyList)$

	ExpMat=[]$
	for(I=0;I<PolyListNum;I++)
		for(Poly=dp_ptod(PolyList[I],Vars);Poly!=0;Poly=dp_rest(Poly))
			ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$

	ExpMat=reverse(ExpMat)$
	ExpMat=newvect(length(ExpMat),ExpMat)$

	ExpMat=qsort(ExpMat,junban2)$
	ExpMat2=[]$
	for(I=0;I<size(ExpMat)[0];I++)
		if(car(ExpMat2)!=ExpMat[I])
			ExpMat2=cons(ExpMat[I],ExpMat2)$

	ExpMat=newvect(length(ExpMat2),ExpMat2)$
	ExpMatRowNum=size(ExpMat)[0]$
	ExpMatColNum=size(ExpMat[0])[0]$

	Res=newvect(ExpMatColNum);
	for(I=0;I<ExpMatColNum;I++)
		Res[I]=newvect(2,[Vars[I]])$

	for(I=0;I<ExpMatRowNum;I++)
		for(J=0;J<ExpMatColNum;J++)
			Res[J][1]+=ExpMat[I][J]$

	for(I=0;I<ExpMatColNum;I++)
		if(Res[I][1]==0)
			Res[I][1]=1$
		else
			Res[I][1]=1/Res[I][1]$

	RET=cons(makeret(vtol(Res),Vars,1),RET)$
	RET=cons(Vars,RET)$

	return RET$
}

end$