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