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