![[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++)
166: if(C[J]==A[I][0])
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.8 kimura 206: ResVec=newvect(ResNum)$
207: for(M=0,I=0;I<ResNum;I++){
208: if(member(Res[I][0],Vars)){
209: ResVec[I]=Res[I][1]$
1.4 kimura 210:
1.10 kimura 211: if(FLAG && type(ResVec[I])==1){
1.4 kimura 212: if(M==0)
1.8 kimura 213: M=ResVec[I]$
1.4 kimura 214: else
1.8 kimura 215: if(ResVec[I]<M)
216: M=ResVec[I]$
1.4 kimura 217: }
218: }
1.8 kimura 219: }
220:
221: if(M!=0)
222: ResVec=ResVec/M;
1.4 kimura 223:
1.8 kimura 224: RET=newvect(VarsNum,Vars)$
1.4 kimura 225:
1.8 kimura 226: for(I=0;I<ResNum;I++){
227: for(J=0;J<VarsNum;J++)
228: if(Vars[J]==Res[I][0])
229: break$
1.4 kimura 230:
1.8 kimura 231: if(J<VarsNum)
232: RET[J]=ResVec[I]$
233: }
234:
1.10 kimura 235:
236: for(J=0;J<length(Vars);J++)
237: RET=map(subst,RET,Vars[J],
238: strtov(rtostr(Vars[J])+"_deg"))$
239:
1.4 kimura 240: for(I=0;I<VarsNum;I++)
1.8 kimura 241: if(type(RET[I])!=1)
242: return [1,RET]$
1.4 kimura 243:
1.8 kimura 244: return [0,RET]$
1.4 kimura 245: }
246:
247: def roundret(V){
248:
1.8 kimura 249: VN=size(V)[0]$
1.4 kimura 250:
1.8 kimura 251: RET0=V$
1.13 kimura 252: for(I=1;I<1000;I++){
1.4 kimura 253: RET1=I*RET0$
254: for(J=0;J<VN;J++){
255: X=drint(RET1[J])$
256: if(dabs(X-RET1[J])<0.2)
257: RET1[J]=X$
258: else
259: break$
260: }
261: if(J==VN)
262: break$
263: }
264:
265: if(I==1000)
266: return []$
267: else
268: return RET1$
269: }
270:
271: def chkou(L,ExpMat,CHAGORD){
272:
1.8 kimura 273: for(P=1,I=0;I<L;I++){
1.4 kimura 274: Q=ExpMat[L][CHAGORD[I]]$
275: for(J=0;J<size(ExpMat[0])[0];J++){
276: ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
277: *ExpMat[L][CHAGORD[J]]-
278: Q*ExpMat[I][CHAGORD[J]])/P)$
279: }
280:
281: P=ExpMat[I][CHAGORD[I]]$
282: }
283:
284: for(J=0;J<size(ExpMat[0])[0];J++)
285: if(ExpMat[L][CHAGORD[J]]!=0)
286: break$
287:
288: if(J==size(ExpMat[0])[0])
289: return L$
290: else{
291: TMP=CHAGORD[L]$
292: CHAGORD[L]=CHAGORD[J]$
293: CHAGORD[J]=TMP$
294: return (L+1)$
295: }
296: }
297:
1.8 kimura 298: def qcheckmain(PolyList,Vars){
1.4 kimura 299:
300: RET=[]$
301: PolyListNum=length(PolyList)$
302: VarsNum=length(Vars)$
303:
304: ExpMat=newvect(VarsNum)$
305: CHAGORD=newvect(VarsNum)$
306: for(I=0;I<VarsNum;I++)
307: CHAGORD[I]=I$
308:
309: L=0$
310: for(I=0;I<PolyListNum;I++){
311: Poly=dp_ptod(PolyList[I],Vars)$
312: BASE0=dp_etov(dp_ht(Poly))$
313: Poly=dp_rest(Poly)$
314: for(;Poly!=0;Poly=dp_rest(Poly)){
315: ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
316: L=chkou(L,ExpMat,CHAGORD)$
1.8 kimura 317: if(L==VarsNum-1)
318: return [L,CHAGORD,ExpMat]$
1.4 kimura 319: }
320: }
321:
1.8 kimura 322: return [L,CHAGORD,ExpMat]$
1.4 kimura 323: }
324:
325: def inner(A,B){
326:
327: SUM=0$
328: for(I=0;I<size(A)[0];I++)
329: SUM+=A[I]*B[I]$
330:
331: return SUM$
332: }
333:
334: def checktd(PolyList,Vars,ResVars){
335:
336: PolyListNum=length(PolyList)$
337: VarsNum=length(Vars)$
338:
339: L=0$
340: for(I=0;I<PolyListNum;I++){
341: Poly=dp_ptod(PolyList[I],Vars)$
342: J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
343: Poly=dp_rest(Poly)$
344: for(;Poly!=0;Poly=dp_rest(Poly))
345: if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
346: return 0$
347: }
348:
349: return 1$
350: }
351:
1.10 kimura 352: def qcheck(PolyList,Vars,FLAG){
1.4 kimura 353:
1.10 kimura 354: RET=[]$
1.8 kimura 355: Res=qcheckmain(PolyList,Vars)$
1.4 kimura 356: VarsNum=length(Vars)$
357:
358: IndNum=Res[0]$
359: CHAGORD=Res[1]$
360: ExpMat=Res[2]$
361:
362: SolveList=[]$
363: for(I=0;I<IndNum;I++){
364: TMP=0$
365: for(J=0;J<VarsNum;J++)
366: TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$
367:
368: SolveList=cons(TMP,SolveList)$
369: }
370:
1.8 kimura 371: Rea=vars(SolveList)$
372:
1.4 kimura 373: VarsList=[]$
374: for(I=0;I<VarsNum;I++)
1.8 kimura 375: if(member(Vars[CHAGORD[I]],Rea))
376: VarsList=cons(Vars[CHAGORD[I]],VarsList)$
1.4 kimura 377:
378: Res=solve(reverse(SolveList),reverse(VarsList))$
1.8 kimura 379: Res=getgcd(Res,Rea)$
1.4 kimura 380:
381: if(nonposdegchk(Res)){
382:
1.10 kimura 383: ResVars=makeret(Res,Vars,0)$
1.4 kimura 384:
1.10 kimura 385: if(checktd(PolyList,Vars,ResVars[1])==1){
386: if(ResVars[0]==0){
387: RET=append(RET,wsort(ResVars[1],Vars,
388: ResVars[1],FLAG))$
389: return RET$
390: }
391: else{
1.11 kimura 392: RET=append(RET,[[Vars,vtol(ResVars[1])]])$
1.10 kimura 393: return RET$
394: }
1.4 kimura 395: }
396: else
1.8 kimura 397: return []$
1.4 kimura 398: }
399: else
1.8 kimura 400: return []$
1.4 kimura 401:
402: }
403:
1.10 kimura 404: def leastsq(NormMat,ExpMat,Vars,FLAG){
1.8 kimura 405:
406: RET=[]$
1.4 kimura 407:
408: ExpMatRowNum=size(ExpMat)[0]$
409: ExpMatColNum=size(ExpMat[0])[0]$
410:
1.8 kimura 411: if(NormMat==0){
412: NormMat=newmat(ExpMatColNum,ExpMatColNum)$
413:
414: for(I=0;I<ExpMatColNum;I++)
415: for(J=I;J<ExpMatColNum;J++)
416: for(K=0;K<ExpMatRowNum;K++)
417: NormMat[I][J]+=
418: ExpMat[K][I]*ExpMat[K][J]$
419: }
1.4 kimura 420:
1.8 kimura 421: BVec=newvect(ExpMatColNum)$
1.4 kimura 422:
423: for(I=0;I<ExpMatColNum;I++)
424: for(J=0;J<ExpMatRowNum;J++)
1.8 kimura 425: BVec[I]+=ExpMat[J][I]$
1.4 kimura 426:
427: SolveList=[]$
428: for(I=0;I<ExpMatColNum;I++){
429: TMP=0$
1.8 kimura 430: for(J=0;J<I;J++)
431: TMP+=NormMat[J][I]*Vars[J]$
432:
433: for(J=I;J<ExpMatColNum;J++)
1.4 kimura 434: TMP+=NormMat[I][J]*Vars[J]$
435:
1.8 kimura 436: TMP-=BVec[I]$
1.4 kimura 437: SolveList=cons(TMP,SolveList)$
438: }
439:
440: Rea=vars(SolveList)$
1.8 kimura 441:
442: VarsList=[]$
443: for(I=0;I<length(Vars);I++)
444: if(member(Vars[I],Rea))
445: VarsList=cons(Vars[I],VarsList)$
446:
447: Res=solve(SolveList,VarsList)$
448: Res=getgcd(Res,Rea)$
1.4 kimura 449:
450: if(nonposdegchk(Res)){
1.10 kimura 451: TMP1=makeret(Res,Vars,1)$
1.8 kimura 452: if(TMP1[0]==0){
1.14 kimura 453: TMP=roundret(TMP1[1])$
1.8 kimura 454: if(TMP!=[])
1.10 kimura 455: RET=append(RET,wsort(TMP1[1],Vars,TMP,FLAG))$
1.8 kimura 456:
1.10 kimura 457: RET=append(RET,wsort(TMP1[1],Vars,
458: map(drint,TMP1[1]*1.0),FLAG))$
1.8 kimura 459:
460: return RET$
461: }
462: else{
1.11 kimura 463: RET=append(RET,[[Vars,vtol(TMP1[1]*1.0)]])$
1.8 kimura 464: return RET$
465: }
466: }
467: else
468: return RET$
469:
470: }
471:
1.10 kimura 472: def weightr(ExpMat,Vars,PolyListNum,OneMat,FLAG){
1.8 kimura 473:
474: RET=[]$
475:
476: ExpMatRowNum=size(ExpMat)[0]$
477: ExpMatColNum=size(ExpMat[0])[0]$
478: ExtMatColNum=ExpMatColNum+PolyListNum$
479:
480: ExtVars=reverse(Vars)$
481: for(I=0;I<PolyListNum;I++)
482: ExtVars=cons(uc(),ExtVars)$
483:
484: ExtVars=reverse(ExtVars)$
485:
486: NormMat=newmat(ExpMatColNum,ExtMatColNum)$
487:
488: for(I=0;I<ExpMatColNum;I++)
489: for(J=I;J<ExpMatColNum;J++)
490: for(K=0;K<ExpMatRowNum;K++)
491: NormMat[I][J]+=
492: ExpMat[K][I]*
493: ExpMat[K][J]$
494:
495: for(I=0;I<ExpMatColNum;I++)
496: for(J=0;J<PolyListNum;J++)
497: for(K=OneMat[J];K<OneMat[J+1];K++)
498: NormMat[I][J+ExpMatColNum]-=
499: ExpMat[K][I]$
500:
501: WVect=newvect(PolyListNum)$
502: for(I=0;I<PolyListNum;I++)
503: WVect[I]=OneMat[I+1]-OneMat[I]$
504:
505: for(F=0;F<ExtMatColNum;F++){
506: SolveList=[]$
507: for(I=0;I<ExpMatColNum;I++){
508: if (F==I)
509: continue$
510:
511: TMP=0$
512:
513: for(J=0;J<I;J++)
514: if(J!=F)
515: TMP+=NormMat[J][I]*ExtVars[J]$
516:
517: for(J=I;J<ExtMatColNum;J++)
518: if(J!=F)
519: TMP+=NormMat[I][J]*ExtVars[J]$
520:
521: if(F<I)
522: TMP+=NormMat[F][I]$
523: else
524: TMP+=NormMat[I][F]$
525:
526: SolveList=cons(TMP,SolveList)$
527: }
528:
529: for(I=0;I<PolyListNum;I++){
530: if(F==(I+ExpMatColNum))
531: continue$
532:
533: TMP=0$
534: for(J=0;J<ExpMatColNum;J++)
535: if(J!=F)
536: TMP+=NormMat[J][I+ExpMatColNum]
537: *ExtVars[J]$
538:
539: TMP+=WVect[I]*ExtVars[I+ExpMatColNum]$
540:
541: if(F<ExpMatColNum)
542: TMP+=NormMat[F][I+ExpMatColNum]$
543:
544: SolveList=cons(TMP,SolveList)$
545: }
546:
547: Rea=vars(SolveList)$
548:
549: SolVars=[]$
550: for(I=0;I<ExtMatColNum;I++)
551: if(I!=F && member(ExtVars[I],Rea))
552: SolVars=cons(ExtVars[I],SolVars)$
553:
554: Res=solve(SolveList,SolVars)$
555: Res=cons([ExtVars[F],1],Res)$
556:
1.9 kimura 557: TMP=[]$
558: for(I=0;I<length(Rea);I++)
559: if(member(Rea[I],Vars))
560: TMP=cons(Rea[I],TMP)$
561:
1.15 ! kimura 562: if(member(ExtVars[F],Vars))
! 563: TMP=cons(ExtVars[F],TMP)$
! 564:
1.9 kimura 565: Res=getgcd(Res,TMP)$
1.8 kimura 566:
567: if(nonposdegchk(Res)){
568:
1.10 kimura 569: TMP1=makeret(Res,Vars,1)$
1.8 kimura 570: if(TMP1[0]==0){
1.14 kimura 571: TMP=roundret(TMP1[1])$
1.8 kimura 572: if(TMP!=[])
1.10 kimura 573: RET=append(RET,wsort(TMP1[1],Vars,
574: TMP,FLAG))$
1.8 kimura 575:
1.10 kimura 576: RET=append(RET,wsort(TMP1[1],Vars,
577: map(drint,TMP1[1]*1.0),FLAG))$
1.8 kimura 578: }
579: else{
1.11 kimura 580: RET=append(RET,[[Vars,vtol(TMP1[1]*1.0)]])$
1.8 kimura 581: }
1.4 kimura 582: }
1.8 kimura 583:
1.4 kimura 584: }
585:
1.8 kimura 586: return [NormMat,RET]$
1.6 kimura 587: }
588:
1.10 kimura 589: def weight(PolyList,Vars,FLAG1,FLAG2){
1.6 kimura 590:
591: Vars0=vars(PolyList)$
592: Vars1=[]$
593: for(I=0;I<length(Vars);I++)
594: if(member(Vars[I],Vars0))
595: Vars1=cons(Vars[I],Vars1)$
596:
597: Vars=reverse(Vars1)$
598:
599: RET=[]$
600:
1.10 kimura 601: TMP=qcheck(PolyList,Vars,FLAG2)$
1.6 kimura 602:
1.8 kimura 603: if(TMP!=[]){
604: RET=append(RET,TMP)$
605: return cons(1,RET)$
1.6 kimura 606: }
1.4 kimura 607:
1.6 kimura 608: dp_ord(2)$
1.4 kimura 609:
1.6 kimura 610: PolyListNum=length(PolyList)$
1.4 kimura 611:
1.10 kimura 612: if(FLAG1){
1.9 kimura 613:
614: OneMat=newvect(PolyListNum+1,[0])$
615: ExpMat=[]$
616: for(I=0;I<PolyListNum;I++){
617: for(Poly=dp_ptod(PolyList[I],Vars);
618: Poly!=0;Poly=dp_rest(Poly)){
619: ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
620: }
621: OneMat[I+1]=length(ExpMat)$
1.8 kimura 622: }
1.4 kimura 623:
1.9 kimura 624: ExpMat=reverse(ExpMat)$
625: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 626:
1.10 kimura 627: TMP=weightr(ExpMat,Vars,PolyListNum,OneMat,FLAG2)$
1.8 kimura 628: RET=append(RET,TMP[1])$
1.10 kimura 629: RET=append(RET,leastsq(TMP[0],ExpMat,Vars,FLAG2))$
1.8 kimura 630: }
1.9 kimura 631: else{
632: ExpMat=[]$
633: for(I=0;I<PolyListNum;I++){
634: for(Poly=dp_ptod(PolyList[I],Vars);
635: Poly!=0;Poly=dp_rest(Poly)){
636: if(nonzerovec(TMP=dp_etov(dp_ht(Poly))))
637: ExpMat=cons(TMP,ExpMat)$
638: }
639: }
1.8 kimura 640:
1.9 kimura 641: ExpMat=reverse(ExpMat)$
642: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 643:
1.10 kimura 644: RET=append(RET,leastsq(0,ExpMat,Vars,FLAG2))$
1.9 kimura 645: }
1.4 kimura 646:
1.6 kimura 647: ExpMat=qsort(ExpMat,junban)$
1.8 kimura 648:
1.6 kimura 649: ExpMat2=[]$
650: for(I=0;I<size(ExpMat)[0];I++)
651: if(car(ExpMat2)!=ExpMat[I])
652: ExpMat2=cons(ExpMat[I],ExpMat2)$
1.4 kimura 653:
1.6 kimura 654: if(size(ExpMat)[0]!=length(ExpMat2)){
655: ExpMat=newvect(length(ExpMat2),ExpMat2)$
1.10 kimura 656: RET=append(RET,leastsq(0,ExpMat,Vars,FLAG2))$
1.4 kimura 657: }
658:
1.8 kimura 659: RET=derase(RET)$
660: return cons(0,RET)$
1.4 kimura 661: }
662:
663: end$