![[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.9 kimura 115: def nonzerovec(A){
116:
117: for(I=0;I<size(A)[0];I++)
118: if(A[I]!=0)
119: return 1$
120:
121: return 0$
122: }
123:
1.8 kimura 124: def junban(A,B){
125: return (A<B ? 1:(A>B ? -1:0))$
126: }
127:
128: def worder(A,B){
129: return (A[0]<B[0] ? 1:(A[0]>B[0] ? -1:0))$
130: }
131:
1.10 kimura 132: def bsort(A){
1.8 kimura 133:
1.10 kimura 134: K=size(A)[0]-1$
135: while(K>=0){
136: J=-1$
137: for(I=1;I<=K;I++)
138: if(A[I-1][0]<A[I][0]){
139: J=I-1$
140: X=A[J]$
141: A[J]=A[I]$
142: A[I]=X$
143: }
144: K=J$
145: }
146: return A$
147: }
148:
149: def perm(I,P,TMP){
150:
151: if(I>0){
152: TMP=perm(I-1,P,TMP)$
153: for(J=I-1;J>=0;J--){
154: T=P[I]$
155: P[I]=P[J]$
156: P[J]=T$
157: TMP=perm(I-1,P,TMP)$
158: T=P[I]$
159: P[I]=P[J]$
160: P[J]=T$
161: }
162:
163: return TMP$
164: }
165: else{
166: for(TMP0=[],K=0;K<size(P)[0];K++)
167: TMP0=cons(P[K],TMP0)$
168:
169: TMP=cons(TMP0,TMP)$
170: return TMP$
171: }
172: }
173:
174: def marge(A,B){
175:
176: RET=[]$
177: for(I=0;I<length(A);I++)
178: for(J=0;J<length(B);J++)
179: RET=cons(append(A[I],B[J]),RET)$
180:
181: return RET$
182: }
183:
1.18 kimura 184: def wsort(A,B,C,FLAG,ID){
1.10 kimura 185:
186: if(FLAG==0){
187: D=newvect(length(B))$
188: for(I=0;I<length(B);I++)
189: D[I]=[A[I],B[I],C[I]]$
190:
191: D=bsort(D)$
192: E=[]$
193: for(I=0;I<length(B);I++)
194: E=cons(D[I][1],E)$
195: E=reverse(E)$
196: F=[]$
197: for(I=0;I<length(B);I++)
198: F=cons(D[I][2],F)$
199: F=reverse(F)$
1.8 kimura 200:
1.18 kimura 201: return [[ID,E,F]]$
1.10 kimura 202: }
203: else{
204: D=newvect(length(B))$
205: for(I=0;I<length(B);I++)
206: D[I]=[A[I],B[I],C[I]]$
207:
208: D=qsort(D,worder)$
209: D0=[]$
210:
211: for(I=0,J=0,TMP=[],X=0;I<size(D)[0];I++){
212: if(X==D[I][0])
213: TMP=cons(cdr(D[I]),TMP)$
214: else{
215: D0=cons(TMP,D0)$
216: TMP=[]$
217: TMP=cons(cdr(D[I]),TMP)$
218: X=car(D[I])$
219: }
220: }
221: D0=cdr(reverse(cons(TMP,D0)))$
222: D0=map(ltov,D0)$
223: for(I=0,TMP=[[]];I<length(D0);I++){
224: TMP0=perm(length(D0[I])-1,D0[I],[])$
225: TMP=marge(TMP,TMP0)$
226: }
227:
228: RET=[]$
229: for(I=0;I<length(TMP);I++){
230: TMP0=[]$
231: TMP1=[]$
232: for(J=0;J<length(TMP[I]);J++){
233: TMP0=cons(TMP[I][J][0],TMP0)$
234: TMP1=cons(TMP[I][J][1],TMP1)$
235: }
236: TMP0=reverse(TMP0)$
237: TMP1=reverse(TMP1)$
238:
1.18 kimura 239: RET=cons([ID,TMP0,TMP1],RET)$
1.10 kimura 240: }
241:
242: return RET$
243: }
1.8 kimura 244: }
245:
1.4 kimura 246: def nonposdegchk(Res){
247:
248: for(I=0;I<length(Res);I++)
249: if(Res[I][1]<=0)
250: return 0$
251:
252: return 1$
253: }
254:
1.8 kimura 255: def getgcd(A,B){
256:
257: VarsNumA=length(A)$
258: VarsNumB=length(B)$
259:
260: C=newvect(VarsNumB,B)$
261:
262: for(I=0;I<VarsNumA;I++){
263:
264: for(J=0;J<VarsNumB;J++)
1.16 kimura 265: if(B[J]==A[I][0])
1.8 kimura 266: break$
267:
268: if(J<VarsNumB)
269: C[J]=A[I][1]$
270: }
271:
272: D=0$
273: for(I=0;I<VarsNumB;I++)
274: D=gcd(D,C[I])$
275:
276: if(D!=0){
277: C=C/D$
278: C=map(red,C)$
279: }
1.6 kimura 280:
1.8 kimura 281: for(L=1,D=0,I=0;I<VarsNumB;I++){
282: if(type(TMP=dn(C[I]))==1)
283: L=ilcm(L,TMP)$
284:
285: if(type(TMP=nm(C[I]))==1)
286: D=igcd(D,TMP)$
287: }
288:
289: C=C*L$
290: if(D!=0)
291: C=C/D$
1.6 kimura 292:
1.8 kimura 293: RET=[]$
294: for(I=0;I<VarsNumB;I++)
295: RET=cons([B[I],C[I]],RET)$
1.6 kimura 296:
1.8 kimura 297: return RET$
1.6 kimura 298: }
299:
1.10 kimura 300: def makeret(Res,Vars,FLAG){
1.4 kimura 301:
1.8 kimura 302: ResNum=length(Res)$
1.4 kimura 303: VarsNum=length(Vars)$
304:
1.17 kimura 305: ResVec=newvect(ResNum)$
1.4 kimura 306:
1.17 kimura 307: for(M=0,I=0;I<ResNum;I++){
308: if(member(Res[I][0],Vars)){
309: ResVec[I]=Res[I][1]$
1.16 kimura 310:
1.17 kimura 311: if(FLAG && type(ResVec[I])==1){
312: if(M==0)
313: M=ResVec[I]$
314: else
315: if(ResVec[I]<M)
316: M=ResVec[I]$
317: }
318: }
319: }
1.16 kimura 320:
1.8 kimura 321: if(M!=0)
322: ResVec=ResVec/M;
1.4 kimura 323:
1.17 kimura 324: RET=newvect(VarsNum,Vars)$
1.4 kimura 325:
1.8 kimura 326: for(I=0;I<ResNum;I++){
327: for(J=0;J<VarsNum;J++)
328: if(Vars[J]==Res[I][0])
329: break$
1.4 kimura 330:
1.8 kimura 331: if(J<VarsNum)
332: RET[J]=ResVec[I]$
333: }
334:
1.10 kimura 335:
336: for(J=0;J<length(Vars);J++)
337: RET=map(subst,RET,Vars[J],
338: strtov(rtostr(Vars[J])+"_deg"))$
339:
1.4 kimura 340: for(I=0;I<VarsNum;I++)
1.8 kimura 341: if(type(RET[I])!=1)
342: return [1,RET]$
1.4 kimura 343:
1.8 kimura 344: return [0,RET]$
1.4 kimura 345: }
346:
347: def roundret(V){
348:
1.8 kimura 349: VN=size(V)[0]$
1.4 kimura 350:
1.8 kimura 351: RET0=V$
1.13 kimura 352: for(I=1;I<1000;I++){
1.4 kimura 353: RET1=I*RET0$
354: for(J=0;J<VN;J++){
355: X=drint(RET1[J])$
1.19 noro 356: if(dabs(X-RET1[J])<ROUND_THRESHOLD)
1.4 kimura 357: RET1[J]=X$
358: else
359: break$
360: }
361: if(J==VN)
362: break$
363: }
364:
365: if(I==1000)
366: return []$
367: else
368: return RET1$
369: }
370:
371: def chkou(L,ExpMat,CHAGORD){
372:
1.8 kimura 373: for(P=1,I=0;I<L;I++){
1.4 kimura 374: Q=ExpMat[L][CHAGORD[I]]$
375: for(J=0;J<size(ExpMat[0])[0];J++){
376: ExpMat[L][CHAGORD[J]]=red((ExpMat[I][CHAGORD[I]]
377: *ExpMat[L][CHAGORD[J]]-
378: Q*ExpMat[I][CHAGORD[J]])/P)$
379: }
380:
381: P=ExpMat[I][CHAGORD[I]]$
382: }
383:
384: for(J=0;J<size(ExpMat[0])[0];J++)
385: if(ExpMat[L][CHAGORD[J]]!=0)
386: break$
387:
388: if(J==size(ExpMat[0])[0])
389: return L$
390: else{
391: TMP=CHAGORD[L]$
392: CHAGORD[L]=CHAGORD[J]$
393: CHAGORD[J]=TMP$
394: return (L+1)$
395: }
396: }
397:
1.8 kimura 398: def qcheckmain(PolyList,Vars){
1.4 kimura 399:
400: RET=[]$
401: PolyListNum=length(PolyList)$
402: VarsNum=length(Vars)$
403:
404: ExpMat=newvect(VarsNum)$
405: CHAGORD=newvect(VarsNum)$
406: for(I=0;I<VarsNum;I++)
407: CHAGORD[I]=I$
408:
409: L=0$
410: for(I=0;I<PolyListNum;I++){
411: Poly=dp_ptod(PolyList[I],Vars)$
412: BASE0=dp_etov(dp_ht(Poly))$
413: Poly=dp_rest(Poly)$
414: for(;Poly!=0;Poly=dp_rest(Poly)){
415: ExpMat[L]=dp_etov(dp_ht(Poly))-BASE0$
416: L=chkou(L,ExpMat,CHAGORD)$
1.8 kimura 417: if(L==VarsNum-1)
418: return [L,CHAGORD,ExpMat]$
1.4 kimura 419: }
420: }
421:
1.8 kimura 422: return [L,CHAGORD,ExpMat]$
1.4 kimura 423: }
424:
425: def inner(A,B){
426:
427: SUM=0$
428: for(I=0;I<size(A)[0];I++)
429: SUM+=A[I]*B[I]$
430:
431: return SUM$
432: }
433:
434: def checktd(PolyList,Vars,ResVars){
435:
436: PolyListNum=length(PolyList)$
437: VarsNum=length(Vars)$
438:
439: L=0$
440: for(I=0;I<PolyListNum;I++){
441: Poly=dp_ptod(PolyList[I],Vars)$
442: J0=inner(dp_etov(dp_ht(Poly)),ResVars)$
443: Poly=dp_rest(Poly)$
444: for(;Poly!=0;Poly=dp_rest(Poly))
445: if(J0!=inner(dp_etov(dp_ht(Poly)),ResVars))
446: return 0$
447: }
448:
449: return 1$
450: }
451:
1.10 kimura 452: def qcheck(PolyList,Vars,FLAG){
1.4 kimura 453:
1.10 kimura 454: RET=[]$
1.8 kimura 455: Res=qcheckmain(PolyList,Vars)$
1.4 kimura 456: VarsNum=length(Vars)$
457:
458: IndNum=Res[0]$
459: CHAGORD=Res[1]$
460: ExpMat=Res[2]$
461:
462: SolveList=[]$
463: for(I=0;I<IndNum;I++){
464: TMP=0$
465: for(J=0;J<VarsNum;J++)
466: TMP+=ExpMat[I][CHAGORD[J]]*Vars[CHAGORD[J]]$
467:
468: SolveList=cons(TMP,SolveList)$
469: }
470:
1.8 kimura 471: Rea=vars(SolveList)$
472:
1.4 kimura 473: VarsList=[]$
474: for(I=0;I<VarsNum;I++)
1.8 kimura 475: if(member(Vars[CHAGORD[I]],Rea))
476: VarsList=cons(Vars[CHAGORD[I]],VarsList)$
1.4 kimura 477:
478: Res=solve(reverse(SolveList),reverse(VarsList))$
1.8 kimura 479: Res=getgcd(Res,Rea)$
1.4 kimura 480:
481: if(nonposdegchk(Res)){
482:
1.10 kimura 483: ResVars=makeret(Res,Vars,0)$
1.4 kimura 484:
1.10 kimura 485: if(checktd(PolyList,Vars,ResVars[1])==1){
486: if(ResVars[0]==0){
487: RET=append(RET,wsort(ResVars[1],Vars,
1.18 kimura 488: ResVars[1],FLAG,0))$
489:
1.10 kimura 490: return RET$
491: }
492: else{
1.18 kimura 493: RET=append(RET,[[0,Vars,vtol(ResVars[1])]])$
1.10 kimura 494: return RET$
495: }
1.4 kimura 496: }
497: else
1.8 kimura 498: return []$
1.4 kimura 499: }
500: else
1.8 kimura 501: return []$
1.4 kimura 502:
503: }
504:
1.18 kimura 505: def leastsq(NormMat,ExpMat,Vars,FLAG,ID){
1.8 kimura 506:
507: RET=[]$
1.4 kimura 508:
509: ExpMatRowNum=size(ExpMat)[0]$
510: ExpMatColNum=size(ExpMat[0])[0]$
511:
1.8 kimura 512: if(NormMat==0){
513: NormMat=newmat(ExpMatColNum,ExpMatColNum)$
514:
515: for(I=0;I<ExpMatColNum;I++)
516: for(J=I;J<ExpMatColNum;J++)
517: for(K=0;K<ExpMatRowNum;K++)
518: NormMat[I][J]+=
519: ExpMat[K][I]*ExpMat[K][J]$
520: }
1.4 kimura 521:
1.8 kimura 522: BVec=newvect(ExpMatColNum)$
1.4 kimura 523:
524: for(I=0;I<ExpMatColNum;I++)
525: for(J=0;J<ExpMatRowNum;J++)
1.8 kimura 526: BVec[I]+=ExpMat[J][I]$
1.4 kimura 527:
528: SolveList=[]$
529: for(I=0;I<ExpMatColNum;I++){
530: TMP=0$
1.8 kimura 531: for(J=0;J<I;J++)
532: TMP+=NormMat[J][I]*Vars[J]$
533:
534: for(J=I;J<ExpMatColNum;J++)
1.4 kimura 535: TMP+=NormMat[I][J]*Vars[J]$
536:
1.8 kimura 537: TMP-=BVec[I]$
1.4 kimura 538: SolveList=cons(TMP,SolveList)$
539: }
540:
541: Rea=vars(SolveList)$
1.8 kimura 542:
543: VarsList=[]$
544: for(I=0;I<length(Vars);I++)
545: if(member(Vars[I],Rea))
546: VarsList=cons(Vars[I],VarsList)$
547:
548: Res=solve(SolveList,VarsList)$
549: Res=getgcd(Res,Rea)$
1.4 kimura 550:
551: if(nonposdegchk(Res)){
1.16 kimura 552:
1.10 kimura 553: TMP1=makeret(Res,Vars,1)$
1.16 kimura 554:
1.8 kimura 555: if(TMP1[0]==0){
1.14 kimura 556: TMP=roundret(TMP1[1])$
1.8 kimura 557:
1.10 kimura 558: RET=append(RET,wsort(TMP1[1],Vars,
1.18 kimura 559: map(drint,TMP1[1]*1.0),FLAG,ID))$
560:
561: if(TMP!=[])
562: RET=append(RET,wsort(TMP1[1],Vars,
563: TMP,FLAG,ID+1))$
1.8 kimura 564:
565: return RET$
566: }
567: else{
1.18 kimura 568: RET=append(RET,[[ID,Vars,vtol(TMP1[1]*1.0)]])$
1.8 kimura 569: return RET$
570: }
571: }
572: else
573: return RET$
574:
575: }
576:
1.18 kimura 577: def unitweight(ExpMat,Vars,PolyListNum,OneMat,FLAG){
1.8 kimura 578:
579: RET=[]$
580:
581: ExpMatRowNum=size(ExpMat)[0]$
582: ExpMatColNum=size(ExpMat[0])[0]$
583: ExtMatColNum=ExpMatColNum+PolyListNum$
584:
585: ExtVars=reverse(Vars)$
586: for(I=0;I<PolyListNum;I++)
587: ExtVars=cons(uc(),ExtVars)$
588:
589: ExtVars=reverse(ExtVars)$
590:
1.18 kimura 591: NormMat0=newvect(ExpMatColNum)$
592: for(I=0;I<ExpMatColNum;I++)
593: NormMat0[I]=newvect(ExpMatColNum)$
1.8 kimura 594:
595: for(I=0;I<ExpMatColNum;I++)
596: for(J=I;J<ExpMatColNum;J++)
597: for(K=0;K<ExpMatRowNum;K++)
1.18 kimura 598: NormMat0[I][J]+=
1.8 kimura 599: ExpMat[K][I]*
600: ExpMat[K][J]$
601:
1.18 kimura 602: NormMat1=newvect(ExtMatColNum)$
603: for(I=0;I<ExtMatColNum;I++)
604: NormMat1[I]=newvect(ExtMatColNum)$
605:
606:
607: WorkMat=newvect(ExtMatColNum)$
608: for(I=0;I<ExtMatColNum;I++)
609: WorkMat[I]=newvect(ExtMatColNum)$
610:
611:
612: for(I=0;I<ExpMatColNum;I++)
613: for(J=I;J<ExpMatColNum;J++)
614: NormMat1[I][J]=NormMat0[I][J]$
615:
1.8 kimura 616: for(I=0;I<ExpMatColNum;I++)
617: for(J=0;J<PolyListNum;J++)
618: for(K=OneMat[J];K<OneMat[J+1];K++)
1.18 kimura 619: NormMat1[I][J+ExpMatColNum]-=
1.8 kimura 620: ExpMat[K][I]$
621:
622: for(I=0;I<PolyListNum;I++)
1.18 kimura 623: NormMat1[I+ExpMatColNum][I+ExpMatColNum]=OneMat[I+1]-OneMat[I]$
624:
625: if(jacobi(ExtMatColNum,NormMat1,WorkMat)){
1.8 kimura 626:
1.18 kimura 627: Res=newvect(ExpMatColNum)$
1.8 kimura 628: for(I=0;I<ExpMatColNum;I++){
1.18 kimura 629: Res[I]=newvect(2)$
630: Res[I][0]=Vars[I]$
631: Res[I][1]=WorkMat[ExtMatColNum-1][I]$
632: }
1.8 kimura 633:
634: if(nonposdegchk(Res)){
635:
1.10 kimura 636: TMP1=makeret(Res,Vars,1)$
1.16 kimura 637:
1.18 kimura 638: TMP=roundret(TMP1[1])$
1.8 kimura 639:
1.18 kimura 640: RET=append(RET,wsort(TMP1[1],Vars,
641: map(drint,TMP1[1]*1.0),FLAG,1))$
642:
643: if(TMP!=[])
1.10 kimura 644: RET=append(RET,wsort(TMP1[1],Vars,
1.18 kimura 645: TMP,FLAG,2))$
1.4 kimura 646: }
1.8 kimura 647:
1.18 kimura 648: }
649:
650: return [NormMat0,RET]$
1.6 kimura 651: }
652:
1.18 kimura 653: def weight(PolyList,Vars,FLAG){
1.6 kimura 654:
655: Vars0=vars(PolyList)$
656: Vars1=[]$
657: for(I=0;I<length(Vars);I++)
658: if(member(Vars[I],Vars0))
659: Vars1=cons(Vars[I],Vars1)$
660:
661: Vars=reverse(Vars1)$
662:
663: RET=[]$
664:
1.18 kimura 665: TMP=qcheck(PolyList,Vars,FLAG)$
1.6 kimura 666:
1.8 kimura 667: if(TMP!=[]){
668: RET=append(RET,TMP)$
1.18 kimura 669: return RET$
1.6 kimura 670: }
1.4 kimura 671:
1.6 kimura 672: dp_ord(2)$
1.4 kimura 673:
1.6 kimura 674: PolyListNum=length(PolyList)$
1.4 kimura 675:
1.18 kimura 676: OneMat=newvect(PolyListNum+1,[0])$
677: ExpMat=[]$
678: for(I=0;I<PolyListNum;I++){
679: for(Poly=dp_ptod(PolyList[I],Vars);
680: Poly!=0;Poly=dp_rest(Poly)){
681: ExpMat=cons(dp_etov(dp_ht(Poly)),ExpMat)$
1.8 kimura 682: }
1.18 kimura 683: OneMat[I+1]=length(ExpMat)$
684: }
1.4 kimura 685:
1.18 kimura 686: ExpMat=reverse(ExpMat)$
687: ExpMat=newvect(length(ExpMat),ExpMat)$
1.4 kimura 688:
1.18 kimura 689: TMP=unitweight(ExpMat,Vars,PolyListNum,OneMat,FLAG)$
1.8 kimura 690:
1.18 kimura 691: RET=append(RET,TMP[1])$
1.4 kimura 692:
1.20 kimura 693: TMP0=leastsq(TMP[0],ExpMat,Vars,FLAG,3)$
694:
695: RET=append(RET,TMP0)$
1.4 kimura 696:
1.6 kimura 697: ExpMat=qsort(ExpMat,junban)$
1.8 kimura 698:
1.6 kimura 699: ExpMat2=[]$
700: for(I=0;I<size(ExpMat)[0];I++)
701: if(car(ExpMat2)!=ExpMat[I])
702: ExpMat2=cons(ExpMat[I],ExpMat2)$
1.4 kimura 703:
1.6 kimura 704: if(size(ExpMat)[0]!=length(ExpMat2)){
705: ExpMat=newvect(length(ExpMat2),ExpMat2)$
1.18 kimura 706: RET=append(RET,leastsq(0,ExpMat,Vars,FLAG,5))$
1.20 kimura 707: }
708: else{
709: TMP0=map(ltov,TMP0)$
710:
711: for(I=0;I<length(TMP0);I++)
712: if(TMP0[I][0]==3)
713: TMP0[I][0]=5$
714: else if(TMP0[I][0]==4)
715: TMP0[I][0]=6$
716:
717: TMP0=map(vtol,TMP0)$
718:
719: RET=append(RET,TMP0)$
1.4 kimura 720: }
721:
1.18 kimura 722: return RET$
1.4 kimura 723: }
724:
725: end$