![[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: }
123:
124: if(length(B)==0){
125: if(A)
126: return [Vars,1]$
127: else
128: return []$
129: }
130: else if(length(B)==1){
131:
132: C=fargs(B[0])$
133: D=vars(C)$
134: E=solve(C,D)$
135:
136: if(fop(B[0])==15)
137: return [Vars,E[0][1]+1]$
138: else if(fop(B[0])==11)
139: return [Vars,E[0][1]-1]$
1.25 ! kimura 140: else if(fop(B[0])==8)
! 141: return [Vars,E[0][1]]$
1.24 kimura 142: else
143: return []$
144: }
145: else{
146:
147: C=fargs(B[0])$
148: D=vars(C)$
149: E=solve(C,D)$
150:
151: C=fargs(B[1])$
152: D=vars(C)$
153: F=solve(C,D)$
154:
155: return [Vars,(E[0][1]+F[0][1])/2]$
156: }
157:
158: }
159:
160: def fixpointmain(F,Vars){
161:
162: RET=[]$
163: for(I=length(Vars)-1;I>=1;I--){
164:
165: for(H=[],J=0;J<I;J++)
166: H=cons(Vars[J],H)$
167:
168: G=interval2value(qe(ex(H,F)),Vars[I])$
169:
170: if(G==[])
171: return RET$
172: else
173: RET=cons(G,RET)$
174:
175: F=subf(F,G[0],G[1])$
176: }
177:
178: G=interval2value(simpl(F),Vars[0])$
179:
180: if(G==[])
181: return RET$
182: else
183: RET=cons(G,RET)$
184:
185: return RET$
186: }
187:
188:
189: def fixedpoint(A,FLAG){
190:
191: Vars=vars(A)$
192:
193: N=length(A)$
194:
195: if (FLAG==0)
196: for(F=@true,I=0;I < N; I++ ) { F = F @&& A[I] @> 0$ }
197: else if (FLAG==1)
198: for(F=@true,I=0;I < N; I++ ) { F = F @&& A[I] @< 0$ }
199:
200: return fixpointmain(F,Vars)$
201: }
202:
1.9 kimura 203: def nonzerovec(A){
204:
205: for(I=0;I<size(A)[0];I++)
206: if(A[I]!=0)
207: return 1$
208:
209: return 0$
210: }
211:
1.8 kimura 212: def junban(A,B){
213: return (A<B ? 1:(A>B ? -1:0))$
214: }
215:
216: def worder(A,B){
217: return (A[0]<B[0] ? 1:(A[0]>B[0] ? -1:0))$
218: }
219:
1.10 kimura 220: def bsort(A){
1.8 kimura 221:
1.10 kimura 222: K=size(A)[0]-1$
223: while(K>=0){
224: J=-1$
225: for(I=1;I<=K;I++)
226: if(A[I-1][0]<A[I][0]){
227: J=I-1$
228: X=A[J]$
229: A[J]=A[I]$
230: A[I]=X$
231: }
232: K=J$
233: }
234: return A$
235: }
236:
237: def perm(I,P,TMP){
238:
239: if(I>0){
240: TMP=perm(I-1,P,TMP)$
241: for(J=I-1;J>=0;J--){
242: T=P[I]$
243: P[I]=P[J]$
244: P[J]=T$
245: TMP=perm(I-1,P,TMP)$
246: T=P[I]$
247: P[I]=P[J]$
248: P[J]=T$
249: }
250:
251: return TMP$
252: }
253: else{
254: for(TMP0=[],K=0;K<size(P)[0];K++)
255: TMP0=cons(P[K],TMP0)$
256:
257: TMP=cons(TMP0,TMP)$
258: return TMP$
259: }
260: }
261:
262: def marge(A,B){
263:
264: RET=[]$
265: for(I=0;I<length(A);I++)
266: for(J=0;J<length(B);J++)
267: RET=cons(append(A[I],B[J]),RET)$
268:
269: return RET$
270: }
271:
1.18 kimura 272: def wsort(A,B,C,FLAG,ID){
1.10 kimura 273:
274: if(FLAG==0){
275: D=newvect(length(B))$
276: for(I=0;I<length(B);I++)
277: D[I]=[A[I],B[I],C[I]]$
278:
279: D=bsort(D)$
280: E=[]$
281: for(I=0;I<length(B);I++)
282: E=cons(D[I][1],E)$
283: E=reverse(E)$
284: F=[]$
285: for(I=0;I<length(B);I++)
286: F=cons(D[I][2],F)$
287: F=reverse(F)$
1.8 kimura 288:
1.18 kimura 289: return [[ID,E,F]]$
1.10 kimura 290: }
291: else{
292: D=newvect(length(B))$
293: for(I=0;I<length(B);I++)
294: D[I]=[A[I],B[I],C[I]]$
295:
296: D=qsort(D,worder)$
297: D0=[]$
298:
299: for(I=0,J=0,TMP=[],X=0;I<size(D)[0];I++){
300: if(X==D[I][0])
301: TMP=cons(cdr(D[I]),TMP)$
302: else{
303: D0=cons(TMP,D0)$
304: TMP=[]$
305: TMP=cons(cdr(D[I]),TMP)$
306: X=car(D[I])$
307: }
308: }
309: D0=cdr(reverse(cons(TMP,D0)))$
310: D0=map(ltov,D0)$
311: for(I=0,TMP=[[]];I<length(D0);I++){
312: TMP0=perm(length(D0[I])-1,D0[I],[])$
313: TMP=marge(TMP,TMP0)$
314: }
315:
316: RET=[]$
317: for(I=0;I<length(TMP);I++){
318: TMP0=[]$
319: TMP1=[]$
320: for(J=0;J<length(TMP[I]);J++){
321: TMP0=cons(TMP[I][J][0],TMP0)$
322: TMP1=cons(TMP[I][J][1],TMP1)$
323: }
324: TMP0=reverse(TMP0)$
325: TMP1=reverse(TMP1)$
326:
1.18 kimura 327: RET=cons([ID,TMP0,TMP1],RET)$
1.10 kimura 328: }
329:
330: return RET$
331: }
1.8 kimura 332: }
333:
1.4 kimura 334: def nonposdegchk(Res){
335:
336: for(I=0;I<length(Res);I++)
337: if(Res[I][1]<=0)
338: return 0$
339:
340: return 1$
341: }
342:
1.8 kimura 343: def getgcd(A,B){
344:
345: VarsNumA=length(A)$
346: VarsNumB=length(B)$
347:
348: C=newvect(VarsNumB,B)$
349:
350: for(I=0;I<VarsNumA;I++){
351:
352: for(J=0;J<VarsNumB;J++)
1.16 kimura 353: if(B[J]==A[I][0])
1.8 kimura 354: break$
355:
356: if(J<VarsNumB)
357: C[J]=A[I][1]$
358: }
359:
360: D=0$
361: for(I=0;I<VarsNumB;I++)
362: D=gcd(D,C[I])$
363:
364: if(D!=0){
365: C=C/D$
366: C=map(red,C)$
367: }
1.6 kimura 368:
1.8 kimura 369: for(L=1,D=0,I=0;I<VarsNumB;I++){
370: if(type(TMP=dn(C[I]))==1)
371: L=ilcm(L,TMP)$
372:
373: if(type(TMP=nm(C[I]))==1)
374: D=igcd(D,TMP)$
375: }
376:
377: C=C*L$
378: if(D!=0)
379: C=C/D$
1.6 kimura 380:
1.8 kimura 381: RET=[]$
382: for(I=0;I<VarsNumB;I++)
383: RET=cons([B[I],C[I]],RET)$
1.6 kimura 384:
1.8 kimura 385: return RET$
1.6 kimura 386: }
387:
1.10 kimura 388: def makeret(Res,Vars,FLAG){
1.4 kimura 389:
1.8 kimura 390: ResNum=length(Res)$
1.4 kimura 391: VarsNum=length(Vars)$
392:
1.17 kimura 393: ResVec=newvect(ResNum)$
1.4 kimura 394:
1.17 kimura 395: for(M=0,I=0;I<ResNum;I++){
396: if(member(Res[I][0],Vars)){
397: ResVec[I]=Res[I][1]$
1.16 kimura 398:
1.17 kimura 399: if(FLAG && type(ResVec[I])==1){
400: if(M==0)
401: M=ResVec[I]$
402: else
403: if(ResVec[I]<M)
404: M=ResVec[I]$
405: }
406: }
407: }
1.16 kimura 408:
1.8 kimura 409: if(M!=0)
410: ResVec=ResVec/M;
1.4 kimura 411:
1.17 kimura 412: RET=newvect(VarsNum,Vars)$
1.4 kimura 413:
1.8 kimura 414: for(I=0;I<ResNum;I++){
415: for(J=0;J<VarsNum;J++)
416: if(Vars[J]==Res[I][0])
417: break$
1.4 kimura 418:
1.8 kimura 419: if(J<VarsNum)
420: RET[J]=ResVec[I]$
421: }
422:
1.10 kimura 423:
424: for(J=0;J<length(Vars);J++)
425: RET=map(subst,RET,Vars[J],
426: strtov(rtostr(Vars[J])+"_deg"))$
427:
1.4 kimura 428: for(I=0;I<VarsNum;I++)
1.8 kimura 429: if(type(RET[I])!=1)
430: return [1,RET]$
1.4 kimura 431:
1.8 kimura 432: return [0,RET]$
1.4 kimura 433: }
434:
435: def roundret(V){
436:
1.8 kimura 437: VN=size(V)[0]$
1.4 kimura 438:
1.8 kimura 439: RET0=V$
1.13 kimura 440: for(I=1;I<1000;I++){
1.4 kimura 441: RET1=I*RET0$
442: for(J=0;J<VN;J++){
443: X=drint(RET1[J])$
1.19 noro 444: if(dabs(X-RET1[J])<ROUND_THRESHOLD)
1.4 kimura 445: RET1[J]=X$
446: else
447: break$
448: }
449: if(J==VN)
450: break$
451: }
452:
453: if(I==1000)
454: return []$
455: else
456: return RET1$
457: }
458:
459: def chkou(L,ExpMat,CHAGORD){
460:
1.8 kimura 461: for(P=1,I=0;I<L;I++){
1.4 kimura 462: Q=ExpMat[L][CHAGORD[I]]$
463: for(J=0;J<size(ExpMat[0])[0];J++){
464: ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
465: *ExpMat[L][CHAGORD[J]]-
466: Q*ExpMat[I][CHAGORD[J]])/P)$
467: }
468:
469: P=ExpMat[I][CHAGORD[I]]$
470: }
471:
472: for(J=0;J<size(ExpMat[0])[0];J++)
473: if(ExpMat[L][CHAGORD[J]]!=0)
474: break$
475:
476: if(J==size(ExpMat[0])[0])
477: return L$
478: else{
479: TMP=CHAGORD[L]$
480: CHAGORD[L]=CHAGORD[J]$
481: CHAGORD[J]=TMP$
482: return (L+1)$
483: }
484: }
485:
1.8 kimura 486: def qcheckmain(PolyList,Vars){
1.4 kimura 487:
488: RET=[]$
489: PolyListNum=length(PolyList)$
490: VarsNum=length(Vars)$
491:
492: ExpMat=newvect(VarsNum)$
493: CHAGORD=newvect(VarsNum)$
494: for(I=0;I<VarsNum;I++)
495: CHAGORD[I]=I$
496:
497: L=0$
498: for(I=0;I<PolyListNum;I++){
499: Poly=dp_ptod(PolyList[I],Vars)$
500: BASE0=dp_etov(dp_ht(Poly))$
501: Poly=dp_rest(Poly)$
502: for(;Poly!=0;Poly=dp_rest(Poly)){
503: ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
504: L=chkou(L,ExpMat,CHAGORD)$
1.8 kimura 505: if(L==VarsNum-1)
506: return [L,CHAGORD,ExpMat]$
1.4 kimura 507: }
508: }
509:
1.8 kimura 510: return [L,CHAGORD,ExpMat]$
1.4 kimura 511: }
512:
513: def inner(A,B){
514:
515: SUM=0$
516: for(I=0;I<size(A)[0];I++)
517: SUM+=A[I]*B[I]$
518:
519: return SUM$
520: }
521:
522: def checktd(PolyList,Vars,ResVars){
523:
524: PolyListNum=length(PolyList)$
525: VarsNum=length(Vars)$
526:
527: L=0$
528: for(I=0;I<PolyListNum;I++){
529: Poly=dp_ptod(PolyList[I],Vars)$
530: J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
531: Poly=dp_rest(Poly)$
532: for(;Poly!=0;Poly=dp_rest(Poly))
533: if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
534: return 0$
535: }
536:
537: return 1$
538: }
539:
1.24 kimura 540: def value2(Vars,Ans){
541:
542: N=length(Vars)$
543: Res=newvect(N)$
544: for(I=0;I<N;I++){
545: Res[I]=newvect(2)$
546: Res[I][0]=Vars[I]$
547: Res[I][1]=Ans[I]$
548: }
549:
550: Res=getgcd(Res,Vars)$
551:
552: if(nonposdegchk(Res)){
553: TMP1=makeret(Res,Vars,1)$
554: return vtol(TMP1[1])$
555: }
556: else
557: return []$
558: }
559:
1.10 kimura 560: def qcheck(PolyList,Vars,FLAG){
1.4 kimura 561:
1.10 kimura 562: RET=[]$
1.8 kimura 563: Res=qcheckmain(PolyList,Vars)$
1.4 kimura 564: VarsNum=length(Vars)$
565:
566: IndNum=Res[0]$
567: CHAGORD=Res[1]$
568: ExpMat=Res[2]$
569:
570: SolveList=[]$
571: for(I=0;I<IndNum;I++){
572: TMP=0$
573: for(J=0;J<VarsNum;J++)
574: TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$
575:
576: SolveList=cons(TMP,SolveList)$
577: }
578:
1.8 kimura 579: Rea=vars(SolveList)$
580:
1.4 kimura 581: VarsList=[]$
582: for(I=0;I<VarsNum;I++)
1.8 kimura 583: if(member(Vars[CHAGORD[I]],Rea))
584: VarsList=cons(Vars[CHAGORD[I]],VarsList)$
1.4 kimura 585:
586: Res=solve(reverse(SolveList),reverse(VarsList))$
1.8 kimura 587: Res=getgcd(Res,Rea)$
1.4 kimura 588:
589: if(nonposdegchk(Res)){
590:
1.10 kimura 591: ResVars=makeret(Res,Vars,0)$
1.4 kimura 592:
1.10 kimura 593: if(checktd(PolyList,Vars,ResVars[1])==1){
594: if(ResVars[0]==0){
595: RET=append(RET,wsort(ResVars[1],Vars,
1.18 kimura 596: ResVars[1],FLAG,0))$
597:
1.10 kimura 598: return RET$
599: }
600: else{
1.24 kimura 601:
602: TMP=vtol(ResVars[1])$
603:
604: /*
605: RET=append(RET,[[0,Vars,TMP]])$
606: */
607:
608: if((TMP0=fixedpoint(TMP,0))!=[]){
609:
610: for(I=0;I<length(TMP0);I++)
611: TMP=map(subst,TMP,TMP0[I][0],
612: TMP0[I][1])$
613:
614: TMP=value2(Vars,TMP)$
615:
616: if(TMP!=[])
617: RET=append(RET,wsort(TMP,Vars,
618: TMP,FLAG,1/10))$
619: }
620: else if((TMP0=fixedpoint(TMP,1))!=[]){
621:
622: for(I=0;I<length(TMP0);I++)
623: TMP=map(subst,TMP,TMP0[I][0],
624: TMP0[I][1])$
625:
626: TMP=value2(Vars,TMP)$
627:
628: if(TMP!=[])
629: RET=append(RET,wsort(TMP,Vars,
630: TMP,FLAG,1/10))$
631: }
632:
1.10 kimura 633: return RET$
634: }
1.4 kimura 635: }
636: else
1.8 kimura 637: return []$
1.4 kimura 638: }
639: else
1.8 kimura 640: return []$
1.4 kimura 641:
642: }
643:
1.18 kimura 644: def leastsq(NormMat,ExpMat,Vars,FLAG,ID){
1.8 kimura 645:
646: RET=[]$
1.4 kimura 647:
648: ExpMatRowNum=size(ExpMat)[0]$
649: ExpMatColNum=size(ExpMat[0])[0]$
650:
1.8 kimura 651: if(NormMat==0){
652: NormMat=newmat(ExpMatColNum,ExpMatColNum)$
653:
654: for(I=0;I<ExpMatColNum;I++)
655: for(J=I;J<ExpMatColNum;J++)
656: for(K=0;K<ExpMatRowNum;K++)
657: NormMat[I][J]+=
658: ExpMat[K][I]*ExpMat[K][J]$
659: }
1.4 kimura 660:
1.8 kimura 661: BVec=newvect(ExpMatColNum)$
1.4 kimura 662:
663: for(I=0;I<ExpMatColNum;I++)
664: for(J=0;J<ExpMatRowNum;J++)
1.8 kimura 665: BVec[I]+=ExpMat[J][I]$
1.4 kimura 666:
667: SolveList=[]$
668: for(I=0;I<ExpMatColNum;I++){
669: TMP=0$
1.8 kimura 670: for(J=0;J<I;J++)
671: TMP+=NormMat[J][I]*Vars[J]$
672:
673: for(J=I;J<ExpMatColNum;J++)
1.4 kimura 674: TMP+=NormMat[I][J]*Vars[J]$
675:
1.8 kimura 676: TMP-=BVec[I]$
1.4 kimura 677: SolveList=cons(TMP,SolveList)$
678: }
679:
680: Rea=vars(SolveList)$
1.8 kimura 681:
682: VarsList=[]$
683: for(I=0;I<length(Vars);I++)
684: if(member(Vars[I],Rea))
685: VarsList=cons(Vars[I],VarsList)$
686:
687: Res=solve(SolveList,VarsList)$
688: Res=getgcd(Res,Rea)$
1.4 kimura 689:
690: if(nonposdegchk(Res)){
1.16 kimura 691:
1.10 kimura 692: TMP1=makeret(Res,Vars,1)$
1.16 kimura 693:
1.8 kimura 694: if(TMP1[0]==0){
1.24 kimura 695:
1.14 kimura 696: TMP=roundret(TMP1[1])$
1.8 kimura 697:
1.10 kimura 698: RET=append(RET,wsort(TMP1[1],Vars,
1.18 kimura 699: map(drint,TMP1[1]*1.0),FLAG,ID))$
700:
701: if(TMP!=[])
702: RET=append(RET,wsort(TMP1[1],Vars,
703: TMP,FLAG,ID+1))$
1.8 kimura 704:
705: return RET$
706: }
707: else{
1.24 kimura 708:
709: TMP=vtol(TMP1[1])$
710:
711: /*
1.23 kimura 712: RET=append(RET,[[ID,Vars,vtol(TMP1[1])]])$
1.24 kimura 713: */
714:
715: if((TMP0=fixedpoint(TMP1[1],0))!=[]){
716:
717: for(I=0;I<length(TMP0);I++)
718: TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
719:
720: TMP=value2(Vars,TMP)$
721:
722: if(TMP!=[])
723: RET=append(RET,
724: wsort(TMP,Vars,TMP,FLAG,ID+1/10))$
725:
726: }
727: else if((TMP0=fixedpoint(TMP1[1],1))!=[]){
728:
729: for(I=0;I<length(TMP0);I++)
730: TMP=map(subst,TMP,TMP0[I][0],TMP0[I][1])$
731:
732: TMP=value2(Vars,TMP)$
733:
734: if(TMP!=[])
735: RET=append(RET,
736: wsort(TMP,Vars,TMP,FLAG,ID+1/10))$
737: }
738:
1.8 kimura 739: return RET$
740: }
741: }
742: else
743: return RET$
744:
745: }
746:
1.18 kimura 747: def unitweight(ExpMat,Vars,PolyListNum,OneMat,FLAG){
1.8 kimura 748:
749: RET=[]$
750:
751: ExpMatRowNum=size(ExpMat)[0]$
752: ExpMatColNum=size(ExpMat[0])[0]$
753: ExtMatColNum=ExpMatColNum+PolyListNum$
754:
755: ExtVars=reverse(Vars)$
756: for(I=0;I<PolyListNum;I++)
757: ExtVars=cons(uc(),ExtVars)$
758:
759: ExtVars=reverse(ExtVars)$
760:
1.18 kimura 761: NormMat0=newvect(ExpMatColNum)$
762: for(I=0;I<ExpMatColNum;I++)
763: NormMat0[I]=newvect(ExpMatColNum)$
1.8 kimura 764:
765: for(I=0;I<ExpMatColNum;I++)
766: for(J=I;J<ExpMatColNum;J++)
767: for(K=0;K<ExpMatRowNum;K++)
1.18 kimura 768: NormMat0[I][J]+=
1.8 kimura 769: ExpMat[K][I]*
770: ExpMat[K][J]$
771:
1.18 kimura 772: NormMat1=newvect(ExtMatColNum)$
773: for(I=0;I<ExtMatColNum;I++)
774: NormMat1[I]=newvect(ExtMatColNum)$
775:
776:
777: WorkMat=newvect(ExtMatColNum)$
778: for(I=0;I<ExtMatColNum;I++)
779: WorkMat[I]=newvect(ExtMatColNum)$
780:
781:
782: for(I=0;I<ExpMatColNum;I++)
783: for(J=I;J<ExpMatColNum;J++)
784: NormMat1[I][J]=NormMat0[I][J]$
785:
1.8 kimura 786: for(I=0;I<ExpMatColNum;I++)
787: for(J=0;J<PolyListNum;J++)
788: for(K=OneMat[J];K<OneMat[J+1];K++)
1.18 kimura 789: NormMat1[I][J+ExpMatColNum]-=
1.8 kimura 790: ExpMat[K][I]$
791:
792: for(I=0;I<PolyListNum;I++)
1.18 kimura 793: NormMat1[I+ExpMatColNum][I+ExpMatColNum]=OneMat[I+1]-OneMat[I]$
794:
795: if(jacobi(ExtMatColNum,NormMat1,WorkMat)){
1.8 kimura 796:
1.18 kimura 797: Res=newvect(ExpMatColNum)$
1.8 kimura 798: for(I=0;I<ExpMatColNum;I++){
1.18 kimura 799: Res[I]=newvect(2)$
800: Res[I][0]=Vars[I]$
801: Res[I][1]=WorkMat[ExtMatColNum-1][I]$
802: }
1.8 kimura 803:
804: if(nonposdegchk(Res)){
805:
1.10 kimura 806: TMP1=makeret(Res,Vars,1)$
1.16 kimura 807:
1.18 kimura 808: TMP=roundret(TMP1[1])$
1.8 kimura 809:
1.18 kimura 810: RET=append(RET,wsort(TMP1[1],Vars,
811: map(drint,TMP1[1]*1.0),FLAG,1))$
812:
813: if(TMP!=[])
1.10 kimura 814: RET=append(RET,wsort(TMP1[1],Vars,
1.18 kimura 815: TMP,FLAG,2))$
1.4 kimura 816: }
1.8 kimura 817:
1.18 kimura 818: }
819:
820: return [NormMat0,RET]$
1.6 kimura 821: }
822:
1.18 kimura 823: def weight(PolyList,Vars,FLAG){
1.6 kimura 824:
825: Vars0=vars(PolyList)$
826: Vars1=[]$
827: for(I=0;I<length(Vars);I++)
828: if(member(Vars[I],Vars0))
829: Vars1=cons(Vars[I],Vars1)$
830:
831: Vars=reverse(Vars1)$
832:
833: RET=[]$
834:
1.18 kimura 835: TMP=qcheck(PolyList,Vars,FLAG)$
1.6 kimura 836:
1.8 kimura 837: if(TMP!=[]){
838: RET=append(RET,TMP)$
1.18 kimura 839: return RET$
1.6 kimura 840: }
1.4 kimura 841:
1.6 kimura 842: dp_ord(2)$
1.4 kimura 843:
1.6 kimura 844: PolyListNum=length(PolyList)$
1.4 kimura 845:
1.18 kimura 846: OneMat=newvect(PolyListNum+1,[0])$
847: ExpMat=[]$
848: for(I=0;I<PolyListNum;I++){
849: for(Poly=dp_ptod(PolyList[I],Vars);
850: Poly!=0;Poly=dp_rest(Poly)){
851: ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
1.8 kimura 852: }
1.18 kimura 853: OneMat[I+1]=length(ExpMat)$
854: }
1.4 kimura 855:
1.18 kimura 856: ExpMat=reverse(ExpMat)$
857: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 858:
1.18 kimura 859: TMP=unitweight(ExpMat,Vars,PolyListNum,OneMat,FLAG)$
1.8 kimura 860:
1.18 kimura 861: RET=append(RET,TMP[1])$
1.4 kimura 862:
1.20 kimura 863: TMP0=leastsq(TMP[0],ExpMat,Vars,FLAG,3)$
864:
865: RET=append(RET,TMP0)$
1.4 kimura 866:
1.6 kimura 867: ExpMat=qsort(ExpMat,junban)$
1.8 kimura 868:
1.6 kimura 869: ExpMat2=[]$
870: for(I=0;I<size(ExpMat)[0];I++)
871: if(car(ExpMat2)!=ExpMat[I])
872: ExpMat2=cons(ExpMat[I],ExpMat2)$
1.4 kimura 873:
1.6 kimura 874: if(size(ExpMat)[0]!=length(ExpMat2)){
875: ExpMat=newvect(length(ExpMat2),ExpMat2)$
1.18 kimura 876: RET=append(RET,leastsq(0,ExpMat,Vars,FLAG,5))$
1.20 kimura 877: }
878: else{
879: TMP0=map(ltov,TMP0)$
880:
881: for(I=0;I<length(TMP0);I++)
1.24 kimura 882: if(TMP0[I][0]==3)
883: TMP0[I][0]=5$
884: else if(TMP0[I][0]==4)
885: TMP0[I][0]=6$
1.20 kimura 886:
887: TMP0=map(vtol,TMP0)$
888:
889: RET=append(RET,TMP0)$
1.4 kimura 890: }
891:
1.18 kimura 892: return RET$
1.4 kimura 893: }
894:
895: end$