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