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