| version 1.9, 2003/11/21 08:07:16 |
version 1.32, 2004/02/14 18:28:39 |
|
|
| load("solve")$ |
load("solve")$ |
| load("gr")$ |
load("gr")$ |
| |
|
| def nonzerovec(A){ |
#define EPS 1E-6 |
| |
#define TINY 1E-20 |
| |
#define MAX_ITER 100 |
| |
|
| for(I=0;I<size(A)[0];I++) |
def rotate(A,I,J,K,L,C,S){ |
| if(A[I]!=0) |
|
| return 1$ |
|
| |
|
| return 0$ |
X=A[I][J]; |
| |
Y=A[K][L]; |
| |
A[I][J]=X*C-Y*S; |
| |
A[K][L]=X*S+Y*C; |
| |
|
| |
return 1; |
| } |
} |
| |
|
| |
def jacobi(N,A,W){ |
| |
|
| |
S=OFFDIAG=0.0; |
| |
|
| |
for(J=0;J<N;J++){ |
| |
|
| |
for(K=0;K<N;K++) |
| |
W[J][K]=0.0; |
| |
|
| |
W[J][J]=1.0; |
| |
S+=A[J][J]*A[J][J]; |
| |
|
| |
for(K=J+1;K<N;K++) |
| |
OFFDIAG+=A[J][K]*A[J][K]; |
| |
} |
| |
|
| |
TOLERANCE=EPS*EPS*(S/2+OFFDIAG); |
| |
|
| |
for(ITER=1;ITER<=MAX_ITER;ITER++){ |
| |
|
| |
OFFDIAG=0.0; |
| |
for(J=0;J<N-1;J++) |
| |
for(K=J+1;K<N;K++) |
| |
OFFDIAG+=A[J][K]*A[J][K]; |
| |
|
| |
if(OFFDIAG < TOLERANCE) |
| |
break; |
| |
|
| |
for(J=0;J<N-1;J++){ |
| |
for(K=J+1;K<N;K++){ |
| |
|
| |
if(dabs(A[J][K])<TINY) |
| |
continue; |
| |
|
| |
T=(A[K][K]-A[J][J])/(2.0*A[J][K]); |
| |
|
| |
if(T>=0.0) |
| |
T=1.0/(T+dsqrt(T*T+1)); |
| |
else |
| |
T=1.0/(T-dsqrt(T*T+1)); |
| |
|
| |
C=1.0/dsqrt(T*T+1); |
| |
|
| |
S=T*C; |
| |
|
| |
T*=A[J][K]; |
| |
|
| |
A[J][J]-=T; |
| |
A[K][K]+=T; |
| |
A[J][K]=0.0; |
| |
|
| |
for(I=0;I<J;I++) |
| |
rotate(A,I,J,I,K,C,S); |
| |
|
| |
for(I=J+1;I<K;I++) |
| |
rotate(A,J,I,I,K,C,S); |
| |
|
| |
for(I=K+1;I<N;I++) |
| |
rotate(A,J,I,K,I,C,S); |
| |
|
| |
for(I=0;I<N;I++) |
| |
rotate(W,J,I,K,I,C,S); |
| |
|
| |
} |
| |
} |
| |
} |
| |
|
| |
if (ITER > MAX_ITER) |
| |
return 0; |
| |
|
| |
for(I=0;I<N-1;I++){ |
| |
|
| |
K=I; |
| |
|
| |
T=A[K][K]; |
| |
|
| |
for(J=I+1;J<N;J++) |
| |
if(A[J][J]>T){ |
| |
K=J; |
| |
T=A[K][K]; |
| |
} |
| |
|
| |
A[K][K]=A[I][I]; |
| |
|
| |
A[I][I]=T; |
| |
|
| |
V=W[K]; |
| |
|
| |
W[K]=W[I]; |
| |
|
| |
W[I]=V; |
| |
} |
| |
|
| |
return 1; |
| |
} |
| |
|
| |
def interval2value(A,Vars){ |
| |
|
| |
B=atl(A)$ |
| |
|
| |
if(length(B)>2){ |
| |
print("bug")$ |
| |
return []$ |
| |
} |
| |
else if(length(B)==0){ |
| |
if(fop(A)==0) |
| |
return [Vars,1]$ |
| |
else |
| |
return []$ |
| |
} |
| |
else if(length(B)==1){ |
| |
|
| |
C=fargs(B[0])$ |
| |
D=vars(C)$ |
| |
E=solve(C,D)$ |
| |
|
| |
if(fop(B[0])==15) |
| |
return [Vars,E[0][1]+1]$ |
| |
else if(fop(B[0])==11) |
| |
return [Vars,E[0][1]-1]$ |
| |
else if(fop(B[0])==8) |
| |
return [Vars,E[0][1]]$ |
| |
else |
| |
return []$ |
| |
} |
| |
else{ |
| |
|
| |
C=fargs(B[0])$ |
| |
D=vars(C)$ |
| |
E=solve(C,D)$ |
| |
|
| |
C=fargs(B[1])$ |
| |
D=vars(C)$ |
| |
F=solve(C,D)$ |
| |
|
| |
return [Vars,(E[0][1]+F[0][1])/2]$ |
| |
} |
| |
|
| |
} |
| |
|
| |
def fixpointmain(F,Vars){ |
| |
|
| |
RET=[]$ |
| |
for(I=length(Vars)-1;I>=1;I--){ |
| |
|
| |
for(H=[],J=0;J<I;J++) |
| |
H=cons(Vars[J],H)$ |
| |
|
| |
G=interval2value(qe(ex(H,F)),Vars[I])$ |
| |
|
| |
if(G==[]) |
| |
return RET$ |
| |
else |
| |
RET=cons(G,RET)$ |
| |
|
| |
F=subf(F,G[0],G[1])$ |
| |
} |
| |
|
| |
G=interval2value(simpl(F),Vars[0])$ |
| |
|
| |
if(G==[]) |
| |
return RET$ |
| |
else |
| |
RET=cons(G,RET)$ |
| |
|
| |
return RET$ |
| |
} |
| |
|
| |
|
| |
def fixedpoint(A,FLAG){ |
| |
|
| |
Vars=vars(A)$ |
| |
|
| |
N=length(A)$ |
| |
|
| |
if (FLAG==0) |
| |
for(F=@true,I=0;I < N; I++ ) { F = F @&& A[I] @> 0$ } |
| |
else if (FLAG==1) |
| |
for(F=@true,I=0;I < N; I++ ) { F = F @&& A[I] @< 0$ } |
| |
|
| |
return fixpointmain(F,Vars)$ |
| |
} |
| |
|
| def junban(A,B){ |
def junban(A,B){ |
| return (A<B ? 1:(A>B ? -1:0))$ |
return (A<B ? 1:(A>B ? -1:0))$ |
| } |
} |
| |
|
| def worder(A,B){ |
def bsort(A){ |
| return (A[0]<B[0] ? 1:(A[0]>B[0] ? -1:0))$ |
|
| |
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 wsort(A,B,C){ |
def wsort(A,B,C,ID){ |
| |
|
| D=newvect(length(B))$ |
D=newvect(length(B))$ |
| for(I=0;I<length(B);I++) |
for(I=0;I<length(B);I++) |
| D[I]=[A[I],B[I],C[I]]$ |
D[I]=[A[I],B[I],C[I]]$ |
| |
|
| D=qsort(D,worder)$ |
D=bsort(D)$ |
| E=[]$ |
|
| for(I=0;I<length(B);I++) |
for(E=[],I=0;I<length(B);I++) |
| E=cons(D[I][1],E)$ |
E=cons(D[I][1],E)$ |
| E=reverse(E)$ |
E=reverse(E)$ |
| F=[]$ |
|
| for(I=0;I<length(B);I++) |
for(F=[],I=0;I<length(B);I++) |
| F=cons(D[I][2],F)$ |
F=cons(D[I][2],F)$ |
| F=reverse(F)$ |
F=reverse(F)$ |
| |
|
| return [E,F]$ |
return [[ID,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){ |
def nonposdegchk(Res){ |
| |
|
| for(I=0;I<length(Res);I++) |
for(I=0;I<length(Res);I++) |
| Line 60 def nonposdegchk(Res){ |
|
| Line 249 def nonposdegchk(Res){ |
|
| |
|
| def getgcd(A,B){ |
def getgcd(A,B){ |
| |
|
| VarsNumA=length(A)$ |
Anum=length(A)$ |
| VarsNumB=length(B)$ |
|
| |
|
| C=newvect(VarsNumB,B)$ |
TMP=[]$ |
| |
for(I=0;I<length(B);I++){ |
| |
|
| for(I=0;I<VarsNumA;I++){ |
for(J=0;J<Anum;J++) |
| |
if(B[I]==A[J][0]) |
| |
break; |
| |
|
| for(J=0;J<VarsNumB;J++) |
if(J==Anum) |
| if(C[J]==A[I][0]) |
TMP=cons([B[I],B[I]],TMP)$ |
| break$ |
|
| |
|
| if(J<VarsNumB) |
|
| C[J]=A[I][1]$ |
|
| } |
} |
| |
A=append(A,TMP)$ |
| |
|
| D=0$ |
Anum=length(A)$ |
| for(I=0;I<VarsNumB;I++) |
A=map(ltov,A)$ |
| D=gcd(D,C[I])$ |
|
| |
|
| |
for(D=0,I=0;I<Anum;I++) |
| |
D=gcd(D,A[I][1])$ |
| |
|
| if(D!=0){ |
if(D!=0){ |
| C=C/D$ |
for(I=0;I<Anum;I++) |
| C=map(red,C)$ |
A[I][1]=red(A[I][1]/D)$ |
| } |
} |
| |
|
| for(L=1,D=0,I=0;I<VarsNumB;I++){ |
for(L=1,D=0,I=0;I<Anum;I++){ |
| if(type(TMP=dn(C[I]))==1) |
if(type(TMP=dn(A[I][1]))==1) |
| L=ilcm(L,TMP)$ |
L=ilcm(L,TMP)$ |
| |
|
| if(type(TMP=nm(C[I]))==1) |
if(type(TMP=nm(A[I][1]))==1) |
| D=igcd(D,TMP)$ |
D=igcd(D,TMP)$ |
| } |
} |
| |
|
| C=C*L$ |
for(I=0;I<Anum;I++) |
| if(D!=0) |
A[I][1]=A[I][1]*L$ |
| C=C/D$ |
|
| |
|
| RET=[]$ |
if(D!=0){ |
| for(I=0;I<VarsNumB;I++) |
|
| RET=cons([B[I],C[I]],RET)$ |
|
| |
|
| return RET$ |
for(I=0;I<Anum;I++) |
| } |
A[I][1]=A[I][1]/D$ |
| |
|
| 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)$ |
|
| |
return map(vtol,A)$ |
| } |
} |
| |
|
| def makeret(Res,Vars){ |
def makeret(Res,Vars,FLAG){ |
| |
|
| ResNum=length(Res)$ |
ResNum=length(Res)$ |
| VarsNum=length(Vars)$ |
VarsNum=length(Vars)$ |
| |
|
| ResVec=newvect(ResNum)$ |
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(FLAG) |
| if(M==0) |
M=0$ |
| M=ResVec[I]$ |
else |
| |
M=-1$ |
| |
|
| |
for(I=0;I<ResNum;I++){ |
| |
if(member(Res[I][0],Vars)){ |
| |
ResVec[I]=Res[I][1]$ |
| |
|
| |
if(FLAG){ |
| |
if(type(ResVec[I])==1){ |
| |
if(M==0) |
| |
M=ResVec[I]$ |
| |
else |
| |
if(ResVec[I]<M) |
| |
M=ResVec[I]$ |
| |
} |
| else |
else |
| if(ResVec[I]<M) |
M=-1$ |
| M=ResVec[I]$ |
|
| } |
} |
| } |
} |
| } |
} |
| |
|
| if(M!=0) |
|
| |
if(M!=-1) |
| ResVec=ResVec/M; |
ResVec=ResVec/M; |
| |
|
| RET=newvect(VarsNum,Vars)$ |
RET=newvect(VarsNum,Vars)$ |
| |
|
| for(I=0;I<ResNum;I++){ |
for(I=0;I<ResNum;I++){ |
| for(J=0;J<VarsNum;J++) |
for(J=0;J<VarsNum;J++) |
| Line 163 def roundret(V){ |
|
| Line 350 def roundret(V){ |
|
| |
|
| VN=size(V)[0]$ |
VN=size(V)[0]$ |
| |
|
| |
K=1$ |
| RET0=V$ |
RET0=V$ |
| for(I=1;I<1000;I++){ |
RET1=map(drint,RET0)$ |
| RET1=I*RET0$ |
S=0$ |
| for(J=0;J<VN;J++){ |
for(J=0;J<VN;J++) |
| X=drint(RET1[J])$ |
S+=(RET0[J]-RET1[J])^2$ |
| if(dabs(X-RET1[J])<0.2) |
|
| RET1[J]=X$ |
for(I=2;I<10;I++){ |
| else |
RET0=I*V$ |
| break$ |
RET1=map(drint,RET0)$ |
| } |
|
| if(J==VN) |
T=0$ |
| break$ |
for(J=0;J<VN;J++) |
| |
T+=(RET0[J]-RET1[J])^2$ |
| |
|
| |
if(T<S){ |
| |
K=I$ |
| |
S=T$ |
| |
} |
| } |
} |
| |
|
| if(I==1000) |
return map(drint,K*V)$ |
| return []$ |
|
| else |
|
| return RET1$ |
|
| } |
} |
| |
|
| def chkou(L,ExpMat,CHAGORD){ |
def chkou(L,ExpMat,CHAGORD){ |
| Line 264 def checktd(PolyList,Vars,ResVars){ |
|
| Line 455 def checktd(PolyList,Vars,ResVars){ |
|
| return 1$ |
return 1$ |
| } |
} |
| |
|
| def qcheck(PolyList,Vars){ |
def value2(Vars,Ans){ |
| |
|
| |
N=length(Vars)$ |
| |
Res=newvect(N)$ |
| |
for(I=0;I<N;I++){ |
| |
Res[I]=newvect(2)$ |
| |
Res[I][0]=Vars[I]$ |
| |
Res[I][1]=Ans[I]$ |
| |
} |
| |
Res=map(vtol,Res)$ |
| |
Res=vtol(Res)$ |
| |
|
| |
Res=getgcd(Res,Vars)$ |
| |
|
| |
if(nonposdegchk(Res)){ |
| |
TMP1=makeret(Res,Vars,1)$ |
| |
return vtol(TMP1[1])$ |
| |
} |
| |
else |
| |
return []$ |
| |
} |
| |
|
| |
def qcheck(PolyList,Vars,FLAG){ |
| |
|
| |
RET=[]$ |
| Res=qcheckmain(PolyList,Vars)$ |
Res=qcheckmain(PolyList,Vars)$ |
| VarsNum=length(Vars)$ |
VarsNum=length(Vars)$ |
| |
|
| Line 286 def qcheck(PolyList,Vars){ |
|
| Line 500 def qcheck(PolyList,Vars){ |
|
| |
|
| VarsList=[]$ |
VarsList=[]$ |
| for(I=0;I<VarsNum;I++) |
for(I=0;I<VarsNum;I++) |
| if(member(Vars[CHAGORD[I]],Rea)) |
if(member(TMP0=Vars[CHAGORD[I]],Rea)) |
| VarsList=cons(Vars[CHAGORD[I]],VarsList)$ |
VarsList=cons(Vars[CHAGORD[I]],VarsList)$ |
| |
|
| Res=solve(reverse(SolveList),reverse(VarsList))$ |
Res=solve(reverse(SolveList),reverse(VarsList))$ |
| Line 294 def qcheck(PolyList,Vars){ |
|
| Line 508 def qcheck(PolyList,Vars){ |
|
| |
|
| if(nonposdegchk(Res)){ |
if(nonposdegchk(Res)){ |
| |
|
| ResVars=resvars(Res,Vars)$ |
TMP1=makeret(Res,Vars,0)$ |
| |
|
| if(checktd(PolyList,Vars,ResVars)==1){ |
if(checktd(PolyList,Vars,TMP1[1])==1){ |
| |
|
| for(J=0;J<length(Vars);J++) |
if(FLAG==0){ |
| ResVars=map(subst,ResVars,Vars[J], |
|
| strtov(rtostr(Vars[J])+"_deg"))$ |
if(TMP1[0]==0) |
| |
RET=append(RET,wsort(TMP1[1],Vars,TMP1[1],0))$ |
| return [wsort(ResVars,Vars,ResVars)]$ |
else{ |
| |
|
| |
TMP=vtol(TMP1[1])$ |
| |
RET0=[]$ |
| |
if((TMP0=fixedpoint(TMP,0))!=[]){ |
| |
|
| |
for(I=0;I<length(TMP0);I++) |
| |
TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$ |
| |
RET0=value2(Vars,TMP)$ |
| |
if(RET0!=[]) |
| |
RET0=wsort(RET0,Vars,RET0,-1)$ |
| |
} |
| |
|
| |
TMP=vtol(TMP1[1])$ |
| |
if(RET0==[] && (TMP0=fixedpoint(TMP,1))!=[]){ |
| |
|
| |
for(I=0;I<length(TMP0);I++) |
| |
TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$ |
| |
RET0=value2(Vars,TMP)$ |
| |
|
| |
if(RET0!=[]) |
| |
RET0=wsort(RET0,Vars,RET0,-1)$ |
| |
} |
| |
RET=append(RET,RET0)$ |
| |
} |
| |
} |
| |
else if(FLAG==1) |
| |
RET=append(RET,[[0,Vars,vtol(TMP1[1])]])$ |
| } |
} |
| else |
|
| return []$ |
|
| } |
} |
| else |
|
| return []$ |
|
| |
|
| |
return RET$ |
| } |
} |
| |
|
| def leastsq(NormMat,ExpMat,Vars){ |
def unitweight2(NormMat0,ExpMat,Vars,FLAG,ID){ |
| |
|
| RET=[]$ |
RET=[]$ |
| |
|
| ExpMatRowNum=size(ExpMat)[0]$ |
ExpMatRowNum=size(ExpMat)[0]$ |
| ExpMatColNum=size(ExpMat[0])[0]$ |
ExpMatColNum=size(ExpMat[0])[0]$ |
| |
ExtMatColNum=ExpMatColNum+1$ |
| |
|
| |
ExtVars=append(Vars,[uc()])$ |
| |
|
| if(NormMat==0){ |
if(NormMat==0){ |
| NormMat=newmat(ExpMatColNum,ExpMatColNum)$ |
|
| |
|
| |
NormMat0=newvect(ExtMatColNum)$ |
| |
for(I=0;I<ExtMatColNum;I++) |
| |
NormMat0[I]=newvect(ExtMatColNum)$ |
| |
|
| for(I=0;I<ExpMatColNum;I++) |
for(I=0;I<ExpMatColNum;I++) |
| for(J=I;J<ExpMatColNum;J++) |
for(J=I;J<ExpMatColNum;J++) |
| for(K=0;K<ExpMatRowNum;K++) |
for(K=0;K<ExpMatRowNum;K++) |
| NormMat[I][J]+= |
NormMat0[I][J]+= |
| ExpMat[K][I]*ExpMat[K][J]$ |
ExpMat[K][I]* |
| |
ExpMat[K][J]$ |
| } |
} |
| |
|
| BVec=newvect(ExpMatColNum)$ |
|
| |
|
| for(I=0;I<ExpMatColNum;I++) |
for(I=0;I<ExpMatColNum;I++) |
| for(J=0;J<ExpMatRowNum;J++) |
for(K=0;K<ExpMatRowNum;K++) |
| BVec[I]+=ExpMat[J][I]$ |
NormMat0[I][ExpMatColNum]-=ExpMat[K][I]$ |
| |
|
| SolveList=[]$ |
NormMat0[ExpMatColNum][ExpMatColNum]=ExpMatRowNum$ |
| 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++) |
WorkMat=newvect(ExtMatColNum)$ |
| TMP+=NormMat[I][J]*Vars[J]$ |
for(I=0;I<ExtMatColNum;I++) |
| |
WorkMat[I]=newvect(ExtMatColNum)$ |
| |
|
| TMP-=BVec[I]$ |
if(jacobi(ExtMatColNum,NormMat0,WorkMat)){ |
| SolveList=cons(TMP,SolveList)$ |
|
| } |
|
| |
|
| Rea=vars(SolveList)$ |
Res=newvect(ExtMatColNum)$ |
| |
for(I=0;I<ExtMatColNum;I++){ |
| |
Res[I]=newvect(2)$ |
| |
Res[I][0]=ExtVars[I]$ |
| |
Res[I][1]=WorkMat[ExtMatColNum-1][I]$ |
| |
} |
| |
|
| VarsList=[]$ |
if(nonposdegchk(Res)){ |
| for(I=0;I<length(Vars);I++) |
|
| if(member(Vars[I],Rea)) |
|
| VarsList=cons(Vars[I],VarsList)$ |
|
| |
|
| Res=solve(SolveList,VarsList)$ |
TMP1=makeret(Res,Vars,1)$ |
| Res=getgcd(Res,Rea)$ |
|
| |
|
| if(nonposdegchk(Res)){ |
if(FLAG==0){ |
| TMP1=makeret(Res,Vars)$ |
TMP=roundret(TMP1[1])$ |
| 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, |
RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]*1.0),ID))$ |
| map(drint,TMP1[1]*1.0)),RET)$ |
|
| |
|
| return RET$ |
if(TMP!=[]) |
| |
RET=append(RET,wsort(TMP1[1],Vars,TMP,ID+1))$ |
| |
} |
| |
else if(FLAG==1) |
| |
RET=append(RET,[[ID,Vars,vtol(TMP1[1])]])$ |
| } |
} |
| else{ |
} |
| RET=cons(wsort(TMP1[1],Vars,TMP1[1]*1.0),RET)$ |
|
| return RET$ |
return RET$ |
| } |
|
| } |
|
| else |
|
| return RET$ |
|
| |
|
| } |
} |
| |
|
| def weightr(ExpMat,Vars,PolyListNum,OneMat){ |
def unitweight1(ExpMat,Vars,PolyListNum,OneMat,FLAG){ |
| |
|
| RET=[]$ |
RET=[]$ |
| |
|
| Line 394 def weightr(ExpMat,Vars,PolyListNum,OneMat){ |
|
| Line 627 def weightr(ExpMat,Vars,PolyListNum,OneMat){ |
|
| |
|
| ExtVars=reverse(ExtVars)$ |
ExtVars=reverse(ExtVars)$ |
| |
|
| NormMat=newmat(ExpMatColNum,ExtMatColNum)$ |
NormMat0=newvect(ExpMatColNum+1)$ |
| |
for(I=0;I<ExpMatColNum;I++) |
| |
NormMat0[I]=newvect(ExpMatColNum+1)$ |
| |
|
| for(I=0;I<ExpMatColNum;I++) |
for(I=0;I<ExpMatColNum;I++) |
| for(J=I;J<ExpMatColNum;J++) |
for(J=I;J<ExpMatColNum;J++) |
| for(K=0;K<ExpMatRowNum;K++) |
for(K=0;K<ExpMatRowNum;K++) |
| NormMat[I][J]+= |
NormMat0[I][J]+= |
| ExpMat[K][I]* |
ExpMat[K][I]* |
| ExpMat[K][J]$ |
ExpMat[K][J]$ |
| |
|
| |
NormMat1=newvect(ExtMatColNum)$ |
| |
for(I=0;I<ExtMatColNum;I++) |
| |
NormMat1[I]=newvect(ExtMatColNum)$ |
| |
|
| |
WorkMat=newvect(ExtMatColNum)$ |
| |
for(I=0;I<ExtMatColNum;I++) |
| |
WorkMat[I]=newvect(ExtMatColNum)$ |
| |
|
| for(I=0;I<ExpMatColNum;I++) |
for(I=0;I<ExpMatColNum;I++) |
| |
for(J=I;J<ExpMatColNum;J++) |
| |
NormMat1[I][J]=NormMat0[I][J]$ |
| |
|
| |
for(I=0;I<ExpMatColNum;I++) |
| for(J=0;J<PolyListNum;J++) |
for(J=0;J<PolyListNum;J++) |
| for(K=OneMat[J];K<OneMat[J+1];K++) |
for(K=OneMat[J];K<OneMat[J+1];K++) |
| NormMat[I][J+ExpMatColNum]-= |
NormMat1[I][J+ExpMatColNum]-=ExpMat[K][I]$ |
| ExpMat[K][I]$ |
|
| |
|
| WVect=newvect(PolyListNum)$ |
|
| for(I=0;I<PolyListNum;I++) |
for(I=0;I<PolyListNum;I++) |
| WVect[I]=OneMat[I+1]-OneMat[I]$ |
NormMat1[I+ExpMatColNum][I+ExpMatColNum]=OneMat[I+1]-OneMat[I]$ |
| |
|
| for(F=0;F<ExtMatColNum;F++){ |
if(jacobi(ExtMatColNum,NormMat1,WorkMat)){ |
| SolveList=[]$ |
|
| for(I=0;I<ExpMatColNum;I++){ |
|
| if (F==I) |
|
| continue$ |
|
| |
|
| TMP=0$ |
Res=newvect(ExtMatColNum)$ |
| |
for(I=0;I<ExtMatColNum;I++){ |
| for(J=0;J<I;J++) |
Res[I]=newvect(2)$ |
| if(J!=F) |
Res[I][0]=ExtVars[I]$ |
| TMP+=NormMat[J][I]*ExtVars[J]$ |
Res[I][1]=WorkMat[ExtMatColNum-1][I]$ |
| |
} |
| |
|
| for(J=I;J<ExtMatColNum;J++) |
if(nonposdegchk(Res)){ |
| 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]$ |
TMP1=makeret(Res,Vars,1)$ |
| |
|
| if(F<ExpMatColNum) |
if(FLAG==0){ |
| TMP+=NormMat[F][I+ExpMatColNum]$ |
TMP=roundret(TMP1[1])$ |
| |
|
| SolveList=cons(TMP,SolveList)$ |
RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]*1.0),1))$ |
| } |
|
| |
|
| 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)$ |
|
| |
|
| TMP=cons(ExtVars[F],TMP)$ |
|
| Res=getgcd(Res,TMP)$ |
|
| |
|
| if(nonposdegchk(Res)){ |
|
| |
|
| TMP1=makeret(Res,Vars)$ |
|
| if(TMP1[0]==0){ |
|
| TMP=roundret(TMP1[1]*1.0)$ |
|
| if(TMP!=[]) |
if(TMP!=[]) |
| RET=cons(wsort(TMP1[1],Vars,TMP),RET)$ |
RET=append(RET,wsort(TMP1[1],Vars,TMP,2))$ |
| |
|
| RET=cons(wsort(TMP1[1],Vars, |
|
| map(drint,TMP1[1]*1.0)),RET)$ |
|
| } |
} |
| else{ |
else if(FLAG==1) |
| RET=cons(wsort(TMP1[1],Vars,TMP1[1]*1.0),RET)$ |
RET=append(RET,[[1,Vars,vtol(TMP1[1])]])$ |
| } |
|
| } |
} |
| |
} |
| } |
|
| |
return [NormMat0,RET]$ |
| return [NormMat,RET]$ |
|
| } |
} |
| |
|
| def weight(PolyList,Vars,FLAG){ |
def weight(PolyList,Vars,FLAG){ |
| Line 506 def weight(PolyList,Vars,FLAG){ |
|
| Line 699 def weight(PolyList,Vars,FLAG){ |
|
| |
|
| RET=[]$ |
RET=[]$ |
| |
|
| TMP=qcheck(PolyList,Vars)$ |
TMP=qcheck(PolyList,Vars,FLAG)$ |
| |
|
| if(TMP!=[]){ |
if(TMP!=[]){ |
| RET=append(RET,TMP)$ |
RET=append(RET,TMP)$ |
| return cons(1,RET)$ |
return RET$ |
| } |
} |
| |
|
| dp_ord(2)$ |
dp_ord(2)$ |
| |
|
| PolyListNum=length(PolyList)$ |
PolyListNum=length(PolyList)$ |
| |
|
| if(FLAG){ |
OneMat=newvect(PolyListNum+1,[0])$ |
| |
ExpMat=[]$ |
| OneMat=newvect(PolyListNum+1,[0])$ |
for(I=0;I<PolyListNum;I++){ |
| ExpMat=[]$ |
for(Poly=dp_ptod(PolyList[I],Vars); |
| for(I=0;I<PolyListNum;I++){ |
Poly!=0;Poly=dp_rest(Poly)){ |
| for(Poly=dp_ptod(PolyList[I],Vars); |
ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$ |
| Poly!=0;Poly=dp_rest(Poly)){ |
|
| ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$ |
|
| } |
|
| OneMat[I+1]=length(ExpMat)$ |
|
| } |
} |
| |
OneMat[I+1]=length(ExpMat)$ |
| |
} |
| |
|
| ExpMat=reverse(ExpMat)$ |
ExpMat=reverse(ExpMat)$ |
| ExpMat=newvect(length(ExpMat),ExpMat)$ |
ExpMat=newvect(length(ExpMat),ExpMat)$ |
| |
|
| TMP=weightr(ExpMat,Vars,PolyListNum,OneMat)$ |
TMP=unitweight1(ExpMat,Vars,PolyListNum,OneMat,FLAG)$ |
| RET=append(RET,TMP[1])$ |
|
| RET=append(RET,leastsq(TMP[0],ExpMat,Vars))$ |
|
| } |
|
| 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)$ |
RET=append(RET,TMP[1])$ |
| ExpMat=newvect(length(ExpMat),ExpMat)$ |
|
| |
|
| RET=append(RET,leastsq(0,ExpMat,Vars))$ |
TMP=unitweight2(TMP[0],ExpMat,Vars,FLAG,3)$ |
| } |
|
| |
|
| |
RET=append(RET,TMP)$ |
| |
|
| ExpMat=qsort(ExpMat,junban)$ |
ExpMat=qsort(ExpMat,junban)$ |
| |
|
| ExpMat2=[]$ |
ExpMat2=[]$ |
| Line 561 def weight(PolyList,Vars,FLAG){ |
|
| Line 740 def weight(PolyList,Vars,FLAG){ |
|
| |
|
| if(size(ExpMat)[0]!=length(ExpMat2)){ |
if(size(ExpMat)[0]!=length(ExpMat2)){ |
| ExpMat=newvect(length(ExpMat2),ExpMat2)$ |
ExpMat=newvect(length(ExpMat2),ExpMat2)$ |
| RET=append(RET,leastsq(0,ExpMat,Vars))$ |
RET=append(RET,unitweight2(0,ExpMat,Vars,FLAG,5))$ |
| } |
} |
| |
else{ |
| |
TMP0=map(ltov,TMP0)$ |
| |
|
| RET=derase(RET)$ |
for(I=0;I<length(TMP0);I++) |
| return cons(0,RET)$ |
if(TMP0[I][0]==3) |
| |
TMP0[I][0]=5$ |
| |
else if(TMP0[I][0]==4) |
| |
TMP0[I][0]=6$ |
| |
|
| |
TMP0=map(vtol,TMP0)$ |
| |
|
| |
RET=append(RET,TMP0)$ |
| |
} |
| |
|
| |
return RET$ |
| } |
} |
| |
|
| end$ |
end$ |