![[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.9 kimura 202: def nonzerovec(A){
203:
204: for(I=0;I<size(A)[0];I++)
205: if(A[I]!=0)
206: return 1$
207:
208: return 0$
209: }
210:
1.8 kimura 211: def junban(A,B){
212: return (A<B ? 1:(A>B ? -1:0))$
213: }
214:
215: def worder(A,B){
216: return (A[0]<B[0] ? 1:(A[0]>B[0] ? -1:0))$
217: }
218:
1.10 kimura 219: def bsort(A){
1.8 kimura 220:
1.10 kimura 221: K=size(A)[0]-1$
222: while(K>=0){
223: J=-1$
224: for(I=1;I<=K;I++)
225: if(A[I-1][0]<A[I][0]){
226: J=I-1$
227: X=A[J]$
228: A[J]=A[I]$
229: A[I]=X$
230: }
231: K=J$
232: }
233: return A$
234: }
235:
1.26 kimura 236: def wsort(A,B,C,ID){
1.10 kimura 237:
1.26 kimura 238: D=newvect(length(B))$
239: for(I=0;I<length(B);I++)
240: D[I]=[A[I],B[I],C[I]]$
1.10 kimura 241:
1.26 kimura 242: D=bsort(D)$
1.10 kimura 243:
1.26 kimura 244: for(E=[],I=0;I<length(B);I++)
245: E=cons(D[I][1],E)$
246: E=reverse(E)$
1.10 kimura 247:
1.26 kimura 248: for(F=[],I=0;I<length(B);I++)
249: F=cons(D[I][2],F)$
250: F=reverse(F)$
1.8 kimura 251:
1.26 kimura 252: return [[ID,E,F]]$
1.8 kimura 253: }
254:
1.4 kimura 255: def nonposdegchk(Res){
256:
257: for(I=0;I<length(Res);I++)
258: if(Res[I][1]<=0)
259: return 0$
260:
261: return 1$
262: }
263:
1.8 kimura 264: def getgcd(A,B){
265:
266: VarsNumA=length(A)$
267: VarsNumB=length(B)$
268:
269: C=newvect(VarsNumB,B)$
270:
271: for(I=0;I<VarsNumA;I++){
272:
273: for(J=0;J<VarsNumB;J++)
1.16 kimura 274: if(B[J]==A[I][0])
1.8 kimura 275: break$
276:
277: if(J<VarsNumB)
278: C[J]=A[I][1]$
279: }
280:
281: D=0$
282: for(I=0;I<VarsNumB;I++)
283: D=gcd(D,C[I])$
284:
285: if(D!=0){
286: C=C/D$
287: C=map(red,C)$
288: }
1.6 kimura 289:
1.8 kimura 290: for(L=1,D=0,I=0;I<VarsNumB;I++){
291: if(type(TMP=dn(C[I]))==1)
292: L=ilcm(L,TMP)$
293:
294: if(type(TMP=nm(C[I]))==1)
295: D=igcd(D,TMP)$
296: }
297:
298: C=C*L$
299: if(D!=0)
300: C=C/D$
1.6 kimura 301:
1.8 kimura 302: RET=[]$
303: for(I=0;I<VarsNumB;I++)
304: RET=cons([B[I],C[I]],RET)$
1.6 kimura 305:
1.8 kimura 306: return RET$
1.6 kimura 307: }
308:
1.10 kimura 309: def makeret(Res,Vars,FLAG){
1.4 kimura 310:
1.8 kimura 311: ResNum=length(Res)$
1.4 kimura 312: VarsNum=length(Vars)$
313:
1.17 kimura 314: ResVec=newvect(ResNum)$
1.4 kimura 315:
1.29 kimura 316: if(FLAG)
317: M=0$
318: else
319: M=-1$
320:
321: for(I=0;I<ResNum;I++){
1.17 kimura 322: if(member(Res[I][0],Vars)){
323: ResVec[I]=Res[I][1]$
1.16 kimura 324:
1.29 kimura 325: if(FLAG){
326: if(type(ResVec[I])==1){
327: if(M==0)
328: M=ResVec[I]$
329: else
330: if(ResVec[I]<M)
331: M=ResVec[I]$
332: }
333: else
334: M=-1$
335: }
336: }
1.17 kimura 337: }
1.16 kimura 338:
1.29 kimura 339:
340: if(M!=-1)
1.8 kimura 341: ResVec=ResVec/M;
1.4 kimura 342:
1.17 kimura 343: RET=newvect(VarsNum,Vars)$
1.4 kimura 344:
1.8 kimura 345: for(I=0;I<ResNum;I++){
346: for(J=0;J<VarsNum;J++)
347: if(Vars[J]==Res[I][0])
348: break$
1.4 kimura 349:
1.8 kimura 350: if(J<VarsNum)
351: RET[J]=ResVec[I]$
352: }
353:
1.10 kimura 354:
355: for(J=0;J<length(Vars);J++)
356: RET=map(subst,RET,Vars[J],
357: strtov(rtostr(Vars[J])+"_deg"))$
358:
1.4 kimura 359: for(I=0;I<VarsNum;I++)
1.8 kimura 360: if(type(RET[I])!=1)
361: return [1,RET]$
1.4 kimura 362:
1.8 kimura 363: return [0,RET]$
1.4 kimura 364: }
365:
366: def roundret(V){
367:
1.8 kimura 368: VN=size(V)[0]$
1.4 kimura 369:
1.8 kimura 370: RET0=V$
1.13 kimura 371: for(I=1;I<1000;I++){
1.4 kimura 372: RET1=I*RET0$
373: for(J=0;J<VN;J++){
374: X=drint(RET1[J])$
1.19 noro 375: if(dabs(X-RET1[J])<ROUND_THRESHOLD)
1.4 kimura 376: RET1[J]=X$
377: else
378: break$
379: }
380: if(J==VN)
381: break$
382: }
383:
384: if(I==1000)
385: return []$
386: else
387: return RET1$
388: }
389:
390: def chkou(L,ExpMat,CHAGORD){
391:
1.8 kimura 392: for(P=1,I=0;I<L;I++){
1.4 kimura 393: Q=ExpMat[L][CHAGORD[I]]$
394: for(J=0;J<size(ExpMat[0])[0];J++){
395: ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
396: *ExpMat[L][CHAGORD[J]]-
397: Q*ExpMat[I][CHAGORD[J]])/P)$
398: }
399:
400: P=ExpMat[I][CHAGORD[I]]$
401: }
402:
403: for(J=0;J<size(ExpMat[0])[0];J++)
404: if(ExpMat[L][CHAGORD[J]]!=0)
405: break$
406:
407: if(J==size(ExpMat[0])[0])
408: return L$
409: else{
410: TMP=CHAGORD[L]$
411: CHAGORD[L]=CHAGORD[J]$
412: CHAGORD[J]=TMP$
413: return (L+1)$
414: }
415: }
416:
1.8 kimura 417: def qcheckmain(PolyList,Vars){
1.4 kimura 418:
419: RET=[]$
420: PolyListNum=length(PolyList)$
421: VarsNum=length(Vars)$
422:
423: ExpMat=newvect(VarsNum)$
424: CHAGORD=newvect(VarsNum)$
425: for(I=0;I<VarsNum;I++)
426: CHAGORD[I]=I$
427:
428: L=0$
429: for(I=0;I<PolyListNum;I++){
430: Poly=dp_ptod(PolyList[I],Vars)$
431: BASE0=dp_etov(dp_ht(Poly))$
432: Poly=dp_rest(Poly)$
433: for(;Poly!=0;Poly=dp_rest(Poly)){
434: ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
435: L=chkou(L,ExpMat,CHAGORD)$
1.8 kimura 436: if(L==VarsNum-1)
437: return [L,CHAGORD,ExpMat]$
1.4 kimura 438: }
439: }
440:
1.8 kimura 441: return [L,CHAGORD,ExpMat]$
1.4 kimura 442: }
443:
444: def inner(A,B){
445:
446: SUM=0$
447: for(I=0;I<size(A)[0];I++)
448: SUM+=A[I]*B[I]$
449:
450: return SUM$
451: }
452:
453: def checktd(PolyList,Vars,ResVars){
454:
455: PolyListNum=length(PolyList)$
456: VarsNum=length(Vars)$
457:
458: L=0$
459: for(I=0;I<PolyListNum;I++){
460: Poly=dp_ptod(PolyList[I],Vars)$
461: J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
462: Poly=dp_rest(Poly)$
463: for(;Poly!=0;Poly=dp_rest(Poly))
464: if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
465: return 0$
466: }
467:
468: return 1$
469: }
470:
1.24 kimura 471: def value2(Vars,Ans){
472:
473: N=length(Vars)$
474: Res=newvect(N)$
475: for(I=0;I<N;I++){
476: Res[I]=newvect(2)$
477: Res[I][0]=Vars[I]$
478: Res[I][1]=Ans[I]$
479: }
480:
481: Res=getgcd(Res,Vars)$
482:
483: if(nonposdegchk(Res)){
484: TMP1=makeret(Res,Vars,1)$
485: return vtol(TMP1[1])$
486: }
487: else
488: return []$
489: }
490:
1.10 kimura 491: def qcheck(PolyList,Vars,FLAG){
1.4 kimura 492:
1.10 kimura 493: RET=[]$
1.8 kimura 494: Res=qcheckmain(PolyList,Vars)$
1.4 kimura 495: VarsNum=length(Vars)$
496:
497: IndNum=Res[0]$
498: CHAGORD=Res[1]$
499: ExpMat=Res[2]$
500:
501: SolveList=[]$
502: for(I=0;I<IndNum;I++){
503: TMP=0$
504: for(J=0;J<VarsNum;J++)
505: TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$
506:
507: SolveList=cons(TMP,SolveList)$
508: }
509:
1.8 kimura 510: Rea=vars(SolveList)$
511:
1.4 kimura 512: VarsList=[]$
513: for(I=0;I<VarsNum;I++)
1.8 kimura 514: if(member(Vars[CHAGORD[I]],Rea))
515: VarsList=cons(Vars[CHAGORD[I]],VarsList)$
1.4 kimura 516:
517: Res=solve(reverse(SolveList),reverse(VarsList))$
1.8 kimura 518: Res=getgcd(Res,Rea)$
1.4 kimura 519:
520: if(nonposdegchk(Res)){
521:
1.26 kimura 522: TMP1=makeret(Res,Vars,0)$
1.4 kimura 523:
1.26 kimura 524: if(checktd(PolyList,Vars,TMP1[1])==1){
1.18 kimura 525:
1.26 kimura 526: if(FLAG==0){
1.24 kimura 527:
1.26 kimura 528: if(TMP1[0]==0)
529: RET=append(RET,wsort(TMP1[1],Vars,TMP1[1],0))$
530: else{
1.24 kimura 531:
1.26 kimura 532: TMP=vtol(TMP1[1])$
1.27 kimura 533: RET0=[]$
1.26 kimura 534: if((TMP0=fixedpoint(TMP,0))!=[]){
1.24 kimura 535:
1.26 kimura 536: for(I=0;I<length(TMP0);I++)
537: TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.27 kimura 538: RET0=value2(Vars,TMP)$
539: if(RET0!=[])
540: RET0=wsort(RET0,Vars,RET0,1/10)$
541: }
542:
543: TMP=vtol(TMP1[1])$
544: if(RET0==[] && (TMP0=fixedpoint(TMP,1))!=[]){
1.24 kimura 545:
1.26 kimura 546: for(I=0;I<length(TMP0);I++)
547: TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.27 kimura 548: RET0=value2(Vars,TMP)$
1.26 kimura 549:
1.27 kimura 550: if(RET0!=[])
551: RET0=wsort(RET0,Vars,RET0,1/10)$
1.26 kimura 552: }
1.27 kimura 553: RET=append(RET,RET0)$
1.26 kimura 554: }
1.10 kimura 555: }
1.26 kimura 556: else if(FLAG==1)
557: RET=append(RET,[[0,Vars,vtol(TMP1[1])]])$
1.4 kimura 558: }
559: }
560:
1.26 kimura 561: return RET$
1.4 kimura 562: }
563:
1.18 kimura 564: def leastsq(NormMat,ExpMat,Vars,FLAG,ID){
1.8 kimura 565:
566: RET=[]$
1.4 kimura 567:
568: ExpMatRowNum=size(ExpMat)[0]$
569: ExpMatColNum=size(ExpMat[0])[0]$
570:
1.8 kimura 571: if(NormMat==0){
572: NormMat=newmat(ExpMatColNum,ExpMatColNum)$
573:
574: for(I=0;I<ExpMatColNum;I++)
575: for(J=I;J<ExpMatColNum;J++)
576: for(K=0;K<ExpMatRowNum;K++)
577: NormMat[I][J]+=
578: ExpMat[K][I]*ExpMat[K][J]$
579: }
1.4 kimura 580:
1.8 kimura 581: BVec=newvect(ExpMatColNum)$
1.4 kimura 582:
583: for(I=0;I<ExpMatColNum;I++)
584: for(J=0;J<ExpMatRowNum;J++)
1.8 kimura 585: BVec[I]+=ExpMat[J][I]$
1.4 kimura 586:
587: SolveList=[]$
588: for(I=0;I<ExpMatColNum;I++){
589: TMP=0$
1.8 kimura 590: for(J=0;J<I;J++)
591: TMP+=NormMat[J][I]*Vars[J]$
592:
593: for(J=I;J<ExpMatColNum;J++)
1.4 kimura 594: TMP+=NormMat[I][J]*Vars[J]$
595:
1.8 kimura 596: TMP-=BVec[I]$
1.4 kimura 597: SolveList=cons(TMP,SolveList)$
598: }
599:
600: Rea=vars(SolveList)$
1.8 kimura 601:
602: VarsList=[]$
603: for(I=0;I<length(Vars);I++)
604: if(member(Vars[I],Rea))
605: VarsList=cons(Vars[I],VarsList)$
606:
607: Res=solve(SolveList,VarsList)$
608: Res=getgcd(Res,Rea)$
1.4 kimura 609:
610: if(nonposdegchk(Res)){
1.16 kimura 611:
1.10 kimura 612: TMP1=makeret(Res,Vars,1)$
1.16 kimura 613:
1.26 kimura 614: if(FLAG==0){
1.24 kimura 615:
1.26 kimura 616: if(TMP1[0]==0){
1.8 kimura 617:
1.26 kimura 618: TMP=roundret(TMP1[1])$
1.18 kimura 619:
1.26 kimura 620: RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]*1.0),ID))$
1.8 kimura 621:
1.26 kimura 622: if(TMP!=[])
623: RET=append(RET,wsort(TMP1[1],Vars,TMP,ID+1))$
624: }
625: else{
1.24 kimura 626:
1.26 kimura 627: TMP=vtol(TMP1[1])$
1.27 kimura 628: RET0=[]$
629: if((TMP0=fixedpoint(TMP,0))!=[]){
1.24 kimura 630:
1.26 kimura 631: for(I=0;I<length(TMP0);I++)
632: TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.27 kimura 633: RET0=value2(Vars,TMP)$
634: if(RET0!=[])
1.28 kimura 635: RET0=wsort(RET0,Vars,RET0,ID+1/10)$
1.27 kimura 636: }
637:
638: TMP=vtol(TMP1[1])$
639: if(RET0==[] && (TMP0=fixedpoint(TMP,1))!=[]){
640:
1.26 kimura 641: for(I=0;I<length(TMP0);I++)
642: TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
1.27 kimura 643: RET0=value2(Vars,TMP)$
644:
645: if(RET0!=[])
1.28 kimura 646: RET0=wsort(RET0,Vars,RET0,ID+1/10)$
1.27 kimura 647: }
1.24 kimura 648:
1.27 kimura 649: RET=append(RET,RET0)$
1.26 kimura 650: }
1.24 kimura 651:
1.8 kimura 652: }
1.26 kimura 653: else if(FLAG==1)
654: RET=append(RET,[[ID,Vars,vtol(TMP1[1])]])$
1.8 kimura 655: }
656:
1.26 kimura 657: return RET$
1.8 kimura 658: }
659:
1.18 kimura 660: def unitweight(ExpMat,Vars,PolyListNum,OneMat,FLAG){
1.8 kimura 661:
662: RET=[]$
663:
664: ExpMatRowNum=size(ExpMat)[0]$
665: ExpMatColNum=size(ExpMat[0])[0]$
666: ExtMatColNum=ExpMatColNum+PolyListNum$
667:
668: ExtVars=reverse(Vars)$
669: for(I=0;I<PolyListNum;I++)
670: ExtVars=cons(uc(),ExtVars)$
671:
672: ExtVars=reverse(ExtVars)$
673:
1.18 kimura 674: NormMat0=newvect(ExpMatColNum)$
675: for(I=0;I<ExpMatColNum;I++)
676: NormMat0[I]=newvect(ExpMatColNum)$
1.8 kimura 677:
678: for(I=0;I<ExpMatColNum;I++)
679: for(J=I;J<ExpMatColNum;J++)
680: for(K=0;K<ExpMatRowNum;K++)
1.18 kimura 681: NormMat0[I][J]+=
1.8 kimura 682: ExpMat[K][I]*
683: ExpMat[K][J]$
684:
1.18 kimura 685: NormMat1=newvect(ExtMatColNum)$
686: for(I=0;I<ExtMatColNum;I++)
687: NormMat1[I]=newvect(ExtMatColNum)$
688:
689:
690: WorkMat=newvect(ExtMatColNum)$
691: for(I=0;I<ExtMatColNum;I++)
692: WorkMat[I]=newvect(ExtMatColNum)$
693:
694:
695: for(I=0;I<ExpMatColNum;I++)
696: for(J=I;J<ExpMatColNum;J++)
697: NormMat1[I][J]=NormMat0[I][J]$
698:
1.8 kimura 699: for(I=0;I<ExpMatColNum;I++)
700: for(J=0;J<PolyListNum;J++)
701: for(K=OneMat[J];K<OneMat[J+1];K++)
1.26 kimura 702: NormMat1[I][J+ExpMatColNum]-=ExpMat[K][I]$
1.8 kimura 703:
704: for(I=0;I<PolyListNum;I++)
1.18 kimura 705: NormMat1[I+ExpMatColNum][I+ExpMatColNum]=OneMat[I+1]-OneMat[I]$
706:
707: if(jacobi(ExtMatColNum,NormMat1,WorkMat)){
1.8 kimura 708:
1.30 ! kimura 709: Res=newvect(ExtMatColNum)$
! 710: for(I=0;I<ExtMatColNum;I++){
1.18 kimura 711: Res[I]=newvect(2)$
1.30 ! kimura 712: Res[I][0]=ExtVars[I]$
1.18 kimura 713: Res[I][1]=WorkMat[ExtMatColNum-1][I]$
714: }
1.8 kimura 715:
716: if(nonposdegchk(Res)){
717:
1.10 kimura 718: TMP1=makeret(Res,Vars,1)$
1.16 kimura 719:
1.26 kimura 720: if(FLAG==0){
721: TMP=roundret(TMP1[1])$
1.8 kimura 722:
1.26 kimura 723: RET=append(RET,wsort(TMP1[1],Vars,map(drint,TMP1[1]*1.0),1))$
1.18 kimura 724:
1.26 kimura 725: if(TMP!=[])
726: RET=append(RET,wsort(TMP1[1],Vars,TMP,2))$
727: }
1.27 kimura 728: else if(FLAG==1)
1.26 kimura 729: RET=append(RET,[[1,Vars,vtol(TMP1[1])]])$
1.4 kimura 730: }
1.18 kimura 731: }
732:
733: return [NormMat0,RET]$
1.6 kimura 734: }
735:
1.18 kimura 736: def weight(PolyList,Vars,FLAG){
1.6 kimura 737:
738: Vars0=vars(PolyList)$
739: Vars1=[]$
740: for(I=0;I<length(Vars);I++)
741: if(member(Vars[I],Vars0))
742: Vars1=cons(Vars[I],Vars1)$
743:
744: Vars=reverse(Vars1)$
745:
746: RET=[]$
747:
1.18 kimura 748: TMP=qcheck(PolyList,Vars,FLAG)$
1.6 kimura 749:
1.8 kimura 750: if(TMP!=[]){
751: RET=append(RET,TMP)$
1.18 kimura 752: return RET$
1.6 kimura 753: }
1.4 kimura 754:
1.6 kimura 755: dp_ord(2)$
1.4 kimura 756:
1.6 kimura 757: PolyListNum=length(PolyList)$
1.4 kimura 758:
1.18 kimura 759: OneMat=newvect(PolyListNum+1,[0])$
760: ExpMat=[]$
761: for(I=0;I<PolyListNum;I++){
762: for(Poly=dp_ptod(PolyList[I],Vars);
763: Poly!=0;Poly=dp_rest(Poly)){
764: ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
1.8 kimura 765: }
1.18 kimura 766: OneMat[I+1]=length(ExpMat)$
767: }
1.4 kimura 768:
1.18 kimura 769: ExpMat=reverse(ExpMat)$
770: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 771:
1.18 kimura 772: TMP=unitweight(ExpMat,Vars,PolyListNum,OneMat,FLAG)$
1.8 kimura 773:
1.18 kimura 774: RET=append(RET,TMP[1])$
1.4 kimura 775:
1.20 kimura 776: TMP0=leastsq(TMP[0],ExpMat,Vars,FLAG,3)$
777:
778: RET=append(RET,TMP0)$
1.4 kimura 779:
1.6 kimura 780: ExpMat=qsort(ExpMat,junban)$
1.8 kimura 781:
1.6 kimura 782: ExpMat2=[]$
783: for(I=0;I<size(ExpMat)[0];I++)
784: if(car(ExpMat2)!=ExpMat[I])
785: ExpMat2=cons(ExpMat[I],ExpMat2)$
1.4 kimura 786:
1.6 kimura 787: if(size(ExpMat)[0]!=length(ExpMat2)){
788: ExpMat=newvect(length(ExpMat2),ExpMat2)$
1.18 kimura 789: RET=append(RET,leastsq(0,ExpMat,Vars,FLAG,5))$
1.20 kimura 790: }
791: else{
792: TMP0=map(ltov,TMP0)$
793:
794: for(I=0;I<length(TMP0);I++)
1.24 kimura 795: if(TMP0[I][0]==3)
796: TMP0[I][0]=5$
797: else if(TMP0[I][0]==4)
798: TMP0[I][0]=6$
1.20 kimura 799:
800: TMP0=map(vtol,TMP0)$
801:
802: RET=append(RET,TMP0)$
1.4 kimura 803: }
804:
1.18 kimura 805: return RET$
1.4 kimura 806: }
807:
808: end$