![[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.18 kimura 4: #define EPS 1E-6
5: #define TINY 1E-20
6: #define MAX_ITER 100
1.19 noro 7: #define ROUND_THRESHOLD 0.4
1.18 kimura 8:
9: def rotate(A,I,J,K,L,C,S){
10:
11: X=A[I][J];
12: Y=A[K][L];
13: A[I][J]=X*C-Y*S;
14: A[K][L]=X*S+Y*C;
15:
16: return 1;
17: }
18:
19: def jacobi(N,A,W){
20:
21: S=OFFDIAG=0.0;
22:
23: for(J=0;J<N;J++){
24:
25: for(K=0;K<N;K++)
26: W[J][K]=0.0;
27:
28: W[J][J]=1.0;
29: S+=A[J][J]*A[J][J];
30:
31: for(K=J+1;K<N;K++)
32: OFFDIAG+=A[J][K]*A[J][K];
33: }
34:
35: TOLERANCE=EPS*EPS*(S/2+OFFDIAG);
36:
37: for(ITER=1;ITER<=MAX_ITER;ITER++){
38:
39: OFFDIAG=0.0;
40: for(J=0;J<N-1;J++)
41: for(K=J+1;K<N;K++)
42: OFFDIAG+=A[J][K]*A[J][K];
43:
44: if(OFFDIAG < TOLERANCE)
45: break;
46:
47: for(J=0;J<N-1;J++){
48: for(K=J+1;K<N;K++){
49:
50: if(dabs(A[J][K])<TINY)
51: continue;
52:
53: T=(A[K][K]-A[J][J])/(2.0*A[J][K]);
54:
55: if(T>=0.0)
56: T=1.0/(T+dsqrt(T*T+1));
57: else
58: T=1.0/(T-dsqrt(T*T+1));
59:
60: C=1.0/dsqrt(T*T+1);
61:
62: S=T*C;
63:
64: T*=A[J][K];
65:
66: A[J][J]-=T;
67: A[K][K]+=T;
68: A[J][K]=0.0;
69:
70: for(I=0;I<J;I++)
71: rotate(A,I,J,I,K,C,S);
72:
73: for(I=J+1;I<K;I++)
74: rotate(A,J,I,I,K,C,S);
75:
76: for(I=K+1;I<N;I++)
77: rotate(A,J,I,K,I,C,S);
78:
79: for(I=0;I<N;I++)
80: rotate(W,J,I,K,I,C,S);
81:
82: }
83: }
84: }
85:
86: if (ITER > MAX_ITER)
87: return 0;
88:
89: for(I=0;I<N-1;I++){
90:
91: K=I;
92:
93: T=A[K][K];
94:
95: for(J=I+1;J<N;J++)
96: if(A[J][J]>T){
97: K=J;
98: T=A[K][K];
99: }
100:
101: A[K][K]=A[I][I];
102:
103: A[I][I]=T;
104:
105: V=W[K];
106:
107: W[K]=W[I];
108:
109: W[I]=V;
110: }
111:
112: return 1;
113: }
114:
1.24 kimura 115: def interval2value(A,Vars){
116:
117: B=atl(A)$
118:
119: if(length(B)>2){
120: print("bug")$
121: return []$
122: }
1.26 kimura 123: else if(length(B)==0){
124: if(fop(A)==0)
1.24 kimura 125: return [Vars,1]$
126: else
127: return []$
128: }
129: else if(length(B)==1){
130:
131: C=fargs(B[0])$
132: D=vars(C)$
133: E=solve(C,D)$
134:
135: if(fop(B[0])==15)
136: return [Vars,E[0][1]+1]$
137: else if(fop(B[0])==11)
138: return [Vars,E[0][1]-1]$
1.25 kimura 139: else if(fop(B[0])==8)
140: return [Vars,E[0][1]]$
1.24 kimura 141: else
142: return []$
143: }
144: else{
145:
146: C=fargs(B[0])$
147: D=vars(C)$
148: E=solve(C,D)$
149:
150: C=fargs(B[1])$
151: D=vars(C)$
152: F=solve(C,D)$
153:
154: return [Vars,(E[0][1]+F[0][1])/2]$
155: }
156:
157: }
158:
159: def fixpointmain(F,Vars){
160:
161: RET=[]$
162: for(I=length(Vars)-1;I>=1;I--){
163:
164: for(H=[],J=0;J<I;J++)
165: H=cons(Vars[J],H)$
166:
167: G=interval2value(qe(ex(H,F)),Vars[I])$
168:
169: if(G==[])
170: return RET$
171: else
172: RET=cons(G,RET)$
173:
174: F=subf(F,G[0],G[1])$
175: }
176:
177: G=interval2value(simpl(F),Vars[0])$
178:
179: if(G==[])
180: return RET$
181: else
182: RET=cons(G,RET)$
183:
184: return RET$
185: }
186:
187:
188: def fixedpoint(A,FLAG){
189:
190: Vars=vars(A)$
191:
192: N=length(A)$
193:
194: if (FLAG==0)
195: for(F=@true,I=0;I < N; I++ ) { F = F @&& A[I] @> 0$ }
196: else if (FLAG==1)
197: for(F=@true,I=0;I < N; I++ ) { F = F @&& A[I] @< 0$ }
198:
199: return fixpointmain(F,Vars)$
200: }
201:
1.8 kimura 202: def junban(A,B){
203: return (A<B ? 1:(A>B ? -1:0))$
204: }
205:
1.10 kimura 206: def bsort(A){
1.8 kimura 207:
1.10 kimura 208: K=size(A)[0]-1$
209: while(K>=0){
210: J=-1$
211: for(I=1;I<=K;I++)
212: if(A[I-1][0]<A[I][0]){
213: J=I-1$
214: X=A[J]$
215: A[J]=A[I]$
216: A[I]=X$
217: }
218: K=J$
219: }
220: return A$
221: }
222:
1.26 kimura 223: def wsort(A,B,C,ID){
1.10 kimura 224:
1.26 kimura 225: D=newvect(length(B))$
226: for(I=0;I<length(B);I++)
227: D[I]=[A[I],B[I],C[I]]$
1.10 kimura 228:
1.26 kimura 229: D=bsort(D)$
1.10 kimura 230:
1.26 kimura 231: for(E=[],I=0;I<length(B);I++)
232: E=cons(D[I][1],E)$
233: E=reverse(E)$
1.10 kimura 234:
1.26 kimura 235: for(F=[],I=0;I<length(B);I++)
236: F=cons(D[I][2],F)$
237: F=reverse(F)$
1.8 kimura 238:
1.26 kimura 239: return [[ID,E,F]]$
1.8 kimura 240: }
241:
1.4 kimura 242: def nonposdegchk(Res){
243:
244: for(I=0;I<length(Res);I++)
245: if(Res[I][1]<=0)
246: return 0$
247:
248: return 1$
249: }
250:
1.8 kimura 251: def getgcd(A,B){
252:
1.31 ! kimura 253: Anum=length(A)$
1.8 kimura 254:
1.31 ! kimura 255: TMP=[]$
! 256: for(I=0;I<length(B);I++){
1.8 kimura 257:
1.31 ! kimura 258: for(J=0;J<Anum;J++)
! 259: if(B[I]==A[J][0])
! 260: break;
1.8 kimura 261:
1.31 ! kimura 262: if(J==Anum)
! 263: TMP=cons([B[I],B[I]],TMP)$
! 264: }
! 265: A=append(A,TMP)$
1.8 kimura 266:
1.31 ! kimura 267: Anum=length(A)$
! 268: A=map(ltov,A)$
1.8 kimura 269:
1.31 ! kimura 270: for(D=0,I=0;I<Anum;I++)
! 271: D=gcd(D,A[I][1])$
1.8 kimura 272:
273: if(D!=0){
1.31 ! kimura 274: for(I=0;I<Anum;I++)
! 275: A[I][1]=red(A[I][1]/D)$
1.8 kimura 276: }
1.6 kimura 277:
1.31 ! kimura 278: for(L=1,D=0,I=0;I<Anum;I++){
! 279: if(type(TMP=dn(A[I][1]))==1)
1.8 kimura 280: L=ilcm(L,TMP)$
281:
1.31 ! kimura 282: if(type(TMP=nm(A[I][1]))==1)
1.8 kimura 283: D=igcd(D,TMP)$
284: }
285:
1.31 ! kimura 286: for(I=0;I<Anum;I++)
! 287: A[I][1]=A[I][1]*L$
! 288:
! 289: if(D!=0){
1.6 kimura 290:
1.31 ! kimura 291: for(I=0;I<Anum;I++)
! 292: A[I][1]=A[I][1]/D$
! 293: }
1.6 kimura 294:
1.31 ! kimura 295: return map(vtol,A)$
1.6 kimura 296: }
297:
1.10 kimura 298: def makeret(Res,Vars,FLAG){
1.4 kimura 299:
1.8 kimura 300: ResNum=length(Res)$
1.4 kimura 301: VarsNum=length(Vars)$
302:
1.17 kimura 303: ResVec=newvect(ResNum)$
1.4 kimura 304:
1.29 kimura 305: if(FLAG)
306: M=0$
307: else
308: M=-1$
309:
310: for(I=0;I<ResNum;I++){
1.17 kimura 311: if(member(Res[I][0],Vars)){
312: ResVec[I]=Res[I][1]$
1.16 kimura 313:
1.29 kimura 314: if(FLAG){
315: if(type(ResVec[I])==1){
316: if(M==0)
317: M=ResVec[I]$
318: else
319: if(ResVec[I]<M)
320: M=ResVec[I]$
321: }
322: else
323: M=-1$
324: }
325: }
1.17 kimura 326: }
1.16 kimura 327:
1.29 kimura 328:
329: if(M!=-1)
1.8 kimura 330: ResVec=ResVec/M;
1.4 kimura 331:
1.17 kimura 332: RET=newvect(VarsNum,Vars)$
1.4 kimura 333:
1.8 kimura 334: for(I=0;I<ResNum;I++){
335: for(J=0;J<VarsNum;J++)
336: if(Vars[J]==Res[I][0])
337: break$
1.4 kimura 338:
1.8 kimura 339: if(J<VarsNum)
340: RET[J]=ResVec[I]$
341: }
342:
1.4 kimura 343: for(I=0;I<VarsNum;I++)
1.8 kimura 344: if(type(RET[I])!=1)
345: return [1,RET]$
1.4 kimura 346:
1.8 kimura 347: return [0,RET]$
1.4 kimura 348: }
349:
350: def roundret(V){
351:
1.8 kimura 352: VN=size(V)[0]$
1.4 kimura 353:
1.8 kimura 354: RET0=V$
1.13 kimura 355: for(I=1;I<1000;I++){
1.4 kimura 356: RET1=I*RET0$
357: for(J=0;J<VN;J++){
358: X=drint(RET1[J])$
1.19 noro 359: if(dabs(X-RET1[J])<ROUND_THRESHOLD)
1.4 kimura 360: RET1[J]=X$
361: else
362: break$
363: }
364: if(J==VN)
365: break$
366: }
367:
368: if(I==1000)
369: return []$
370: else
371: return RET1$
372: }
373:
374: def chkou(L,ExpMat,CHAGORD){
375:
1.8 kimura 376: for(P=1,I=0;I<L;I++){
1.4 kimura 377: Q=ExpMat[L][CHAGORD[I]]$
378: for(J=0;J<size(ExpMat[0])[0];J++){
379: ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
380: *ExpMat[L][CHAGORD[J]]-
381: Q*ExpMat[I][CHAGORD[J]])/P)$
382: }
383:
384: P=ExpMat[I][CHAGORD[I]]$
385: }
386:
387: for(J=0;J<size(ExpMat[0])[0];J++)
388: if(ExpMat[L][CHAGORD[J]]!=0)
389: break$
390:
391: if(J==size(ExpMat[0])[0])
392: return L$
393: else{
394: TMP=CHAGORD[L]$
395: CHAGORD[L]=CHAGORD[J]$
396: CHAGORD[J]=TMP$
397: return (L+1)$
398: }
399: }
400:
1.8 kimura 401: def qcheckmain(PolyList,Vars){
1.4 kimura 402:
403: RET=[]$
404: PolyListNum=length(PolyList)$
405: VarsNum=length(Vars)$
406:
407: ExpMat=newvect(VarsNum)$
408: CHAGORD=newvect(VarsNum)$
409: for(I=0;I<VarsNum;I++)
410: CHAGORD[I]=I$
411:
412: L=0$
413: for(I=0;I<PolyListNum;I++){
414: Poly=dp_ptod(PolyList[I],Vars)$
415: BASE0=dp_etov(dp_ht(Poly))$
416: Poly=dp_rest(Poly)$
417: for(;Poly!=0;Poly=dp_rest(Poly)){
418: ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
419: L=chkou(L,ExpMat,CHAGORD)$
1.8 kimura 420: if(L==VarsNum-1)
421: return [L,CHAGORD,ExpMat]$
1.4 kimura 422: }
423: }
424:
1.8 kimura 425: return [L,CHAGORD,ExpMat]$
1.4 kimura 426: }
427:
428: def inner(A,B){
429:
430: SUM=0$
431: for(I=0;I<size(A)[0];I++)
432: SUM+=A[I]*B[I]$
433:
434: return SUM$
435: }
436:
437: def checktd(PolyList,Vars,ResVars){
438:
439: PolyListNum=length(PolyList)$
440: VarsNum=length(Vars)$
441:
442: L=0$
443: for(I=0;I<PolyListNum;I++){
444: Poly=dp_ptod(PolyList[I],Vars)$
445: J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
446: Poly=dp_rest(Poly)$
447: for(;Poly!=0;Poly=dp_rest(Poly))
448: if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
449: return 0$
450: }
451:
452: return 1$
453: }
454:
1.24 kimura 455: def value2(Vars,Ans){
456:
457: N=length(Vars)$
458: Res=newvect(N)$
459: for(I=0;I<N;I++){
460: Res[I]=newvect(2)$
461: Res[I][0]=Vars[I]$
462: Res[I][1]=Ans[I]$
463: }
1.31 ! kimura 464: Res=map(vtol,Res)$
! 465: Res=vtol(Res)$
1.24 kimura 466:
467: Res=getgcd(Res,Vars)$
468:
469: if(nonposdegchk(Res)){
470: TMP1=makeret(Res,Vars,1)$
471: return vtol(TMP1[1])$
472: }
473: else
474: return []$
475: }
476:
1.10 kimura 477: def qcheck(PolyList,Vars,FLAG){
1.4 kimura 478:
1.10 kimura 479: RET=[]$
1.8 kimura 480: Res=qcheckmain(PolyList,Vars)$
1.4 kimura 481: VarsNum=length(Vars)$
482:
483: IndNum=Res[0]$
484: CHAGORD=Res[1]$
485: ExpMat=Res[2]$
486:
487: SolveList=[]$
488: for(I=0;I<IndNum;I++){
489: TMP=0$
490: for(J=0;J<VarsNum;J++)
491: TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$
492:
493: SolveList=cons(TMP,SolveList)$
494: }
495:
1.8 kimura 496: Rea=vars(SolveList)$
497:
1.4 kimura 498: VarsList=[]$
499: for(I=0;I<VarsNum;I++)
1.31 ! kimura 500: if(member(TMP0=Vars[CHAGORD[I]],Rea))
1.8 kimura 501: VarsList=cons(Vars[CHAGORD[I]],VarsList)$
1.4 kimura 502:
503: Res=solve(reverse(SolveList),reverse(VarsList))$
1.8 kimura 504: Res=getgcd(Res,Rea)$
1.4 kimura 505:
506: if(nonposdegchk(Res)){
507:
1.26 kimura 508: TMP1=makeret(Res,Vars,0)$
1.4 kimura 509:
1.26 kimura 510: if(checktd(PolyList,Vars,TMP1[1])==1){
1.18 kimura 511:
1.26 kimura 512: if(FLAG==0){
1.24 kimura 513:
1.26 kimura 514: if(TMP1[0]==0)
515: RET=append(RET,wsort(TMP1[1],Vars,TMP1[1],0))$
516: else{
1.24 kimura 517:
1.26 kimura 518: TMP=vtol(TMP1[1])$
1.27 kimura 519: RET0=[]$
1.26 kimura 520: if((TMP0=fixedpoint(TMP,0))!=[]){
1.24 kimura 521:
1.26 kimura 522: for(I=0;I<length(TMP0);I++)
523: TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.27 kimura 524: RET0=value2(Vars,TMP)$
525: if(RET0!=[])
526: RET0=wsort(RET0,Vars,RET0,1/10)$
527: }
528:
529: TMP=vtol(TMP1[1])$
530: if(RET0==[] && (TMP0=fixedpoint(TMP,1))!=[]){
1.24 kimura 531:
1.26 kimura 532: for(I=0;I<length(TMP0);I++)
533: TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.27 kimura 534: RET0=value2(Vars,TMP)$
1.26 kimura 535:
1.27 kimura 536: if(RET0!=[])
537: RET0=wsort(RET0,Vars,RET0,1/10)$
1.26 kimura 538: }
1.27 kimura 539: RET=append(RET,RET0)$
1.26 kimura 540: }
1.10 kimura 541: }
1.26 kimura 542: else if(FLAG==1)
543: RET=append(RET,[[0,Vars,vtol(TMP1[1])]])$
1.4 kimura 544: }
545: }
546:
1.26 kimura 547: return RET$
1.4 kimura 548: }
549:
1.18 kimura 550: def leastsq(NormMat,ExpMat,Vars,FLAG,ID){
1.8 kimura 551:
552: RET=[]$
1.4 kimura 553:
554: ExpMatRowNum=size(ExpMat)[0]$
555: ExpMatColNum=size(ExpMat[0])[0]$
556:
1.8 kimura 557: if(NormMat==0){
558: NormMat=newmat(ExpMatColNum,ExpMatColNum)$
559:
560: for(I=0;I<ExpMatColNum;I++)
561: for(J=I;J<ExpMatColNum;J++)
562: for(K=0;K<ExpMatRowNum;K++)
563: NormMat[I][J]+=
564: ExpMat[K][I]*ExpMat[K][J]$
565: }
1.4 kimura 566:
1.8 kimura 567: BVec=newvect(ExpMatColNum)$
1.4 kimura 568:
569: for(I=0;I<ExpMatColNum;I++)
570: for(J=0;J<ExpMatRowNum;J++)
1.8 kimura 571: BVec[I]+=ExpMat[J][I]$
1.4 kimura 572:
573: SolveList=[]$
574: for(I=0;I<ExpMatColNum;I++){
575: TMP=0$
1.8 kimura 576: for(J=0;J<I;J++)
577: TMP+=NormMat[J][I]*Vars[J]$
578:
579: for(J=I;J<ExpMatColNum;J++)
1.4 kimura 580: TMP+=NormMat[I][J]*Vars[J]$
581:
1.8 kimura 582: TMP-=BVec[I]$
1.4 kimura 583: SolveList=cons(TMP,SolveList)$
584: }
585:
586: Rea=vars(SolveList)$
1.8 kimura 587:
588: VarsList=[]$
589: for(I=0;I<length(Vars);I++)
590: if(member(Vars[I],Rea))
591: VarsList=cons(Vars[I],VarsList)$
592:
593: Res=solve(SolveList,VarsList)$
594: Res=getgcd(Res,Rea)$
1.4 kimura 595:
596: if(nonposdegchk(Res)){
1.16 kimura 597:
1.10 kimura 598: TMP1=makeret(Res,Vars,1)$
1.16 kimura 599:
1.26 kimura 600: if(FLAG==0){
1.24 kimura 601:
1.26 kimura 602: if(TMP1[0]==0){
1.8 kimura 603:
1.26 kimura 604: TMP=roundret(TMP1[1])$
1.18 kimura 605:
1.26 kimura 606: RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]*1.0),ID))$
1.8 kimura 607:
1.26 kimura 608: if(TMP!=[])
609: RET=append(RET,wsort(TMP1[1],Vars,TMP,ID+1))$
610: }
611: else{
1.24 kimura 612:
1.26 kimura 613: TMP=vtol(TMP1[1])$
1.27 kimura 614: RET0=[]$
615: if((TMP0=fixedpoint(TMP,0))!=[]){
1.24 kimura 616:
1.26 kimura 617: for(I=0;I<length(TMP0);I++)
618: TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.27 kimura 619: RET0=value2(Vars,TMP)$
620: if(RET0!=[])
1.28 kimura 621: RET0=wsort(RET0,Vars,RET0,ID+1/10)$
1.27 kimura 622: }
623:
624: TMP=vtol(TMP1[1])$
625: if(RET0==[] && (TMP0=fixedpoint(TMP,1))!=[]){
626:
1.26 kimura 627: for(I=0;I<length(TMP0);I++)
628: TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.27 kimura 629: RET0=value2(Vars,TMP)$
630:
631: if(RET0!=[])
1.28 kimura 632: RET0=wsort(RET0,Vars,RET0,ID+1/10)$
1.27 kimura 633: }
1.24 kimura 634:
1.27 kimura 635: RET=append(RET,RET0)$
1.26 kimura 636: }
1.24 kimura 637:
1.8 kimura 638: }
1.26 kimura 639: else if(FLAG==1)
640: RET=append(RET,[[ID,Vars,vtol(TMP1[1])]])$
1.8 kimura 641: }
642:
1.26 kimura 643: return RET$
1.8 kimura 644: }
645:
1.18 kimura 646: def unitweight(ExpMat,Vars,PolyListNum,OneMat,FLAG){
1.8 kimura 647:
648: RET=[]$
649:
650: ExpMatRowNum=size(ExpMat)[0]$
651: ExpMatColNum=size(ExpMat[0])[0]$
652: ExtMatColNum=ExpMatColNum+PolyListNum$
653:
654: ExtVars=reverse(Vars)$
655: for(I=0;I<PolyListNum;I++)
656: ExtVars=cons(uc(),ExtVars)$
657:
658: ExtVars=reverse(ExtVars)$
659:
1.18 kimura 660: NormMat0=newvect(ExpMatColNum)$
661: for(I=0;I<ExpMatColNum;I++)
662: NormMat0[I]=newvect(ExpMatColNum)$
1.8 kimura 663:
664: for(I=0;I<ExpMatColNum;I++)
665: for(J=I;J<ExpMatColNum;J++)
666: for(K=0;K<ExpMatRowNum;K++)
1.18 kimura 667: NormMat0[I][J]+=
1.8 kimura 668: ExpMat[K][I]*
669: ExpMat[K][J]$
670:
1.18 kimura 671: NormMat1=newvect(ExtMatColNum)$
672: for(I=0;I<ExtMatColNum;I++)
673: NormMat1[I]=newvect(ExtMatColNum)$
674:
675:
676: WorkMat=newvect(ExtMatColNum)$
677: for(I=0;I<ExtMatColNum;I++)
678: WorkMat[I]=newvect(ExtMatColNum)$
679:
680:
681: for(I=0;I<ExpMatColNum;I++)
682: for(J=I;J<ExpMatColNum;J++)
683: NormMat1[I][J]=NormMat0[I][J]$
684:
1.8 kimura 685: for(I=0;I<ExpMatColNum;I++)
686: for(J=0;J<PolyListNum;J++)
687: for(K=OneMat[J];K<OneMat[J+1];K++)
1.26 kimura 688: NormMat1[I][J+ExpMatColNum]-=ExpMat[K][I]$
1.8 kimura 689:
690: for(I=0;I<PolyListNum;I++)
1.18 kimura 691: NormMat1[I+ExpMatColNum][I+ExpMatColNum]=OneMat[I+1]-OneMat[I]$
692:
693: if(jacobi(ExtMatColNum,NormMat1,WorkMat)){
1.8 kimura 694:
1.30 kimura 695: Res=newvect(ExtMatColNum)$
696: for(I=0;I<ExtMatColNum;I++){
1.18 kimura 697: Res[I]=newvect(2)$
1.30 kimura 698: Res[I][0]=ExtVars[I]$
1.18 kimura 699: Res[I][1]=WorkMat[ExtMatColNum-1][I]$
700: }
1.8 kimura 701:
702: if(nonposdegchk(Res)){
703:
1.10 kimura 704: TMP1=makeret(Res,Vars,1)$
1.16 kimura 705:
1.26 kimura 706: if(FLAG==0){
707: TMP=roundret(TMP1[1])$
1.8 kimura 708:
1.26 kimura 709: RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]*1.0),1))$
1.18 kimura 710:
1.26 kimura 711: if(TMP!=[])
712: RET=append(RET,wsort(TMP1[1],Vars,TMP,2))$
713: }
1.27 kimura 714: else if(FLAG==1)
1.26 kimura 715: RET=append(RET,[[1,Vars,vtol(TMP1[1])]])$
1.4 kimura 716: }
1.18 kimura 717: }
718:
719: return [NormMat0,RET]$
1.6 kimura 720: }
721:
1.18 kimura 722: def weight(PolyList,Vars,FLAG){
1.6 kimura 723:
724: Vars0=vars(PolyList)$
725: Vars1=[]$
726: for(I=0;I<length(Vars);I++)
727: if(member(Vars[I],Vars0))
728: Vars1=cons(Vars[I],Vars1)$
729:
730: Vars=reverse(Vars1)$
731:
732: RET=[]$
733:
1.18 kimura 734: TMP=qcheck(PolyList,Vars,FLAG)$
1.6 kimura 735:
1.8 kimura 736: if(TMP!=[]){
737: RET=append(RET,TMP)$
1.18 kimura 738: return RET$
1.6 kimura 739: }
1.4 kimura 740:
1.6 kimura 741: dp_ord(2)$
1.4 kimura 742:
1.6 kimura 743: PolyListNum=length(PolyList)$
1.4 kimura 744:
1.18 kimura 745: OneMat=newvect(PolyListNum+1,[0])$
746: ExpMat=[]$
747: for(I=0;I<PolyListNum;I++){
748: for(Poly=dp_ptod(PolyList[I],Vars);
749: Poly!=0;Poly=dp_rest(Poly)){
750: ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
1.8 kimura 751: }
1.18 kimura 752: OneMat[I+1]=length(ExpMat)$
753: }
1.4 kimura 754:
1.18 kimura 755: ExpMat=reverse(ExpMat)$
756: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 757:
1.18 kimura 758: TMP=unitweight(ExpMat,Vars,PolyListNum,OneMat,FLAG)$
1.8 kimura 759:
1.18 kimura 760: RET=append(RET,TMP[1])$
1.4 kimura 761:
1.20 kimura 762: TMP0=leastsq(TMP[0],ExpMat,Vars,FLAG,3)$
763:
764: RET=append(RET,TMP0)$
1.4 kimura 765:
1.6 kimura 766: ExpMat=qsort(ExpMat,junban)$
1.8 kimura 767:
1.6 kimura 768: ExpMat2=[]$
769: for(I=0;I<size(ExpMat)[0];I++)
770: if(car(ExpMat2)!=ExpMat[I])
771: ExpMat2=cons(ExpMat[I],ExpMat2)$
1.4 kimura 772:
1.6 kimura 773: if(size(ExpMat)[0]!=length(ExpMat2)){
774: ExpMat=newvect(length(ExpMat2),ExpMat2)$
1.18 kimura 775: RET=append(RET,leastsq(0,ExpMat,Vars,FLAG,5))$
1.20 kimura 776: }
777: else{
778: TMP0=map(ltov,TMP0)$
779:
780: for(I=0;I<length(TMP0);I++)
1.24 kimura 781: if(TMP0[I][0]==3)
782: TMP0[I][0]=5$
783: else if(TMP0[I][0]==4)
784: TMP0[I][0]=6$
1.20 kimura 785:
786: TMP0=map(vtol,TMP0)$
787:
788: RET=append(RET,TMP0)$
1.4 kimura 789: }
790:
1.18 kimura 791: return RET$
1.4 kimura 792: }
793:
794: end$