![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to weight CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / lib|
Return to weight CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
1.4 kimura 1: load("solve")$
2: load("gr")$
3:
4: def nonposdegchk(Res){
5:
6: for(I=0;I<length(Res);I++)
7: if(Res[I][1]<=0)
8: return 0$
9:
10: return 1$
11: }
12:
1.6 ! kimura 13:
! 14: def notzerovec(Vec){
! 15:
! 16: for(I=0;I<size(Vec)[0];I++)
! 17: if(Vec[I]!=0)
! 18: return 1$
! 19:
! 20: return 0$
! 21: }
! 22:
1.4 kimura 23: def resvars(Res,Vars){
24:
25: ResVars=newvect(length(Vars),Vars)$
26: for(I=0;I<length(Res);I++){
27:
28: for(J=0;J<size(ResVars)[0];J++)
29: if(Res[I][0]==ResVars[J])
30: break$
31:
32: if(J<size(ResVars)[0])
33: ResVars[J]=Res[I][1]$
34: }
35: return(ResVars)$
36: }
37:
38: def makeret1(Res,Vars){
39:
40: VarsNum=length(Vars)$
41:
42: ResVec=newvect(VarsNum,Vars)$
43:
1.6 ! kimura 44: for(F=0,I=0,M=0;I<length(Res);I++){
1.4 kimura 45:
46: for(J=0;J<VarsNum;J++)
47: if(Res[I][0]==Vars[J])
48: break$
49:
50: if(J<VarsNum){
51: ResVec[J]=Res[I][1]$
52:
1.6 ! kimura 53: if(F==0 && type(ResVec[J])==1){
1.4 kimura 54: if(M==0)
55: M=ResVec[J]$
56: else
57: if(ResVec[J]<M)
58: M=ResVec[J]$
59: }
1.6 ! kimura 60: else
! 61: F=1$
1.4 kimura 62: }
63:
64: }
65:
66: if(F==0)
67: for(I=0;I<VarsNum;I++)
68: ResVec[I]=ResVec[I]/M*1.0$
69:
70: for(I=0;I<VarsNum;I++)
71: for(J=0;J<length(Vars);J++)
72: ResVec[I]=subst(ResVec[I],Vars[J],
73: strtov(rtostr(Vars[J])+"_deg"))$
74:
75: ResVec=cons(F,vtol(ResVec))$
76: return ResVec$
77: }
78:
1.5 kimura 79: def junban(A,B){
1.4 kimura 80:
81: for(I=0;I<size(A)[0];I++){
82: if(A[I]<B[I])
83: return 1$
84:
85: if(A[I]>B[I])
86: return -1$
87: }
88:
89: return 0$
90: }
91:
92: def roundret(V){
93:
94: VN=length(V)$
95: RET0=newvect(VN,V)$
96:
97: for(I=1;I<1000;I++){
98: RET1=I*RET0$
99: for(J=0;J<VN;J++){
100: X=drint(RET1[J])$
101: if(dabs(X-RET1[J])<0.2)
102: RET1[J]=X$
103: else
104: break$
105: }
106: if(J==VN)
107: break$
108: }
109:
110: if(I==1000)
111: return []$
112: else
113: return RET1$
114: }
115:
116: def chkou(L,ExpMat,CHAGORD){
117:
118: P=1$
119: F=ExpMat[L]$
120:
121: for(I=0;I<L;I++){
122: Q=ExpMat[L][CHAGORD[I]]$
123: for(J=0;J<size(ExpMat[0])[0];J++){
124: ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
125: *ExpMat[L][CHAGORD[J]]-
126: Q*ExpMat[I][CHAGORD[J]])/P)$
127: }
128:
129: P=ExpMat[I][CHAGORD[I]]$
130: }
131:
132: for(J=0;J<size(ExpMat[0])[0];J++)
133: if(ExpMat[L][CHAGORD[J]]!=0)
134: break$
135:
136: if(J==size(ExpMat[0])[0])
137: return L$
138: else{
139: TMP=CHAGORD[L]$
140: CHAGORD[L]=CHAGORD[J]$
141: CHAGORD[J]=TMP$
142: return (L+1)$
143: }
144: }
145:
146: def qcheck0(PolyList,Vars){
147:
148: RET=[]$
149: PolyListNum=length(PolyList)$
150: VarsNum=length(Vars)$
151:
152: ExpMat=newvect(VarsNum)$
153: CHAGORD=newvect(VarsNum)$
154: for(I=0;I<VarsNum;I++)
155: CHAGORD[I]=I$
156:
157: L=0$
158: for(I=0;I<PolyListNum;I++){
159: Poly=dp_ptod(PolyList[I],Vars)$
160: BASE0=dp_etov(dp_ht(Poly))$
161: Poly=dp_rest(Poly)$
162: for(;Poly!=0;Poly=dp_rest(Poly)){
163: ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
164: L=chkou(L,ExpMat,CHAGORD)$
165: if(L==VarsNum-1){
166: RET=cons(ExpMat,RET)$
167: RET=cons(CHAGORD,RET)$
168: RET=cons(L,RET)$
169: return RET$
170: }
171: }
172: }
173:
174: RET=cons(ExpMat,RET)$
175: RET=cons(CHAGORD,RET)$
176: RET=cons(L,RET)$
177: return RET$
178: }
179:
180: def inner(A,B){
181:
182: SUM=0$
183: for(I=0;I<size(A)[0];I++)
184: SUM+=A[I]*B[I]$
185:
186: return SUM$
187: }
188:
189: def checktd(PolyList,Vars,ResVars){
190:
191: PolyListNum=length(PolyList)$
192: VarsNum=length(Vars)$
193:
194: L=0$
195: for(I=0;I<PolyListNum;I++){
196: Poly=dp_ptod(PolyList[I],Vars)$
197: J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
198: Poly=dp_rest(Poly)$
199: for(;Poly!=0;Poly=dp_rest(Poly))
200: if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
201: return 0$
202: }
203:
204: return 1$
205: }
206:
207: def getgcd(A,B){
208:
209: VarsNumA=length(A)$
210: VarsNumB=length(B)$
211:
212: C=newvect(VarsNumB,B)$
213:
214: for(I=0;I<VarsNumA;I++){
215:
216: for(J=0;J<VarsNumB;J++)
217: if(C[J]==A[I][0])
218: break$
219:
220: C[J]=A[I][1]$
221: }
222:
223: D=0$
224: for(I=0;I<VarsNumB;I++)
225: D=gcd(D,C[I])$
226:
227: if(D!=0){
228:
229: for(I=0;I<VarsNumB;I++)
230: C[I]=red(C[I]/D)$
231:
232: }
233:
1.6 ! kimura 234: for(L=1,D=0,I=0;I<VarsNumB;I++){
! 235: if(type(TMP=dn(C[I]))==1)
! 236: L=ilcm(L,TMP)$
! 237:
! 238: if(type(TMP=nm(C[I]))==1)
! 239: D=igcd(D,TMP)$
! 240: }
1.4 kimura 241:
242: for(I=0;I<VarsNumB;I++)
243: C[I]=C[I]*L$
244:
245: if(D!=0)
246: for(I=0;I<VarsNumB;I++)
247: C[I]=C[I]/D$
248:
249:
250: RET=newvect(VarsNumB)$
251: for(I=0;I<VarsNumB;I++){
252: RET[I]=newvect(2)$
253: RET[I][0]=B[I]$
254: RET[I][1]=C[I]$
255: }
256:
257: return vtol(map(vtol,RET))$
258: }
259:
260: def qcheck(PolyList,Vars){
261:
262: RET=[]$
263: Res=qcheck0(PolyList,Vars)$
264: VarsNum=length(Vars)$
265:
266: IndNum=Res[0]$
267: CHAGORD=Res[1]$
268: ExpMat=Res[2]$
269:
270: SolveList=[]$
271: for(I=0;I<IndNum;I++){
272: TMP=0$
273: for(J=0;J<VarsNum;J++)
274: TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$
275:
276: SolveList=cons(TMP,SolveList)$
277: }
278:
279: VarsList=[]$
280: for(I=0;I<VarsNum;I++)
281: VarsList=cons(Vars[CHAGORD[I]],VarsList)$
282:
283: Rea=vars(SolveList)$
284: Res=solve(reverse(SolveList),reverse(VarsList))$
285:
286: if(nonposdegchk(Res)){
287:
288: Res=getgcd(Res,Rea)$
289: ResVars=resvars(Res,Vars)$
290:
291: if(checktd(PolyList,Vars,ResVars)==1){
292:
293: for(J=0;J<length(Vars);J++)
294: ResVars=map(subst,ResVars,Vars[J],
295: strtov(rtostr(Vars[J])+"_deg"))$
296:
297: RET=cons([vtol(ResVars),ResVars,[]],RET)$
298: return cons(1,RET)$
299: }
300: else
301: return cons(0,RET)$
302: }
303: else
304: return cons(0,RET)$
305:
306: }
307:
1.6 ! kimura 308: def leastsq(ExpMat,Vars){
1.4 kimura 309:
310: ExpMatRowNum=size(ExpMat)[0]$
311: ExpMatColNum=size(ExpMat[0])[0]$
312:
313: NormMat=newmat(ExpMatColNum,ExpMatColNum+1)$
314:
315: for(I=0;I<ExpMatColNum;I++)
316: for(J=0;J<ExpMatColNum;J++)
317: for(K=0;K<ExpMatRowNum;K++)
318: NormMat[I][J]+=ExpMat[K][I]*ExpMat[K][J]$
319:
320: for(I=0;I<ExpMatColNum;I++)
321: for(J=0;J<ExpMatRowNum;J++)
322: NormMat[I][ExpMatColNum]+=ExpMat[J][I]$
323:
324: SolveList=[]$
325: for(I=0;I<ExpMatColNum;I++){
326: TMP=0$
327: for(J=0;J<ExpMatColNum;J++)
328: TMP+=NormMat[I][J]*Vars[J]$
329:
330: TMP-=NormMat[I][ExpMatColNum]$
331: SolveList=cons(TMP,SolveList)$
332: }
333:
334: Rea=vars(SolveList)$
335: Res=solve(SolveList,Vars)$
336:
337: if(nonposdegchk(Res)){
338: Res=getgcd(Res,Rea)$
339: TMP1=makeret1(Res,Vars);
340: if(car(TMP1)==0){
341: TMP2=roundret(cdr(TMP1));
342: TMP3=map(drint,cdr(TMP1))$
1.6 ! kimura 343: return([cdr(TMP1),newvect(length(TMP3),TMP3),TMP2])$
1.4 kimura 344: }
345: else
1.6 ! kimura 346: return([cdr(TMP1),[],[]])$
1.4 kimura 347: }
348:
1.6 ! kimura 349: }
! 350:
! 351: def weight(PolyList,Vars){
! 352:
! 353: Vars0=vars(PolyList)$
! 354: Vars1=[]$
! 355: for(I=0;I<length(Vars);I++)
! 356: if(member(Vars[I],Vars0))
! 357: Vars1=cons(Vars[I],Vars1)$
! 358:
! 359: Vars=reverse(Vars1)$
! 360:
! 361: RET=[]$
! 362:
! 363: TMP=qcheck(PolyList,Vars)$
! 364:
! 365: if(car(TMP)==1){
! 366: RET=cdr(TMP)$
! 367: RET=cons(Vars,RET)$
! 368: RET=cons(1,RET)$
! 369: return RET$
! 370: }
1.4 kimura 371:
1.6 ! kimura 372: dp_ord(2)$
1.4 kimura 373:
1.6 ! kimura 374: PolyListNum=length(PolyList)$
! 375: VPolyList=newvect(PolyListNum,PolyList)$
1.4 kimura 376:
1.6 ! kimura 377: ExpMat=[]$
! 378: for(I=0;I<PolyListNum;I++)
! 379: for(Poly=dp_ptod(VPolyList[I],Vars);
! 380: Poly!=0;Poly=dp_rest(Poly)){
! 381: Exp=dp_etov(dp_ht(Poly))$
! 382: if(notzerovec(Exp))
! 383: ExpMat=cons(Exp,ExpMat)$
! 384: }
1.4 kimura 385:
1.6 ! kimura 386: ExpMat=reverse(ExpMat)$
! 387: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 388:
1.6 ! kimura 389: /* first */
1.4 kimura 390:
1.6 ! kimura 391: RET=cons(leastsq(ExpMat,Vars),RET)$
1.4 kimura 392:
1.6 ! kimura 393: /* second */
1.4 kimura 394:
1.6 ! kimura 395: ExpMat=qsort(ExpMat,junban)$
! 396: ExpMat2=[]$
! 397: for(I=0;I<size(ExpMat)[0];I++)
! 398: if(car(ExpMat2)!=ExpMat[I])
! 399: ExpMat2=cons(ExpMat[I],ExpMat2)$
1.4 kimura 400:
1.6 ! kimura 401: if(size(ExpMat)[0]!=length(ExpMat2)){
! 402: ExpMat=newvect(length(ExpMat2),ExpMat2)$
! 403: RET=cons(leastsq(ExpMat,Vars),RET)$
1.4 kimura 404: }
405:
406: RET=cons(Vars,reverse(RET))$
407: RET=cons(0,RET)$
408: return RET$
409: }
410:
411: end$