load("solve")$
load("gr")$
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 wsort(A,B,C){
D=newvect(length(B))$
for(I=0;I<length(B);I++)
D[I]=[A[I],B[I],C[I]]$
D=qsort(D,worder)$
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]$
}
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(C[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 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$
if(J<size(ResVars)[0])
ResVars[J]=Res[I][1]$
}
return(ResVars)$
}
def makeret(Res,Vars){
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(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(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){
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=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"))$
return [wsort(ResVars,Vars,ResVars)]$
}
else
return []$
}
else
return []$
}
def leastsq(NormMat,ExpMat,Vars){
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)$
if(TMP1[0]==0){
TMP=roundret(TMP1[1]*1.0)$
if(TMP!=[])
RET=cons(wsort(TMP1[1],Vars,TMP),RET)$
RET=cons(wsort(TMP1[1],Vars,
map(drint,TMP1[1]*1.0)),RET)$
return RET$
}
else{
RET=cons(wsort(TMP1[1],Vars,TMP1[1]*1.0),RET)$
return RET$
}
}
else
return RET$
}
def weightr(ExpMat,Vars,PolyListNum,OneMat){
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)$
Rea=cons(ExtVars[F],Rea)$
Res=getgcd(Res,Rea)$
if(nonposdegchk(Res)){
TMP1=makeret(Res,Vars)$
if(TMP1[0]==0){
TMP=roundret(TMP1[1]*1.0)$
if(TMP!=[])
RET=cons(wsort(TMP1[1],Vars,TMP),RET)$
RET=cons(wsort(TMP1[1],Vars,
map(drint,TMP1[1]*1.0)),RET)$
}
else{
RET=cons(wsort(TMP1[1],Vars,TMP1[1]*1.0),RET)$
}
}
}
return [NormMat,RET]$
}
def weight(PolyList,Vars,FLAG){
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=[]$
/* first */
TMP=qcheck(PolyList,Vars)$
if(TMP!=[]){
RET=append(RET,TMP)$
return cons(1,RET)$
}
dp_ord(2)$
PolyListNum=length(PolyList)$
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)$
/* second */
if(FLAG){
TMP=weightr(ExpMat,Vars,PolyListNum,OneMat)$
RET=append(RET,TMP[1])$
}
/* third */
if(FLAG)
RET=append(RET,leastsq(TMP[0],ExpMat,Vars))$
else
RET=append(RET,leastsq(0,ExpMat,Vars))$
/* forth */
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))$
}
RET=derase(RET)$
return cons(0,RET)$
}
end$