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