![[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:
1.9 kimura 4: def nonzerovec(A){
5:
6: for(I=0;I<size(A)[0];I++)
7: if(A[I]!=0)
8: return 1$
9:
10: return 0$
11: }
12:
1.8 kimura 13: def junban(A,B){
14: return (A<B ? 1:(A>B ? -1:0))$
15: }
16:
17: def worder(A,B){
18: return (A[0]<B[0] ? 1:(A[0]>B[0] ? -1:0))$
19: }
20:
1.10 kimura 21: def bsort(A){
1.8 kimura 22:
1.10 kimura 23: K=size(A)[0]-1$
24: while(K>=0){
25: J=-1$
26: for(I=1;I<=K;I++)
27: if(A[I-1][0]<A[I][0]){
28: J=I-1$
29: X=A[J]$
30: A[J]=A[I]$
31: A[I]=X$
32: }
33: K=J$
34: }
35: return A$
36: }
37:
38: def perm(I,P,TMP){
39:
40: if(I>0){
41: TMP=perm(I-1,P,TMP)$
42: for(J=I-1;J>=0;J--){
43: T=P[I]$
44: P[I]=P[J]$
45: P[J]=T$
46: TMP=perm(I-1,P,TMP)$
47: T=P[I]$
48: P[I]=P[J]$
49: P[J]=T$
50: }
51:
52: return TMP$
53: }
54: else{
55: for(TMP0=[],K=0;K<size(P)[0];K++)
56: TMP0=cons(P[K],TMP0)$
57:
58: TMP=cons(TMP0,TMP)$
59: return TMP$
60: }
61: }
62:
63: def marge(A,B){
64:
65: RET=[]$
66: for(I=0;I<length(A);I++)
67: for(J=0;J<length(B);J++)
68: RET=cons(append(A[I],B[J]),RET)$
69:
70: return RET$
71: }
72:
73: def wsort(A,B,C,FLAG){
74:
75: if(FLAG==0){
76: D=newvect(length(B))$
77: for(I=0;I<length(B);I++)
78: D[I]=[A[I],B[I],C[I]]$
79:
80: D=bsort(D)$
81: E=[]$
82: for(I=0;I<length(B);I++)
83: E=cons(D[I][1],E)$
84: E=reverse(E)$
85: F=[]$
86: for(I=0;I<length(B);I++)
87: F=cons(D[I][2],F)$
88: F=reverse(F)$
1.8 kimura 89:
1.10 kimura 90: return [[E,F]]$
91: }
92: else{
93: D=newvect(length(B))$
94: for(I=0;I<length(B);I++)
95: D[I]=[A[I],B[I],C[I]]$
96:
97: D=qsort(D,worder)$
98: D0=[]$
99:
100: for(I=0,J=0,TMP=[],X=0;I<size(D)[0];I++){
101: if(X==D[I][0])
102: TMP=cons(cdr(D[I]),TMP)$
103: else{
104: D0=cons(TMP,D0)$
105: TMP=[]$
106: TMP=cons(cdr(D[I]),TMP)$
107: X=car(D[I])$
108: }
109: }
110: D0=cdr(reverse(cons(TMP,D0)))$
111: D0=map(ltov,D0)$
112: for(I=0,TMP=[[]];I<length(D0);I++){
113: TMP0=perm(length(D0[I])-1,D0[I],[])$
114: TMP=marge(TMP,TMP0)$
115: }
116:
117: RET=[]$
118: for(I=0;I<length(TMP);I++){
119: TMP0=[]$
120: TMP1=[]$
121: for(J=0;J<length(TMP[I]);J++){
122: TMP0=cons(TMP[I][J][0],TMP0)$
123: TMP1=cons(TMP[I][J][1],TMP1)$
124: }
125: TMP0=reverse(TMP0)$
126: TMP1=reverse(TMP1)$
127:
128: RET=cons([TMP0,TMP1],RET)$
129: }
130:
131: return RET$
132: }
1.8 kimura 133: }
134:
135: def derase(A){
136:
137: B=newvect(length(A),A)$
138: B=qsort(B,junban)$
139: C=[]$
140: for(I=0;I<size(B)[0];I++)
141: if(car(C)!=B[I])
142: C=cons(B[I],C)$
143:
144: return reverse(C)$
145: }
146:
1.4 kimura 147: def nonposdegchk(Res){
148:
149: for(I=0;I<length(Res);I++)
150: if(Res[I][1]<=0)
151: return 0$
152:
153: return 1$
154: }
155:
1.8 kimura 156: def getgcd(A,B){
157:
158: VarsNumA=length(A)$
159: VarsNumB=length(B)$
160:
161: C=newvect(VarsNumB,B)$
162:
163: for(I=0;I<VarsNumA;I++){
164:
165: for(J=0;J<VarsNumB;J++)
1.16 kimura 166: if(B[J]==A[I][0])
1.8 kimura 167: break$
168:
169: if(J<VarsNumB)
170: C[J]=A[I][1]$
171: }
172:
173: D=0$
174: for(I=0;I<VarsNumB;I++)
175: D=gcd(D,C[I])$
176:
177: if(D!=0){
178: C=C/D$
179: C=map(red,C)$
180: }
1.6 kimura 181:
1.8 kimura 182: for(L=1,D=0,I=0;I<VarsNumB;I++){
183: if(type(TMP=dn(C[I]))==1)
184: L=ilcm(L,TMP)$
185:
186: if(type(TMP=nm(C[I]))==1)
187: D=igcd(D,TMP)$
188: }
189:
190: C=C*L$
191: if(D!=0)
192: C=C/D$
1.6 kimura 193:
1.8 kimura 194: RET=[]$
195: for(I=0;I<VarsNumB;I++)
196: RET=cons([B[I],C[I]],RET)$
1.6 kimura 197:
1.8 kimura 198: return RET$
1.6 kimura 199: }
200:
1.10 kimura 201: def makeret(Res,Vars,FLAG){
1.4 kimura 202:
1.8 kimura 203: ResNum=length(Res)$
1.4 kimura 204: VarsNum=length(Vars)$
205:
1.17 ! kimura 206: ResVec=newvect(ResNum)$
1.4 kimura 207:
1.17 ! kimura 208: for(M=0,I=0;I<ResNum;I++){
! 209: if(member(Res[I][0],Vars)){
! 210: ResVec[I]=Res[I][1]$
1.16 kimura 211:
1.17 ! kimura 212: if(FLAG && type(ResVec[I])==1){
! 213: if(M==0)
! 214: M=ResVec[I]$
! 215: else
! 216: if(ResVec[I]<M)
! 217: M=ResVec[I]$
! 218: }
! 219: }
! 220: }
1.16 kimura 221:
1.8 kimura 222: if(M!=0)
223: ResVec=ResVec/M;
1.4 kimura 224:
1.17 ! kimura 225: RET=newvect(VarsNum,Vars)$
1.4 kimura 226:
1.8 kimura 227: for(I=0;I<ResNum;I++){
228: for(J=0;J<VarsNum;J++)
229: if(Vars[J]==Res[I][0])
230: break$
1.4 kimura 231:
1.8 kimura 232: if(J<VarsNum)
233: RET[J]=ResVec[I]$
234: }
235:
1.10 kimura 236:
237: for(J=0;J<length(Vars);J++)
238: RET=map(subst,RET,Vars[J],
239: strtov(rtostr(Vars[J])+"_deg"))$
240:
1.4 kimura 241: for(I=0;I<VarsNum;I++)
1.8 kimura 242: if(type(RET[I])!=1)
243: return [1,RET]$
1.4 kimura 244:
1.8 kimura 245: return [0,RET]$
1.4 kimura 246: }
247:
248: def roundret(V){
249:
1.8 kimura 250: VN=size(V)[0]$
1.4 kimura 251:
1.8 kimura 252: RET0=V$
1.13 kimura 253: for(I=1;I<1000;I++){
1.4 kimura 254: RET1=I*RET0$
255: for(J=0;J<VN;J++){
256: X=drint(RET1[J])$
257: if(dabs(X-RET1[J])<0.2)
258: RET1[J]=X$
259: else
260: break$
261: }
262: if(J==VN)
263: break$
264: }
265:
266: if(I==1000)
267: return []$
268: else
269: return RET1$
270: }
271:
272: def chkou(L,ExpMat,CHAGORD){
273:
1.8 kimura 274: for(P=1,I=0;I<L;I++){
1.4 kimura 275: Q=ExpMat[L][CHAGORD[I]]$
276: for(J=0;J<size(ExpMat[0])[0];J++){
277: ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
278: *ExpMat[L][CHAGORD[J]]-
279: Q*ExpMat[I][CHAGORD[J]])/P)$
280: }
281:
282: P=ExpMat[I][CHAGORD[I]]$
283: }
284:
285: for(J=0;J<size(ExpMat[0])[0];J++)
286: if(ExpMat[L][CHAGORD[J]]!=0)
287: break$
288:
289: if(J==size(ExpMat[0])[0])
290: return L$
291: else{
292: TMP=CHAGORD[L]$
293: CHAGORD[L]=CHAGORD[J]$
294: CHAGORD[J]=TMP$
295: return (L+1)$
296: }
297: }
298:
1.8 kimura 299: def qcheckmain(PolyList,Vars){
1.4 kimura 300:
301: RET=[]$
302: PolyListNum=length(PolyList)$
303: VarsNum=length(Vars)$
304:
305: ExpMat=newvect(VarsNum)$
306: CHAGORD=newvect(VarsNum)$
307: for(I=0;I<VarsNum;I++)
308: CHAGORD[I]=I$
309:
310: L=0$
311: for(I=0;I<PolyListNum;I++){
312: Poly=dp_ptod(PolyList[I],Vars)$
313: BASE0=dp_etov(dp_ht(Poly))$
314: Poly=dp_rest(Poly)$
315: for(;Poly!=0;Poly=dp_rest(Poly)){
316: ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
317: L=chkou(L,ExpMat,CHAGORD)$
1.8 kimura 318: if(L==VarsNum-1)
319: return [L,CHAGORD,ExpMat]$
1.4 kimura 320: }
321: }
322:
1.8 kimura 323: return [L,CHAGORD,ExpMat]$
1.4 kimura 324: }
325:
326: def inner(A,B){
327:
328: SUM=0$
329: for(I=0;I<size(A)[0];I++)
330: SUM+=A[I]*B[I]$
331:
332: return SUM$
333: }
334:
335: def checktd(PolyList,Vars,ResVars){
336:
337: PolyListNum=length(PolyList)$
338: VarsNum=length(Vars)$
339:
340: L=0$
341: for(I=0;I<PolyListNum;I++){
342: Poly=dp_ptod(PolyList[I],Vars)$
343: J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
344: Poly=dp_rest(Poly)$
345: for(;Poly!=0;Poly=dp_rest(Poly))
346: if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
347: return 0$
348: }
349:
350: return 1$
351: }
352:
1.10 kimura 353: def qcheck(PolyList,Vars,FLAG){
1.4 kimura 354:
1.10 kimura 355: RET=[]$
1.8 kimura 356: Res=qcheckmain(PolyList,Vars)$
1.4 kimura 357: VarsNum=length(Vars)$
358:
359: IndNum=Res[0]$
360: CHAGORD=Res[1]$
361: ExpMat=Res[2]$
362:
363: SolveList=[]$
364: for(I=0;I<IndNum;I++){
365: TMP=0$
366: for(J=0;J<VarsNum;J++)
367: TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$
368:
369: SolveList=cons(TMP,SolveList)$
370: }
371:
1.8 kimura 372: Rea=vars(SolveList)$
373:
1.4 kimura 374: VarsList=[]$
375: for(I=0;I<VarsNum;I++)
1.8 kimura 376: if(member(Vars[CHAGORD[I]],Rea))
377: VarsList=cons(Vars[CHAGORD[I]],VarsList)$
1.4 kimura 378:
379: Res=solve(reverse(SolveList),reverse(VarsList))$
1.8 kimura 380: Res=getgcd(Res,Rea)$
1.4 kimura 381:
382: if(nonposdegchk(Res)){
383:
1.10 kimura 384: ResVars=makeret(Res,Vars,0)$
1.4 kimura 385:
1.10 kimura 386: if(checktd(PolyList,Vars,ResVars[1])==1){
387: if(ResVars[0]==0){
388: RET=append(RET,wsort(ResVars[1],Vars,
389: ResVars[1],FLAG))$
390: return RET$
391: }
392: else{
1.11 kimura 393: RET=append(RET,[[Vars,vtol(ResVars[1])]])$
1.10 kimura 394: return RET$
395: }
1.4 kimura 396: }
397: else
1.8 kimura 398: return []$
1.4 kimura 399: }
400: else
1.8 kimura 401: return []$
1.4 kimura 402:
403: }
404:
1.10 kimura 405: def leastsq(NormMat,ExpMat,Vars,FLAG){
1.8 kimura 406:
407: RET=[]$
1.4 kimura 408:
409: ExpMatRowNum=size(ExpMat)[0]$
410: ExpMatColNum=size(ExpMat[0])[0]$
411:
1.8 kimura 412: if(NormMat==0){
413: NormMat=newmat(ExpMatColNum,ExpMatColNum)$
414:
415: for(I=0;I<ExpMatColNum;I++)
416: for(J=I;J<ExpMatColNum;J++)
417: for(K=0;K<ExpMatRowNum;K++)
418: NormMat[I][J]+=
419: ExpMat[K][I]*ExpMat[K][J]$
420: }
1.4 kimura 421:
1.8 kimura 422: BVec=newvect(ExpMatColNum)$
1.4 kimura 423:
424: for(I=0;I<ExpMatColNum;I++)
425: for(J=0;J<ExpMatRowNum;J++)
1.8 kimura 426: BVec[I]+=ExpMat[J][I]$
1.4 kimura 427:
428: SolveList=[]$
429: for(I=0;I<ExpMatColNum;I++){
430: TMP=0$
1.8 kimura 431: for(J=0;J<I;J++)
432: TMP+=NormMat[J][I]*Vars[J]$
433:
434: for(J=I;J<ExpMatColNum;J++)
1.4 kimura 435: TMP+=NormMat[I][J]*Vars[J]$
436:
1.8 kimura 437: TMP-=BVec[I]$
1.4 kimura 438: SolveList=cons(TMP,SolveList)$
439: }
440:
441: Rea=vars(SolveList)$
1.8 kimura 442:
443: VarsList=[]$
444: for(I=0;I<length(Vars);I++)
445: if(member(Vars[I],Rea))
446: VarsList=cons(Vars[I],VarsList)$
447:
448: Res=solve(SolveList,VarsList)$
449: Res=getgcd(Res,Rea)$
1.4 kimura 450:
451: if(nonposdegchk(Res)){
1.16 kimura 452:
1.10 kimura 453: TMP1=makeret(Res,Vars,1)$
1.16 kimura 454:
1.8 kimura 455: if(TMP1[0]==0){
1.14 kimura 456: TMP=roundret(TMP1[1])$
1.8 kimura 457: if(TMP!=[])
1.10 kimura 458: RET=append(RET,wsort(TMP1[1],Vars,TMP,FLAG))$
1.8 kimura 459:
1.10 kimura 460: RET=append(RET,wsort(TMP1[1],Vars,
461: map(drint,TMP1[1]*1.0),FLAG))$
1.8 kimura 462:
463: return RET$
464: }
465: else{
1.11 kimura 466: RET=append(RET,[[Vars,vtol(TMP1[1]*1.0)]])$
1.8 kimura 467: return RET$
468: }
469: }
470: else
471: return RET$
472:
473: }
474:
1.10 kimura 475: def weightr(ExpMat,Vars,PolyListNum,OneMat,FLAG){
1.8 kimura 476:
477: RET=[]$
478:
479: ExpMatRowNum=size(ExpMat)[0]$
480: ExpMatColNum=size(ExpMat[0])[0]$
481: ExtMatColNum=ExpMatColNum+PolyListNum$
482:
483: ExtVars=reverse(Vars)$
484: for(I=0;I<PolyListNum;I++)
485: ExtVars=cons(uc(),ExtVars)$
486:
487: ExtVars=reverse(ExtVars)$
488:
489: NormMat=newmat(ExpMatColNum,ExtMatColNum)$
490:
491: for(I=0;I<ExpMatColNum;I++)
492: for(J=I;J<ExpMatColNum;J++)
493: for(K=0;K<ExpMatRowNum;K++)
494: NormMat[I][J]+=
495: ExpMat[K][I]*
496: ExpMat[K][J]$
497:
498: for(I=0;I<ExpMatColNum;I++)
499: for(J=0;J<PolyListNum;J++)
500: for(K=OneMat[J];K<OneMat[J+1];K++)
501: NormMat[I][J+ExpMatColNum]-=
502: ExpMat[K][I]$
503:
504: WVect=newvect(PolyListNum)$
505: for(I=0;I<PolyListNum;I++)
506: WVect[I]=OneMat[I+1]-OneMat[I]$
507:
508: for(F=0;F<ExtMatColNum;F++){
509: SolveList=[]$
510: for(I=0;I<ExpMatColNum;I++){
511: if (F==I)
512: continue$
513:
514: TMP=0$
515:
516: for(J=0;J<I;J++)
517: if(J!=F)
518: TMP+=NormMat[J][I]*ExtVars[J]$
519:
520: for(J=I;J<ExtMatColNum;J++)
521: if(J!=F)
522: TMP+=NormMat[I][J]*ExtVars[J]$
523:
524: if(F<I)
525: TMP+=NormMat[F][I]$
526: else
527: TMP+=NormMat[I][F]$
528:
529: SolveList=cons(TMP,SolveList)$
530: }
531:
532: for(I=0;I<PolyListNum;I++){
533: if(F==(I+ExpMatColNum))
534: continue$
535:
536: TMP=0$
537: for(J=0;J<ExpMatColNum;J++)
538: if(J!=F)
539: TMP+=NormMat[J][I+ExpMatColNum]
540: *ExtVars[J]$
541:
542: TMP+=WVect[I]*ExtVars[I+ExpMatColNum]$
543:
544: if(F<ExpMatColNum)
545: TMP+=NormMat[F][I+ExpMatColNum]$
546:
547: SolveList=cons(TMP,SolveList)$
548: }
549:
550: Rea=vars(SolveList)$
551:
552: SolVars=[]$
553: for(I=0;I<ExtMatColNum;I++)
554: if(I!=F && member(ExtVars[I],Rea))
555: SolVars=cons(ExtVars[I],SolVars)$
556:
557: Res=solve(SolveList,SolVars)$
558: Res=cons([ExtVars[F],1],Res)$
559:
1.9 kimura 560: TMP=[]$
561: for(I=0;I<length(Rea);I++)
562: if(member(Rea[I],Vars))
563: TMP=cons(Rea[I],TMP)$
564:
1.15 kimura 565: if(member(ExtVars[F],Vars))
566: TMP=cons(ExtVars[F],TMP)$
567:
1.9 kimura 568: Res=getgcd(Res,TMP)$
1.8 kimura 569:
570: if(nonposdegchk(Res)){
571:
1.10 kimura 572: TMP1=makeret(Res,Vars,1)$
1.16 kimura 573:
1.8 kimura 574: if(TMP1[0]==0){
1.14 kimura 575: TMP=roundret(TMP1[1])$
1.8 kimura 576: if(TMP!=[])
1.10 kimura 577: RET=append(RET,wsort(TMP1[1],Vars,
578: TMP,FLAG))$
1.8 kimura 579:
1.10 kimura 580: RET=append(RET,wsort(TMP1[1],Vars,
581: map(drint,TMP1[1]*1.0),FLAG))$
1.8 kimura 582: }
583: else{
1.11 kimura 584: RET=append(RET,[[Vars,vtol(TMP1[1]*1.0)]])$
1.8 kimura 585: }
1.4 kimura 586: }
1.8 kimura 587:
1.4 kimura 588: }
589:
1.8 kimura 590: return [NormMat,RET]$
1.6 kimura 591: }
592:
1.10 kimura 593: def weight(PolyList,Vars,FLAG1,FLAG2){
1.6 kimura 594:
595: Vars0=vars(PolyList)$
596: Vars1=[]$
597: for(I=0;I<length(Vars);I++)
598: if(member(Vars[I],Vars0))
599: Vars1=cons(Vars[I],Vars1)$
600:
601: Vars=reverse(Vars1)$
602:
603: RET=[]$
604:
1.10 kimura 605: TMP=qcheck(PolyList,Vars,FLAG2)$
1.6 kimura 606:
1.8 kimura 607: if(TMP!=[]){
608: RET=append(RET,TMP)$
609: return cons(1,RET)$
1.6 kimura 610: }
1.4 kimura 611:
1.6 kimura 612: dp_ord(2)$
1.4 kimura 613:
1.6 kimura 614: PolyListNum=length(PolyList)$
1.4 kimura 615:
1.10 kimura 616: if(FLAG1){
1.9 kimura 617:
618: OneMat=newvect(PolyListNum+1,[0])$
619: ExpMat=[]$
620: for(I=0;I<PolyListNum;I++){
621: for(Poly=dp_ptod(PolyList[I],Vars);
622: Poly!=0;Poly=dp_rest(Poly)){
623: ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
624: }
625: OneMat[I+1]=length(ExpMat)$
1.8 kimura 626: }
1.4 kimura 627:
1.9 kimura 628: ExpMat=reverse(ExpMat)$
629: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 630:
1.10 kimura 631: TMP=weightr(ExpMat,Vars,PolyListNum,OneMat,FLAG2)$
1.8 kimura 632: RET=append(RET,TMP[1])$
1.10 kimura 633: RET=append(RET,leastsq(TMP[0],ExpMat,Vars,FLAG2))$
1.8 kimura 634: }
1.9 kimura 635: else{
636: ExpMat=[]$
637: for(I=0;I<PolyListNum;I++){
638: for(Poly=dp_ptod(PolyList[I],Vars);
639: Poly!=0;Poly=dp_rest(Poly)){
640: if(nonzerovec(TMP=dp_etov(dp_ht(Poly))))
641: ExpMat=cons(TMP,ExpMat)$
642: }
643: }
1.8 kimura 644:
1.9 kimura 645: ExpMat=reverse(ExpMat)$
646: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 647:
1.10 kimura 648: RET=append(RET,leastsq(0,ExpMat,Vars,FLAG2))$
1.9 kimura 649: }
1.4 kimura 650:
1.6 kimura 651: ExpMat=qsort(ExpMat,junban)$
1.8 kimura 652:
1.6 kimura 653: ExpMat2=[]$
654: for(I=0;I<size(ExpMat)[0];I++)
655: if(car(ExpMat2)!=ExpMat[I])
656: ExpMat2=cons(ExpMat[I],ExpMat2)$
1.4 kimura 657:
1.6 kimura 658: if(size(ExpMat)[0]!=length(ExpMat2)){
659: ExpMat=newvect(length(ExpMat2),ExpMat2)$
1.10 kimura 660: RET=append(RET,leastsq(0,ExpMat,Vars,FLAG2))$
1.4 kimura 661: }
662:
1.8 kimura 663: RET=derase(RET)$
664: return cons(0,RET)$
1.4 kimura 665: }
666:
667: end$