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