![[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.2 kimura 1: load("solve")$
1.3 ! kimura 2: load("gr")$
1.2 kimura 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:
17: for(I=0;I<length(Res);I++){
18:
19: for(J=0;J<size(ResVars)[0];J++)
20: if(Res[I][0]==ResVars[J])
21: break$
22:
23: ResVars[J]=Res[I][1]$
24: }
25:
26: return(ResVars)$
27: }
28:
1.3 ! kimura 29: def makeret1(Res,Vars){
1.2 kimura 30:
31: VarsNum=length(Vars)$
32:
1.3 ! kimura 33: ResVec=newvect(VarsNum,Vars)$
1.2 kimura 34:
1.3 ! kimura 35: for(I=0,M=0;I<length(Res);I++){
1.2 kimura 36:
1.3 ! kimura 37: for(J=0;J<VarsNum;J++)
! 38: if(Res[I][0]==Vars[J])
1.2 kimura 39: break$
40:
41: if(J<VarsNum){
1.3 ! kimura 42: ResVec[J]=Res[I][1]$
1.2 kimura 43:
1.3 ! kimura 44: if(type(ResVec[J])==1){
! 45: if(M==0)
! 46: M=ResVec[J]$
! 47: else
! 48: if(ResVec[J]<M)
! 49: M=ResVec[J]$
1.2 kimura 50: }
1.3 ! kimura 51: }
! 52:
! 53: }
1.2 kimura 54:
1.3 ! kimura 55: for(F=0,I=0;I<VarsNum;I++)
! 56: if(type(ResVec[I])!=1){
! 57: F=1$
! 58: break$
1.2 kimura 59: }
60:
61: if(F==0)
1.3 ! kimura 62: for(I=0;I<VarsNum;I++)
! 63: ResVec[I]=ResVec[I]/M*1.0$
1.2 kimura 64:
65: for(I=0;I<VarsNum;I++)
66: for(J=0;J<length(Vars);J++)
1.3 ! kimura 67: ResVec[I]=subst(ResVec[I],Vars[J],
1.2 kimura 68: strtov(rtostr(Vars[J])+"_deg"))$
69:
1.3 ! kimura 70: ResVec=cons(F,vtol(ResVec))$
! 71: return ResVec$
! 72: }
1.2 kimura 73:
1.3 ! kimura 74: def junban1(A,B){
! 75: return (nmono(A)<nmono(B) ? -1:(nmono(A)>nmono(B) ? 1:0))$
1.2 kimura 76: }
77:
1.3 ! kimura 78: def junban2(A,B){
1.2 kimura 79:
80: for(I=0;I<size(A)[0];I++){
81: if(A[I]<B[I])
82: return 1$
83:
84: if(A[I]>B[I])
85: return -1$
86: }
87:
88: return 0$
89: }
90:
1.3 ! kimura 91: def roundret(V){
! 92:
! 93: VN=length(V)$
! 94: RET0=newvect(VN,V)$
! 95:
! 96: for(I=1;I<1000;I++){
! 97: RET1=I*RET0$
! 98: for(J=0;J<VN;J++){
! 99: X=drint(RET1[J])$
! 100: if(dabs(X-RET1[J])<0.2)
! 101: RET1[J]=X$
! 102: else
! 103: break$
! 104: }
! 105: if(J==VN)
! 106: break$
! 107: }
! 108:
! 109: if(I==1000)
! 110: return []$
! 111: else
! 112: return RET1$
! 113: }
! 114:
! 115: def chkou(L,ExpMat,CHAGORD){
! 116:
! 117: P=1$
! 118: F=ExpMat[L]$
! 119:
! 120: for(I=0;I<L;I++){
! 121: Q=ExpMat[L][CHAGORD[I]]$
! 122: for(J=0;J<size(ExpMat[0])[0];J++){
! 123: ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
! 124: *ExpMat[L][CHAGORD[J]]-
! 125: Q*ExpMat[I][CHAGORD[J]])/P)$
! 126: }
! 127:
! 128: P=ExpMat[I][CHAGORD[I]]$
! 129: }
! 130:
! 131: for(J=0;J<size(ExpMat[0])[0];J++)
! 132: if(ExpMat[L][CHAGORD[J]]!=0)
! 133: break$
! 134:
! 135: if(J==size(ExpMat[0])[0])
! 136: return L$
! 137: else{
! 138: TMP=CHAGORD[L]$
! 139: CHAGORD[L]=CHAGORD[J]$
! 140: CHAGORD[J]=TMP$
! 141: return (L+1)$
! 142: }
! 143: }
! 144:
! 145: def qcheck0(PolyList,Vars){
! 146:
! 147: RET=[]$
! 148: PolyListNum=length(PolyList)$
! 149: VarsNum=length(Vars)$
! 150:
! 151: ExpMat=newvect(VarsNum)$
! 152: CHAGORD=newvect(VarsNum)$
! 153: for(I=0;I<VarsNum;I++)
! 154: CHAGORD[I]=I$
! 155:
! 156: L=0$
! 157: for(I=0;I<PolyListNum;I++){
! 158: Poly=dp_ptod(PolyList[I],Vars)$
! 159: BASE0=dp_etov(dp_ht(Poly))$
! 160: Poly=dp_rest(Poly)$
! 161: for(;Poly!=0;Poly=dp_rest(Poly)){
! 162: ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
! 163: L=chkou(L,ExpMat,CHAGORD)$
! 164: if(L==VarsNum-1){
! 165: RET=cons(ExpMat,RET)$
! 166: RET=cons(CHAGORD,RET)$
! 167: RET=cons(L,RET)$
! 168: return RET$
! 169: }
! 170: }
! 171: }
! 172:
! 173: RET=cons(ExpMat,RET)$
! 174: RET=cons(CHAGORD,RET)$
! 175: RET=cons(L,RET)$
! 176: return RET$
! 177: }
! 178:
! 179: def inner(A,B){
! 180:
! 181: SUM=0$
! 182: for(I=0;I<size(A)[0];I++)
! 183: SUM+=A[I]*B[I]$
! 184:
! 185: return SUM$
! 186: }
! 187:
! 188: def checktd(PolyList,Vars,ResVars){
! 189:
! 190: PolyListNum=length(PolyList)$
! 191: VarsNum=length(Vars)$
! 192:
! 193: L=0$
! 194: for(I=0;I<PolyListNum;I++){
! 195: Poly=dp_ptod(PolyList[I],Vars)$
! 196: J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
! 197: Poly=dp_rest(Poly)$
! 198: for(;Poly!=0;Poly=dp_rest(Poly))
! 199: if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
! 200: return 0$
! 201: }
! 202:
! 203: return 1$
! 204: }
! 205:
! 206: def getgcd(A,B){
! 207:
! 208: VarsNumA=length(A)$
! 209: VarsNumB=length(B)$
! 210:
! 211: C=newvect(VarsNumB,B)$
! 212:
! 213: for(I=0;I<VarsNumA;I++){
! 214:
! 215: for(J=0;J<VarsNumB;J++)
! 216: if(C[J]==A[I][0])
! 217: break$
! 218:
! 219: C[J]=A[I][1]$
! 220: }
! 221:
! 222: D=0$
! 223: for(I=0;I<VarsNumB;I++)
! 224: D=gcd(D,C[I])$
! 225:
! 226: if(D!=0){
! 227:
! 228: for(I=0;I<VarsNumB;I++)
! 229: C[I]=red(C[I]/D)$
! 230:
! 231: }
! 232:
! 233: for(L=1,D=0,I=0;I<VarsNumB;I++){
! 234:
! 235: if(type(C[I])==1){
! 236: L=ilcm(L,dn(C[I]))$
! 237: D=igcd(D,nm(C[I]))$
! 238: }
! 239: else
! 240: break$
! 241:
! 242: }
! 243:
! 244: if(I==VarsNumB)
! 245: for(I=0;I<VarsNumB;I++)
! 246: C[I]=C[I]*L/D$
! 247:
! 248: RET=newvect(VarsNumB)$
! 249: for(I=0;I<VarsNumB;I++){
! 250: RET[I]=newvect(2)$
! 251: RET[I][0]=B[I]$
! 252: RET[I][1]=C[I]$
! 253: }
! 254:
! 255: return vtol(map(vtol,RET))$
! 256: }
! 257:
! 258: def qcheck(PolyList,Vars){
! 259:
! 260: RET=[]$
! 261: Res=qcheck0(PolyList,Vars)$
! 262: VarsNum=length(Vars)$
! 263:
! 264: IndNum=Res[0]$
! 265: CHAGORD=Res[1]$
! 266: ExpMat=Res[2]$
! 267:
! 268: SolveList=[]$
! 269: for(I=0;I<IndNum;I++){
! 270: TMP=0$
! 271: for(J=0;J<VarsNum;J++)
! 272: TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$
! 273:
! 274: SolveList=cons(TMP,SolveList)$
! 275: }
! 276:
! 277: VarsList=[]$
! 278: for(I=0;I<VarsNum;I++)
! 279: VarsList=cons(Vars[CHAGORD[I]],VarsList)$
! 280:
! 281: Rea=vars(SolveList)$
! 282: Res=solve(reverse(SolveList),reverse(VarsList))$
! 283:
! 284: if(nonposdegchk(Res)){
! 285:
! 286: Res=getgcd(Res,Rea)$
! 287: ResVars=resvars(Res,Vars)$
! 288:
! 289: if(checktd(PolyList,Vars,ResVars)==1){
! 290:
! 291: for(J=0;J<length(Vars);J++)
! 292: ResVars=map(subst,ResVars,Vars[J],
! 293: strtov(rtostr(Vars[J])+"_deg"))$
! 294:
! 295: RET=cons([vtol(ResVars),ResVars,[]],RET)$
! 296: return cons(1,RET)$
! 297: }
! 298: else
! 299: return cons(0,RET)$
! 300: }
! 301: else
! 302: return cons(0,RET)$
! 303:
! 304: }
! 305:
1.2 kimura 306: def weight(PolyList,Vars){
307:
1.3 ! kimura 308: Vars0=vars(PolyList)$
! 309: Vars1=[]$
! 310: for(I=0;I<length(Vars);I++)
! 311: if(member(Vars[I],Vars0))
! 312: Vars1=cons(Vars[I],Vars1)$
! 313:
! 314: Vars=reverse(Vars1)$
! 315:
! 316: RET=[]$
! 317:
! 318: TMP=qcheck(PolyList,Vars)$
! 319:
! 320: if(car(TMP)==1){
! 321: RET=cdr(TMP)$
! 322: RET=cons(Vars,RET)$
! 323: RET=cons(1,RET)$
! 324: return RET$
! 325: }
! 326:
1.2 kimura 327: dp_ord(2)$
328:
329: PolyListNum=length(PolyList)$
1.3 ! kimura 330: VPolyList=qsort(newvect(PolyListNum,PolyList),junban1)$
! 331: VPolyList=vtol(VPolyList)$
1.2 kimura 332:
333: ExpMat=[]$
334: for(I=0;I<PolyListNum;I++)
1.3 ! kimura 335: for(Poly=dp_ptod(VPolyList[I],Vars);Poly!=0;Poly=dp_rest(Poly))
1.2 kimura 336: ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
337:
338: ExpMat=reverse(ExpMat)$
339: ExpMat=newvect(length(ExpMat),ExpMat)$
340:
1.3 ! kimura 341:
! 342: /* first */
! 343:
1.2 kimura 344: ExpMatRowNum=size(ExpMat)[0]$
345: ExpMatColNum=size(ExpMat[0])[0]$
346: ExtMatColNum=ExpMatColNum+PolyListNum$
347:
348: OneMat=newvect(PolyListNum+1,[0])$
1.3 ! kimura 349: for(I=0,SUM=0;I<PolyListNum;I++){
! 350: SUM+=nmono(VPolyList[I])$
1.2 kimura 351: OneMat[I+1]=SUM$
352: }
353:
354: RevOneMat=newvect(ExpMatRowNum)$
355: for(I=0;I<PolyListNum;I++)
356: for(J=OneMat[I];J<OneMat[I+1];J++)
357: RevOneMat[J]=I$
358:
359: NormMat=newmat(ExpMatColNum,ExtMatColNum)$
360:
361: for(I=0;I<ExpMatColNum;I++)
362: for(J=0;J<ExpMatColNum;J++)
363: for(K=0;K<ExpMatRowNum;K++)
364: NormMat[I][J]+=ExpMat[K][I]*ExpMat[K][J]$
365:
366: for(I=0;I<ExpMatColNum;I++)
367: for(J=0;J<PolyListNum-1;J++)
368: for(K=OneMat[J];K<OneMat[J+1];K++)
369: NormMat[I][J+ExpMatColNum]-=ExpMat[K][I]$
370:
371: for(I=0;I<ExpMatColNum;I++)
372: for(J=OneMat[PolyListNum-1];J<OneMat[PolyListNum];J++)
373: NormMat[I][ExtMatColNum-1]+=ExpMat[J][I]$
374:
375: NormMat2=newmat(PolyListNum-1,ExpMatColNum+1)$
376:
377: for(I=0;I<PolyListNum-1;I++)
378: for(J=0;J<ExpMatColNum;J++)
379: for(K=OneMat[I];K<OneMat[I+1];K++)
380: NormMat2[I][J]-=ExpMat[K][J]$
381:
382: for(I=0;I<PolyListNum-1;I++)
383: NormMat2[I][ExpMatColNum]=OneMat[I+1]-OneMat[I]$
384:
385: ExtVars=Vars$
386: for(I=0;I<PolyListNum-1;I++)
387: ExtVars=append(ExtVars,[uc()])$
388:
389: SolveList=[]$
390: for(I=0;I<ExpMatColNum;I++){
391: TMP=0$
392: for(J=0;J<ExtMatColNum-1;J++)
393: TMP+=NormMat[I][J]*ExtVars[J]$
394:
395: TMP-=NormMat[I][ExtMatColNum-1]$
396: SolveList=cons(TMP,SolveList)$
397: }
398:
399: for(I=0;I<PolyListNum-1;I++){
400: TMP=0$
401: for(J=0;J<ExpMatColNum;J++)
402: TMP+=NormMat2[I][J]*ExtVars[J]$
403:
404: TMP+=NormMat2[I][ExpMatColNum]*ExtVars[I+ExpMatColNum]$
405:
406: SolveList=cons(TMP,SolveList)$
407: }
408:
1.3 ! kimura 409: Rea=vars(SolveList)$
1.2 kimura 410: Res=solve(SolveList,reverse(ExtVars))$
411:
412: if(nonposdegchk(Res)){
1.3 ! kimura 413: Res=getgcd(Res,Rea)$
! 414: TMP1=makeret1(Res,Vars);
! 415: if(car(TMP1)==0){
! 416: TMP2=roundret(cdr(TMP1));
! 417: TMP3=map(drint,cdr(TMP1))$
! 418: RET=cons([cdr(TMP1),newvect(length(TMP3),TMP3),TMP2],RET)$
! 419: }
! 420: else
! 421: RET=cons([cdr(TMP1),[],[]],RET)$
! 422: }
! 423:
! 424: /* second */
! 425:
! 426: NormMat=newmat(ExpMatColNum,ExpMatColNum+1)$
! 427:
! 428: for(I=0;I<ExpMatColNum;I++)
! 429: for(J=0;J<ExpMatColNum;J++)
! 430: for(K=0;K<ExpMatRowNum;K++)
! 431: NormMat[I][J]+=ExpMat[K][I]*ExpMat[K][J]$
1.2 kimura 432:
1.3 ! kimura 433: for(I=0;I<ExpMatColNum;I++)
! 434: for(J=0;J<ExpMatRowNum;J++)
! 435: NormMat[I][ExpMatColNum]+=ExpMat[J][I]$
1.2 kimura 436:
1.3 ! kimura 437: SolveList=[]$
! 438: for(I=0;I<ExpMatColNum;I++){
! 439: TMP=0$
! 440: for(J=0;J<ExpMatColNum;J++)
! 441: TMP+=NormMat[I][J]*Vars[J]$
1.2 kimura 442:
1.3 ! kimura 443: TMP-=NormMat[I][ExpMatColNum]$
! 444: SolveList=cons(TMP,SolveList)$
! 445: }
1.2 kimura 446:
1.3 ! kimura 447: Rea=vars(SolveList)$
! 448: Res=solve(SolveList,Vars)$
1.2 kimura 449:
1.3 ! kimura 450: if(nonposdegchk(Res)){
! 451: Res=getgcd(Res,Rea)$
! 452: TMP1=makeret1(Res,Vars);
! 453: if(car(TMP1)==0){
! 454: TMP2=roundret(cdr(TMP1));
! 455: TMP3=map(drint,cdr(TMP1))$
! 456: RET=cons([cdr(TMP1),newvect(length(TMP3),TMP3),TMP2],RET)$
1.2 kimura 457: }
458: else
1.3 ! kimura 459: RET=cons([cdr(TMP1),[],[]],RET)$
1.2 kimura 460: }
461:
1.3 ! kimura 462: /* third */
! 463:
! 464: ExpMat=qsort(ExpMat,junban2)$
1.2 kimura 465: ExpMat2=[]$
466: for(I=0;I<size(ExpMat)[0];I++)
467: if(car(ExpMat2)!=ExpMat[I])
468: ExpMat2=cons(ExpMat[I],ExpMat2)$
469:
470: ExpMat=newvect(length(ExpMat2),ExpMat2)$
471: ExpMatRowNum=size(ExpMat)[0]$
472: ExpMatColNum=size(ExpMat[0])[0]$
473:
474: NormMat=newmat(ExpMatColNum,ExpMatColNum+1)$
475:
476: for(I=0;I<ExpMatColNum;I++)
477: for(J=0;J<ExpMatColNum;J++)
478: for(K=0;K<ExpMatRowNum;K++)
479: NormMat[I][J]+=ExpMat[K][I]*ExpMat[K][J]$
480:
481: for(I=0;I<ExpMatColNum;I++)
482: for(J=0;J<ExpMatRowNum;J++)
483: NormMat[I][ExpMatColNum]+=ExpMat[J][I]$
484:
485: SolveList=[]$
486: for(I=0;I<ExpMatColNum;I++){
487: TMP=0$
488: for(J=0;J<ExpMatColNum;J++)
489: TMP+=NormMat[I][J]*Vars[J]$
490:
491: TMP-=NormMat[I][ExpMatColNum]$
492: SolveList=cons(TMP,SolveList)$
493: }
494:
1.3 ! kimura 495: Rea=vars(SolveList)$
1.2 kimura 496: Res=solve(SolveList,Vars)$
497:
1.3 ! kimura 498: if(nonposdegchk(Res)){
! 499: Res=getgcd(Res,Rea)$
! 500: TMP1=makeret1(Res,Vars);
! 501: if(car(TMP1)==0){
! 502: TMP2=roundret(cdr(TMP1));
! 503: TMP3=map(drint,cdr(TMP1))$
! 504: RET=cons([cdr(TMP1),newvect(length(TMP3),TMP3),TMP2],RET)$
! 505: }
! 506: else
! 507: RET=cons([cdr(TMP1),[],[]],RET)$
! 508: }
! 509:
! 510: RET=cons(Vars,reverse(RET))$
! 511: RET=cons(0,RET)$
! 512: return RET$
! 513: }
! 514:
! 515: def average(PolyList,Vars){
! 516:
! 517: RET=[]$
! 518: dp_ord(2)$
! 519:
! 520: PolyListNum=length(PolyList)$
! 521:
! 522: ExpMat=[]$
! 523: for(I=0;I<PolyListNum;I++)
! 524: for(Poly=dp_ptod(PolyList[I],Vars);Poly!=0;Poly=dp_rest(Poly))
! 525: ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
! 526:
! 527: ExpMat=reverse(ExpMat)$
! 528: ExpMat=newvect(length(ExpMat),ExpMat)$
! 529:
! 530: ExpMat=qsort(ExpMat,junban2)$
! 531: ExpMat2=[]$
! 532: for(I=0;I<size(ExpMat)[0];I++)
! 533: if(car(ExpMat2)!=ExpMat[I])
! 534: ExpMat2=cons(ExpMat[I],ExpMat2)$
! 535:
! 536: ExpMat=newvect(length(ExpMat2),ExpMat2)$
! 537: ExpMatRowNum=size(ExpMat)[0]$
! 538: ExpMatColNum=size(ExpMat[0])[0]$
! 539:
! 540: Res=newvect(ExpMatColNum);
! 541: for(I=0;I<ExpMatColNum;I++)
! 542: Res[I]=newvect(2,[Vars[I]])$
! 543:
! 544: for(I=0;I<ExpMatRowNum;I++)
! 545: for(J=0;J<ExpMatColNum;J++)
! 546: Res[J][1]+=ExpMat[I][J]$
! 547:
! 548: for(I=0;I<ExpMatColNum;I++)
! 549: if(Res[I][1]==0)
! 550: Res[I][1]=1$
! 551: else
! 552: Res[I][1]=1/Res[I][1]$
! 553:
! 554: RET=cons(makeret(vtol(Res),Vars,1),RET)$
! 555: RET=cons(Vars,RET)$
1.2 kimura 556:
1.3 ! kimura 557: return RET$
1.2 kimura 558: }
559:
560: end$