![[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.4 kimura 252: for(I=1;I<1000;I++){
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{
! 392: RET=append(RET,[[Vars,ResVars[1]]])$
! 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){
453: TMP=roundret(TMP1[1]*1.0)$
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.10 ! kimura 463: RET=append(RET,[[Vars,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:
562: TMP=cons(ExtVars[F],TMP)$
563: Res=getgcd(Res,TMP)$
1.8 kimura 564:
565: if(nonposdegchk(Res)){
566:
1.10 ! kimura 567: TMP1=makeret(Res,Vars,1)$
1.8 kimura 568: if(TMP1[0]==0){
569: TMP=roundret(TMP1[1]*1.0)$
570: if(TMP!=[])
1.10 ! kimura 571: RET=append(RET,wsort(TMP1[1],Vars,
! 572: TMP,FLAG))$
1.8 kimura 573:
1.10 ! kimura 574: RET=append(RET,wsort(TMP1[1],Vars,
! 575: map(drint,TMP1[1]*1.0),FLAG))$
1.8 kimura 576: }
577: else{
1.10 ! kimura 578: RET=append(RET,[[Vars,TMP1[1]*1.0]])$
1.8 kimura 579: }
1.4 kimura 580: }
1.8 kimura 581:
1.4 kimura 582: }
583:
1.8 kimura 584: return [NormMat,RET]$
1.6 kimura 585: }
586:
1.10 ! kimura 587: def weight(PolyList,Vars,FLAG1,FLAG2){
1.6 kimura 588:
589: Vars0=vars(PolyList)$
590: Vars1=[]$
591: for(I=0;I<length(Vars);I++)
592: if(member(Vars[I],Vars0))
593: Vars1=cons(Vars[I],Vars1)$
594:
595: Vars=reverse(Vars1)$
596:
597: RET=[]$
598:
1.10 ! kimura 599: TMP=qcheck(PolyList,Vars,FLAG2)$
1.6 kimura 600:
1.8 kimura 601: if(TMP!=[]){
602: RET=append(RET,TMP)$
603: return cons(1,RET)$
1.6 kimura 604: }
1.4 kimura 605:
1.6 kimura 606: dp_ord(2)$
1.4 kimura 607:
1.6 kimura 608: PolyListNum=length(PolyList)$
1.4 kimura 609:
1.10 ! kimura 610: if(FLAG1){
1.9 kimura 611:
612: OneMat=newvect(PolyListNum+1,[0])$
613: ExpMat=[]$
614: for(I=0;I<PolyListNum;I++){
615: for(Poly=dp_ptod(PolyList[I],Vars);
616: Poly!=0;Poly=dp_rest(Poly)){
617: ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
618: }
619: OneMat[I+1]=length(ExpMat)$
1.8 kimura 620: }
1.4 kimura 621:
1.9 kimura 622: ExpMat=reverse(ExpMat)$
623: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 624:
1.10 ! kimura 625: TMP=weightr(ExpMat,Vars,PolyListNum,OneMat,FLAG2)$
1.8 kimura 626: RET=append(RET,TMP[1])$
1.10 ! kimura 627: RET=append(RET,leastsq(TMP[0],ExpMat,Vars,FLAG2))$
1.8 kimura 628: }
1.9 kimura 629: else{
630: ExpMat=[]$
631: for(I=0;I<PolyListNum;I++){
632: for(Poly=dp_ptod(PolyList[I],Vars);
633: Poly!=0;Poly=dp_rest(Poly)){
634: if(nonzerovec(TMP=dp_etov(dp_ht(Poly))))
635: ExpMat=cons(TMP,ExpMat)$
636: }
637: }
1.8 kimura 638:
1.9 kimura 639: ExpMat=reverse(ExpMat)$
640: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 641:
1.10 ! kimura 642: RET=append(RET,leastsq(0,ExpMat,Vars,FLAG2))$
1.9 kimura 643: }
1.4 kimura 644:
1.6 kimura 645: ExpMat=qsort(ExpMat,junban)$
1.8 kimura 646:
1.6 kimura 647: ExpMat2=[]$
648: for(I=0;I<size(ExpMat)[0];I++)
649: if(car(ExpMat2)!=ExpMat[I])
650: ExpMat2=cons(ExpMat[I],ExpMat2)$
1.4 kimura 651:
1.6 kimura 652: if(size(ExpMat)[0]!=length(ExpMat2)){
653: ExpMat=newvect(length(ExpMat2),ExpMat2)$
1.10 ! kimura 654: RET=append(RET,leastsq(0,ExpMat,Vars,FLAG2))$
1.4 kimura 655: }
656:
1.8 kimura 657: RET=derase(RET)$
658: return cons(0,RET)$
1.4 kimura 659: }
660:
661: end$