Return to weight CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / lib

Changed weight

load("solve")$

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 makeret(Res,Vars,B){

	VarsNum=length(Vars)$

	ResMat=newvect(VarsNum)$
	for(I=0;I<VarsNum;I++)
		ResMat[I]=newvect(2)$

	for(I=0;I<VarsNum;I++){
		ResMat[I][0]=Vars[I]$
		ResMat[I][1]=Vars[I]$
	}
	
	for(F=0,L=1,D=0,M=1,I=0;I<length(Res);I++){

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

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

			if(F==0 && type(K)==1){
				if(B==2){
					L=ilcm(L,dn(K))$
					D=igcd(D,nm(K))$
				}
				else{
					if(K<M)
						M=K$
				}
			}
			else
				F=1$

		}
	}

	if(F==0)
		if(B==2)
			for(I=0;I<VarsNum;I++)
				ResMat[I][1]=ResMat[I][1]*L/D$
		else
			for(I=0;I<VarsNum;I++)
				ResMat[I][1]=ResMat[I][1]/M*1.0$

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

	ResMat=map(vtol,ResMat)$
	return(vtol(ResMat))$

}

def afo(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 weight(PolyList,Vars){

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

	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(PolyList[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)$
	}

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

	RET=[]$
	if(nonposdegchk(Res)){

		ResVars=resvars(Res,ExtVars)$
		ResVars=append(vtol(ResVars),[1])$

		for(I=0;I<ExpMatRowNum;I++){
			TMP=0$
			for(J=0;J<ExpMatColNum;J++)
				TMP+=ExpMat[I][J]*ResVars[J]$

			TMP-=ResVars[RevOneMat[I]+ExpMatColNum]$

			if(TMP!=0)
				break$
		}

		if(I==ExpMatRowNum){
			print("complitely homogenized")$
			return([makeret(Res,Vars,2)])$
		}
		else
			RET=cons(makeret(Res,Vars,1),RET)$
	}

	ExpMat=qsort(ExpMat,afo)$
	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)$
	}

	Res=solve(SolveList,Vars)$

	if(nonposdegchk(Res))
		RET=cons(makeret(Res,Vars,1),RET)$

	return reverse(RET)$
}

end$