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