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

def nonzerovec(A){

	for(I=0;I<size(A)[0];I++)
		if(A[I]!=0)
			return 1$

	return 0$
}

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

def worder(A,B){
	return (A[0]<B[0] ? 1:(A[0]>B[0] ? -1:0))$
}

def bsort(A){

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

def perm(I,P,TMP){

	if(I>0){
		TMP=perm(I-1,P,TMP)$
		for(J=I-1;J>=0;J--){
			T=P[I]$
			P[I]=P[J]$
			P[J]=T$
			TMP=perm(I-1,P,TMP)$
			T=P[I]$
			P[I]=P[J]$
			P[J]=T$
		}

		return TMP$
	}
	else{
		for(TMP0=[],K=0;K<size(P)[0];K++)
			TMP0=cons(P[K],TMP0)$
				
		TMP=cons(TMP0,TMP)$
		return TMP$
	}
}

def marge(A,B){

	RET=[]$
	for(I=0;I<length(A);I++)
		for(J=0;J<length(B);J++)
			RET=cons(append(A[I],B[J]),RET)$

	return RET$
}

def wsort(A,B,C,FLAG){

	if(FLAG==0){
		D=newvect(length(B))$
		for(I=0;I<length(B);I++)
			D[I]=[A[I],B[I],C[I]]$

		D=bsort(D)$
		E=[]$
		for(I=0;I<length(B);I++)
			E=cons(D[I][1],E)$
		E=reverse(E)$
		F=[]$
		for(I=0;I<length(B);I++)
			F=cons(D[I][2],F)$
		F=reverse(F)$
	
		return [[E,F]]$
	}
	else{
		D=newvect(length(B))$
		for(I=0;I<length(B);I++)
			D[I]=[A[I],B[I],C[I]]$

		D=qsort(D,worder)$
		D0=[]$

		for(I=0,J=0,TMP=[],X=0;I<size(D)[0];I++){
			if(X==D[I][0])
				TMP=cons(cdr(D[I]),TMP)$
			else{
				D0=cons(TMP,D0)$	
				TMP=[]$
				TMP=cons(cdr(D[I]),TMP)$
				X=car(D[I])$
			}
		}
		D0=cdr(reverse(cons(TMP,D0)))$
		D0=map(ltov,D0)$
		for(I=0,TMP=[[]];I<length(D0);I++){
			TMP0=perm(length(D0[I])-1,D0[I],[])$
			TMP=marge(TMP,TMP0)$
		}
		
		RET=[]$
		for(I=0;I<length(TMP);I++){
			TMP0=[]$
			TMP1=[]$
			for(J=0;J<length(TMP[I]);J++){
				TMP0=cons(TMP[I][J][0],TMP0)$
				TMP1=cons(TMP[I][J][1],TMP1)$
			}	
			TMP0=reverse(TMP0)$
			TMP1=reverse(TMP1)$

			RET=cons([TMP0,TMP1],RET)$
		}
		
		return RET$
	}	
}

def derase(A){

	B=newvect(length(A),A)$
	B=qsort(B,junban)$
	C=[]$
	for(I=0;I<size(B)[0];I++)
		if(car(C)!=B[I])
			C=cons(B[I],C)$

	return reverse(C)$
}

def nonposdegchk(Res){

	for(I=0;I<length(Res);I++)
		if(Res[I][1]<=0)
			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(B[J]==A[I][0])
				break$

		if(J<VarsNumB)
			C[J]=A[I][1]$
	}

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

	if(D!=0){
		C=C/D$
		C=map(red,C)$
	}

	for(L=1,D=0,I=0;I<VarsNumB;I++){
		if(type(TMP=dn(C[I]))==1)
			L=ilcm(L,TMP)$

		if(type(TMP=nm(C[I]))==1)
			D=igcd(D,TMP)$
	}

	C=C*L$
	if(D!=0)
		C=C/D$

	RET=[]$
	for(I=0;I<VarsNumB;I++)
		RET=cons([B[I],C[I]],RET)$

	return RET$
}

def makeret(Res,Vars,FLAG){

	ResNum=length(Res)$
	VarsNum=length(Vars)$

	ResVec=newvect(ResNum)$

        for(M=0,I=0;I<ResNum;I++){
                if(member(Res[I][0],Vars)){
                        ResVec[I]=Res[I][1]$

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

	if(M!=0)
		ResVec=ResVec/M;

 	RET=newvect(VarsNum,Vars)$

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

		if(J<VarsNum)
			RET[J]=ResVec[I]$
	}
			

	for(J=0;J<length(Vars);J++)
		RET=map(subst,RET,Vars[J],
			strtov(rtostr(Vars[J])+"_deg"))$
			
	for(I=0;I<VarsNum;I++)
		if(type(RET[I])!=1)
			return [1,RET]$

	return [0,RET]$
}

def roundret(V){

	VN=size(V)[0]$

	RET0=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){

	for(P=1,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 qcheckmain(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)
				return [L,CHAGORD,ExpMat]$
		}	
	}
	
	return [L,CHAGORD,ExpMat]$
}

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 qcheck(PolyList,Vars,FLAG){

	RET=[]$
	Res=qcheckmain(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)$
	}

	Rea=vars(SolveList)$

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

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

	if(nonposdegchk(Res)){

	   	ResVars=makeret(Res,Vars,0)$

		if(checktd(PolyList,Vars,ResVars[1])==1){
			if(ResVars[0]==0){
				RET=append(RET,wsort(ResVars[1],Vars,
					ResVars[1],FLAG))$
				return RET$
			}
			else{
				RET=append(RET,[[Vars,vtol(ResVars[1])]])$
				return RET$
			}
		}
		else
			return []$
	}
	else
		return []$

}

def leastsq(NormMat,ExpMat,Vars,FLAG){

	RET=[]$

	ExpMatRowNum=size(ExpMat)[0]$
	ExpMatColNum=size(ExpMat[0])[0]$

	if(NormMat==0){
		NormMat=newmat(ExpMatColNum,ExpMatColNum)$

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

	BVec=newvect(ExpMatColNum)$

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

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

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

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

	Rea=vars(SolveList)$

	VarsList=[]$
	for(I=0;I<length(Vars);I++)
		if(member(Vars[I],Rea))
			VarsList=cons(Vars[I],VarsList)$

	Res=solve(SolveList,VarsList)$
	Res=getgcd(Res,Rea)$

	if(nonposdegchk(Res)){

		TMP1=makeret(Res,Vars,1)$

		if(TMP1[0]==0){	
			TMP=roundret(TMP1[1])$
			if(TMP!=[])
				RET=append(RET,wsort(TMP1[1],Vars,TMP,FLAG))$

			RET=append(RET,wsort(TMP1[1],Vars,
				map(drint,TMP1[1]*1.0),FLAG))$

			return RET$
		}
		else{
			RET=append(RET,[[Vars,vtol(TMP1[1]*1.0)]])$
			return RET$
		}
	}
	else
		return RET$

}

def weightr(ExpMat,Vars,PolyListNum,OneMat,FLAG){

	RET=[]$

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

	ExtVars=reverse(Vars)$
	for(I=0;I<PolyListNum;I++)
		ExtVars=cons(uc(),ExtVars)$

	ExtVars=reverse(ExtVars)$

	NormMat=newmat(ExpMatColNum,ExtMatColNum)$

	for(I=0;I<ExpMatColNum;I++)
		for(J=I;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;J++)
			for(K=OneMat[J];K<OneMat[J+1];K++)
				NormMat[I][J+ExpMatColNum]-=
					ExpMat[K][I]$

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

	for(F=0;F<ExtMatColNum;F++){
		SolveList=[]$
		for(I=0;I<ExpMatColNum;I++){
			if (F==I) 
				continue$

			TMP=0$
		
			for(J=0;J<I;J++)
				if(J!=F)
					TMP+=NormMat[J][I]*ExtVars[J]$

			for(J=I;J<ExtMatColNum;J++)
				if(J!=F)	
					TMP+=NormMat[I][J]*ExtVars[J]$
			
			if(F<I)
				TMP+=NormMat[F][I]$
			else
				TMP+=NormMat[I][F]$
	
			SolveList=cons(TMP,SolveList)$
		}		
		
		for(I=0;I<PolyListNum;I++){
			if(F==(I+ExpMatColNum))
				continue$
			
			TMP=0$
			for(J=0;J<ExpMatColNum;J++)
				if(J!=F)
					TMP+=NormMat[J][I+ExpMatColNum]
						*ExtVars[J]$

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

			if(F<ExpMatColNum)
				TMP+=NormMat[F][I+ExpMatColNum]$

			SolveList=cons(TMP,SolveList)$
		}	

		Rea=vars(SolveList)$

		SolVars=[]$
		for(I=0;I<ExtMatColNum;I++)
			if(I!=F && member(ExtVars[I],Rea))
				SolVars=cons(ExtVars[I],SolVars)$

		Res=solve(SolveList,SolVars)$
		Res=cons([ExtVars[F],1],Res)$

		TMP=[]$
		for(I=0;I<length(Rea);I++)
			if(member(Rea[I],Vars))
				TMP=cons(Rea[I],TMP)$
	
		if(member(ExtVars[F],Vars))
			TMP=cons(ExtVars[F],TMP)$

		Res=getgcd(Res,TMP)$
	
		if(nonposdegchk(Res)){

			TMP1=makeret(Res,Vars,1)$

			if(TMP1[0]==0){	
				TMP=roundret(TMP1[1])$
				if(TMP!=[])
					RET=append(RET,wsort(TMP1[1],Vars,
						TMP,FLAG))$

				RET=append(RET,wsort(TMP1[1],Vars,
					map(drint,TMP1[1]*1.0),FLAG))$
			}
			else{
				RET=append(RET,[[Vars,vtol(TMP1[1]*1.0)]])$
			}
		}

	}

	return [NormMat,RET]$
}

def weight(PolyList,Vars,FLAG1,FLAG2){

	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,FLAG2)$

	if(TMP!=[]){
		RET=append(RET,TMP)$
		return cons(1,RET)$
	}

	dp_ord(2)$

	PolyListNum=length(PolyList)$

	if(FLAG1){

		OneMat=newvect(PolyListNum+1,[0])$
		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)$
			}
			OneMat[I+1]=length(ExpMat)$
		}

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

		TMP=weightr(ExpMat,Vars,PolyListNum,OneMat,FLAG2)$
		RET=append(RET,TMP[1])$
		RET=append(RET,leastsq(TMP[0],ExpMat,Vars,FLAG2))$
	}
	else{
		ExpMat=[]$
		for(I=0;I<PolyListNum;I++){
			for(Poly=dp_ptod(PolyList[I],Vars);
				Poly!=0;Poly=dp_rest(Poly)){
				if(nonzerovec(TMP=dp_etov(dp_ht(Poly))))
					ExpMat=cons(TMP,ExpMat)$
			}
		}

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

		RET=append(RET,leastsq(0,ExpMat,Vars,FLAG2))$
	}

	ExpMat=qsort(ExpMat,junban)$

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

	if(size(ExpMat)[0]!=length(ExpMat2)){
		ExpMat=newvect(length(ExpMat2),ExpMat2)$
		RET=append(RET,leastsq(0,ExpMat,Vars,FLAG2))$
	}

	RET=derase(RET)$
	return cons(0,RET)$
}

end$