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