=================================================================== RCS file: /home/cvs/OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v retrieving revision 1.50 retrieving revision 1.69 diff -u -p -r1.50 -r1.69 --- OpenXM/src/asir-contrib/packages/src/os_muldif.rr 2019/06/27 02:53:26 1.50 +++ OpenXM/src/asir-contrib/packages/src/os_muldif.rr 2020/05/17 23:15:26 1.69 @@ -1,12 +1,12 @@ -/* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.49 2019/05/23 01:47:53 takayama Exp $ */ -/* The latest version will be at ftp://akagi.ms.u-tokyo.ac.jp/pub/math/muldif +/* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.68 2020/05/16 05:40:32 takayama Exp $ */ +/* The latest version will be at https://www.ms.u-tokyo.ac.jp/~oshima/index-j.html scp os_muldif.[dp]* ${USER}@lemon.math.kobe-u.ac.jp:/home/web/OpenXM/Current/doc/other-docs */ #define USEMODULE 1 /* #undef USEMODULE */ /* os_muldif.rr (Library for Risa/Asir) - * Toshio Oshima (Nov. 2007 - Feb. 2019) + * Toshio Oshima (Nov. 2007 - May. 2020) * * For polynomials and differential operators with coefficients * in rational funtions (See os_muldif.pdf) @@ -59,6 +59,9 @@ localf countin$ localf mycoef$ localf mydiff$ localf myediff$ +localf mypdiff$ +localf pTaylor$ +localf pwTaylor$ localf m2l$ localf m2ll$ localf mydeg$ @@ -116,6 +119,7 @@ localf ladd$ localf lchange$ localf llsize$ localf llbase$ +localf llget$ localf lsort$ localf rsort$ localf lpair$ @@ -139,6 +143,7 @@ localf mpower$ localf mrot$ localf texlen$ localf isdif$ +localf isfctr$ localf fctrtos$ localf texlim$ localf fmult$ @@ -147,6 +152,8 @@ localf getel$ localf ptol$ localf rmul$ localf mtransbys$ +localf trcolor$ +localf mcolor$ localf drawopt$ localf execdraw$ localf execproc$ @@ -171,6 +178,7 @@ localf myasin$ localf myacos$ localf myatan$ localf mylog$ +localf nlog$ localf mypow$ localf scale$ localf arg$ @@ -260,6 +268,7 @@ localf expat$ localf polbyroot$ localf polbyvalue$ localf pcoef$ +localf pmaj$ localf prehombf$ localf prehombfold$ localf sub3e$ @@ -297,6 +306,7 @@ localf iscoef$ localf iscombox$ localf sproot$ localf spgen$ +localf spbasic$ localf chkspt$ localf cterm$ localf terms$ @@ -337,6 +347,7 @@ localf divmattex$ localf dviout0$ localf myhelp$ localf isMs$ +localf getline$ localf showbyshell$ localf readcsv$ localf tocsv$ @@ -352,6 +363,7 @@ localf texsp$ localf getbygrs$ localf mcop$ localf shiftop$ +localf shiftPfaff; localf conf1sp$ localf confexp$ localf confspt$ @@ -391,10 +403,12 @@ localf primroot$ localf varargs$ localf ptype$ localf pfargs$ +localf regress$ localf average$ localf tobig$ localf sint$ localf frac2n$ +localf openGlib$ localf xyproc$ localf xypos$ localf xyput$ @@ -415,6 +429,8 @@ localf periodicf$ localf cmpf$ localf areabezier$ localf saveproc$ +localf xyplot$ +localf xyaxis$ localf xygraph$ localf xy2graph$ localf addIL$ @@ -470,12 +486,15 @@ extern Canvas$ extern ID_PLOT$ extern Rand$ extern LQS$ -extern SV=SVORG$ +extern SV=SVORG$ #endif static S_Fc,S_Dc,S_Ic,S_Ec,S_EC,S_Lc$ static S_FDot$ extern AMSTeX$ -Muldif.rr="00190620"$ +extern Glib_math_coordinate$ +extern Glib_canvas_x$ +extern Glib_canvas_y$ +Muldif.rr="00200515"$ AMSTeX=1$ TeXEq=5$ TeXLim=80$ @@ -783,6 +802,61 @@ def myediff(P,X) return red(X*diff(P,X)); } +def mypdiff(P,L) +{ + if(type(P)>3) return map(os_md.mypdiff,P,L); + for(Q=0;L!=[];L=cdr(L)){ + Q+=mydiff(P,car(L))*L[1]; + L=cdr(L); + } + return red(Q); +} + +def pTaylor(S,X,N) +{ + if(!isvar(T=getopt(time))) T=t; + if(type(S)<4) S=[S]; + if(type(X)<4) X=[X]; + if(findin(T,varargs(S|all=2))>=0){ + S=cons(z_z,S);X=cons(z_z,X);FT=1; + }else FT=0; + LS=length(S); + FR=(getopt(raw)==1)?1:0; + if(!FR) R=newvect(LS); + else R=R1=[]; + for(L=[],I=0,TS=S,TX=X;I1;N--){ + S=mypdiff(S,L); + K*=++M; + for(TS=S,I=0,R1=[];TS!=[];TS=cdr(TS),I++){ + if(!FR) R[I]+=car(TS)*t^M/K; + else R1=cons(car(TS)/K,R1); + } + if(FR) R=cons(reverse(R1),R); + } + if(FT){ + if(!FR){ + S=newvect(LS-1); + for(I=1;I= 4){ S=(type(P) == 6)?size(P)[0]:0; P = m2l(P); - for(I = 0, Deg = -3; P != []; P = cdr(P), I++){ - if( (DT = mydeg(car(P),X)) == -2) + for(I = 0, Deg = -100000; P != []; P = cdr(P), I++){ + if( (DT = mydeg(car(P),X)) == -2&&type(X)!=4) return -2; if(DT > Deg){ Deg = DT; @@ -823,8 +897,19 @@ def mydeg(P,X) return (Opt==1)?([Deg,(S==0)?II:[idiv(II,S),irem(II,S)]]):Deg; } P = red(P); - if(deg(dn(P),X) == 0) - return deg(nm(P),X); + if(type(X)==2){ + if(deg(dn(P),X) == 0) + return deg(nm(P),X); + }else{ + P=nm(red(P)); + for(D=-100000,I=deg(P,X[1]);I>=0;I--){ + if(TP=mycoef(P,I,X[1])){ + TD=mydeg(TP,X[0])-I; + if(DS) S=T; + } + } + return S; + } R=0; if(type(V)!=4) - for (I = length(V)-1; I >= 0; I--) R+= V[I]^2; + for (I = length(V)-1; I >= 0; I--) R+= real(V[I])^2+imag(V[I])^2; else{ if(type(V[0])>3){ V=ltov(V[0])-ltov(V[1]); return dnorm(V); } - for(;V!=[]; V=cdr(V)) R+=car(V)^2; + for(;V!=[]; V=cdr(V)) R+=real(car(V))^2+imag(car(V))^2; } - return dsqrt(R); + return ctrl("bigfloat")?pari(sqrt,R):dsqrt(R); } def dvprod(V1,V2) @@ -3373,7 +3468,45 @@ def rsort(L,T,K) return (iand(K,1))? reverse(R):R; } +def llget(L,LL,LC) +{ + if(type(LL)==4){ + LM=length(L); + for(R=[];LL!=[];LL=cdr(LL)){ + if(isint(TL=car(LL))) R=cons(TL,R); + else{ + IM=(length(TL)==1)?(LM-1):TL[1]; + for(I=car(TL);I<=IM;I++) R=cons(I,R); + } + } + LL=reverse(R); + if(LC==-1){ + LL=lsort(LL,[],1); + return lsort(L,"num",["sub"]|c1=LL); + } + L=lsort(L,"num",["get"]|c1=LL); + } + if(type(LC)==4){ + LM=length(L[0]); + for(R=[];LC!=[];LC=cdr(LC)){ + if(isint(TL=car(LC))) R=cons(TL,R); + else{ + IM=(length(TL)==1)?(LM-1):TL[1]; + for(I>=car(TL);I<=IM;I++) R=cons(I,R); + } + } + LC=reverse(R); + if(LL==-1){ + LC=lsort(LC,[],1); + return lsort(L,"col",["setminus"]|c1=LC); + } + L=lsort(L,"col",["put"]|c1=LC); + } + if(getopt(flat)==1) L=m2l(L|flat=1); + return L; +} + def lsort(L1,L2,T) { C1=getopt(c1);C2=getopt(c2); @@ -3418,9 +3551,8 @@ def lsort(L1,L2,T) }else{ for(I=0;LT!=[];I++,LT=cdr(LT)) if(findin(I,C1)<0) RT=cons(car(LT),RT); - RT=reverse(RT); } - R=cons(RT,R); + R=cons(reverse(RT),R); } return reverse(R); } @@ -3776,24 +3908,49 @@ def lgcd(L) return []; } -def llcm(L) +def llcm(R) { - if(type(L)==4){ - F=getopt(poly); - V=car(L); - while((L=cdr(L))!=[]){ - if(V!=0){ - if((V0=car(L))!=0) - V=(F==1)?red(V*V0/gcd(V,V0)):ilcm(V,V0); + if(type(R)==5||type(R)==6) R=m2l(R); + if(type(R)<4) R=[R]; + if(type(R)!=4) return 0; + V=getopt(poly); + if(type(V)<1){ + for(L=R;L!=[];L=cdr(L)){ + if(type(car(L))>1){ + V=1; break; } - else V=car(L); } - if(F!=1&&V<0) V=-V; - return V; } - else if(type(L)==5||type(L)==6) - return llcm(m2l(L)|option_list=getopt()); - return []; + if(getopt(dn)!=1){ + for(L=[];R!=[];R=cdr(R)) if(R!=0) L=cons(1/car(R),L); + R=L; + } + P=1; + if(type(V)<1){ + for(;R!=[];R=cdr(R)){ + if(!(TL=car(R))) continue; + else P=ilcm(P,dn(TL)); + } + return P; + } + for(;R!=[];R=cdr(R)){ + if(!car(R)) continue; + D=dn(red(car(R))); + N=red(P/D); + if(type(V)<2){ + if(type(N)!=3) continue; + P*=dn(N); + continue; + } + if(ptype(N,V)>2){ + L=fctr(dn(N)); + for(;L!=[];L=cdr(L)){ + if(ptype(car(L)[0],V)<2) continue; + P*=car(L)[0]^car(L)[1]; + } + } + } + return P; } def ldev(L,S) @@ -3870,6 +4027,9 @@ def lnsol(VV,L) def ladd(X,Y,M) { + if(Y==0){ + Y=X[1];X=X[0]; + } if(type(Y)==4) Y=ltov(Y); if(type(X)==4) X=ltov(X); return vtol(X+M*Y); @@ -4276,10 +4436,12 @@ def fctrtos(P) if(D>1) Op+="^"+((D>9)?(Pre+rtostr(D)+Post):rtostr(D)); PW=Op+Add+"}{"+PW+"}"; }else if(Add!=0) PW=PW+Add; + CD=0; if(TeX>=1){ if(type(PC)==1 && ntype(PC)==0 && PC<0) OC="-"+my_tex_form(-PC); else OC=fctrtos(PC|TeX=1,br=1); + if(isint(PC)&&(PC<-1||PC>1)) CD=1; }else OC=fctrtos(PC|br=1); if(PW!=""){ if(OC == "1") OC = ""; @@ -4310,7 +4472,13 @@ def fctrtos(P) }else L=LL; }else L+=LL; }else if(length(Tm)!=1) PW += CR; /* not final term */ - if(TeX) OC=texsp(OC); + if(TeX){ + OC=texsp(OC); + if(CD){ /* 2*3^x */ + CD=strtoascii(str_cut(PW,0,1)); + if(length(CD)==2&&car(CD)==123&&isnum(CD[1])) OC+="\\cdots"; + } + } if(str_chr(OC,0,"-") == 0 || C==0) str_tb([OC,PW], Out); else{ str_tb(["+",OC,PW],Out); @@ -4340,14 +4508,14 @@ def fctrtos(P) if(imag(P)==0) P = fctr(P); /* usual polynomial */ else P=[[P,1]]; S = str_tb(0,0); - for(J = N = 0; J < length(P); J++){ - if(type(P[J][0]) <= 1){ - if(P[J][0] == -1){ + for(J = N = CD = 0; J < length(P); J++){ + if(type(V=P[J][0]) <= 1){ + if(V == -1){ write_to_tb("-",S); if(length(P) == 1) str_tb("1", S); - }else if(P[J][0] != 1){ - str_tb((TeX>=1)?my_tex_form(P[J][0]):rtostr(P[J][0]), S); + }else if(V != 1){ + str_tb((TeX>=1)?my_tex_form(V):rtostr(V), S); N++; }else if(length(P) == 1) str_tb("1", S); @@ -4355,6 +4523,7 @@ def fctrtos(P) str_tb((TeX>=1)?my_tex_form(P[1][0]):rtostr(P[1][0]), S); J++; } + if(J==0&&isint(V=P[J][0])&&(V<-1||V>1)) CD=1; continue; } if(N > 0 && TeX != 1 && TeX != 2 && TeX != 3) @@ -4368,6 +4537,10 @@ def fctrtos(P) str_tb(["^", (TeX>1)?rtotex(P[J][1]):monotos(P[J][1])],S); }else{ if(nmono(P[J][0])>1) SS="("+SS+")"; + else if(CD&&J==1){ /* 2*3^x */ + CD=strtoascii(str_cut(SS,0,1)); + if(length(CD)==2&&car(CD)==123&&isnum(CD[1])) SS="\\cdot"+SS; + } write_to_tb(SS,S); } } @@ -4579,6 +4752,35 @@ def mtransbys(FN,F,LL) return call(FN, cons(F,LL)|option_list=Opt); } +def trcolor(S) +{ + if(type(S)!=7) return S; + return ((I=findin(S,LCOPT))>=0)?COLOPT[I]:0; +} + +def mcolor(L,P) +{ + if(type(L)!=4) return L; + if(!P||(S=length(L))==1){ + if(type(V=car(L))!=7) return V; + return trcolor(V); + } + P-=ceil(P)-1; + if(P==1){ + if(type(V=L[S-1])!=7) return V; + return trcolor(V); + } + for(S=P*(S-1);S>1;S--,L=cdr(L)); + if(getopt(disc)==1) S=0; + if(type(L0=L[0])==7) L0=trcolor(L0); + if(type(L1=L[1])==7) L1=trcolor(L1); + T=rint(iand(L0,0xff)*(1-S)+iand(L1,0xff)*S); + TT=iand(L0,0xff00)*(1-S)+iand(L1,0xff00)*S; + T+=rint(TT/0x100)*0x100; + TT=iand(L0,0xff0000)*(1-S)+iand(L1,0xff0000)*S; + return T+rint(TT/0x10000)*0x10000; +} + def drawopt(S,T) { if(type(S)!=7) return -1; @@ -4610,6 +4812,24 @@ def drawopt(S,T) return -1; } +def openGlib(W) +{ + extern Glib_canvas_x; + extern Glib_canvas_y; + extern Glib_math_coordinate; + + if(W==0){ + glib_clear(); + return; + } + if(type(W)==4&&length(W)==2){ + Glib_canvas_x=W[0]; + Glib_canvas_y=W[1]; + } + Glib_math_coordinate=1; + if(getopt(null)!=1) return glib_open(); +} + def execdraw(L,P) { if((Proc=getopt(proc))!=1) Proc=0; @@ -4862,6 +5082,12 @@ def execdraw(L,P) LOut=cons(T[2],Out); } } + }else if(T[0]==6){ /* plot */ + F++; + if((T1=findin(T[1],LCOPT))>-1) T1=COLOPT(T1); + else if(type(T1)!=1 && T1!=0) T1=0xffffff; + for(T2=ptaffine(M,T[2]|option_list=Org);T2!=[];T2=cdr(T2)) + draw_obj(Id,Ind,[rint(car(T2)[0]),rint(car(T2)[1])],T1); }else if(Proc==1&&type(T[0])==2){ if(length(T)<3) call(T[0],T[1]); else call(T[0],T[1]|option_list=T[2]); @@ -4941,9 +5167,15 @@ def execdraw(L,P) if(P[0]==2) dviout(T[2]|option_list=T[1]); else LOut=cons(T[2],Out); } + }else if(T[0]==6){ /* plot */ + F++; + if(type(T[1])==7) T1=[T[1],"."]; + else T1="."; + for(T2=ptaffine(M,T[2]|option_list=Org);T2!=[];T2=cdr(T2)) + str_tb(xypos([car(T2)[0],car(T2)[1],T1]),Out); }else if(T[0]==-2) str_tb(["%",T[1],"\n"],Out); - else if(Proc==1&&type(T[0])==2){ + else if(Proc==1&&type(T[0])==2){ if(length(T)<3) call(T[0],T[1]); else call(T[0],T[1]|option_list=T[2]); } @@ -5110,6 +5342,15 @@ def mmulbys(FN,P,F,L) def appldo(P,F,L) { + if(getopt(Pfaff)==1){ + L = vweyl(L); + X = L[0]; DX = L[1]; + for(I=mydeg(P,DX);I>0;I--){ + if(!(TP=mycoef(P,D,DX))) continue; + P=red(P+TP*(muldo(D^(I-1),F,L)-D^I)); + } + return P; + } if(type(F) <= 3){ if(type(L) == 4 && type(L[0]) == 4) return applpdo(P,F,L); @@ -5274,13 +5515,13 @@ def mce(P,L,V,R) { L = vweyl(L); X = L[0]; DX = L[1]; - P = sftexp(laplace1(P,L),L,V,R); + P = sftexp(laplace1(P,L),L,V,R|option_list=getopt()); return laplace(P,L); } def mc(P,L,R) { - return mce(P,L,0,R); + return mce(P,L,0,R|option_list=getopt()); } def rede(P,L) @@ -6521,27 +6762,20 @@ def toeul(F,L,V) L = vweyl(L); X = L[0]; DX = L[1]; I = mydeg(F,DX); - if(V == "infty"){ + if(getopt(raw)!=1){ for(II=I; II>=0; II--){ - J = mydeg(P=mycoef(F,I,DX),X); + J = mydeg(P=mycoef(F,II,DX),X); if(II==I) S=II-J; else if(P!=0 && II-J>S) S=II-J; } F *= X^S; - R = 0; - for( ; I >= 0; I--) + } + if(V == "infty"){ + for(R=0; I >= 0; I--) R += red((mysubst(mycoef(F,I,DX),[X,1/X])*(x*DX)^I)); return(subst(pol2sft(R,DX),DX,-DX)); } - F = subst(F,X,X+V); - for(II=I; II>=0; II--){ - J = mymindeg(P=mycoef(F,II,DX),X); - if(II==I) S=II-J; - else if(P!=0 && II-J>S) S=II-J; - } - F *= X^S; - R = 0; - for( ; I >= 0; I--) + for(R=0; I >= 0; I--) R += (red(mycoef(F,I,DX)/X^I))*DX^I; return pol2sft(R,DX); } @@ -6577,8 +6811,10 @@ def fromeul(P,L,V) S = DX*(S*X + mydiff(S,DX)); R += mycoef(P,J,DX)*S; } - while(mycoef(R,0,X) == 0) - R = tdiv(R,X); + if(getopt(raw)!=1){ + while(mycoef(R,0,X) == 0) + R = tdiv(R,X); + } if(V != "infty" && V != 0) R = mysubst(R,[X,X-V]); return R; @@ -6587,8 +6823,8 @@ def fromeul(P,L,V) def sftexp(P,L,V,N) { L = vweyl(L); DX = L[1]; - P = mysubst(toeul(P,L,V),[DX,DX+N]); - return fromeul(P,L,V); + P = mysubst(toeul(P,L,V|opt_list=getpt()),[DX,DX+N]); + return fromeul(P,L,V|opt_list=getopt()); } @@ -7175,6 +7411,28 @@ def pcoef(P,L,Q) return Coef; } +def pmaj(P) +{ + if(type(P)==4){ + Opt=getopt(var); + Opt=(isvar(Opt))?[["var",Opt]]:[]; + for(Q=[];P!=[];P=cdr(P)) Q=cons(pmaj(car(P)|option_list=Opt),Q); + if(Opt==[]) return reverse(Q); + X=Opt[0][1]; + D=mydeg(Q,X); + for(S=0;D>=0;D--) S+=lmax(mycoef(Q,D,X))*X^D; + return S; + } + V=vars(P); + if(!(K=length(V))) return abs(P); + for(R=0,D=deg(P,X=V[0]);D>=0;D--){ + Q=coef(P,D,X); + if(Q!=0) R+=((type(Q)>1)?pmaj(Q):abs(Q))*X^D; + } + if(isvar(Y=getopt(var))) for(;V!=[];V=cdr(V)) R=subst(R,car(V),Y); + return R; +} + def prehombf(P,Q) { if((Mem=getopt(mem))!=1 && Mem!=-1) @@ -8312,156 +8570,15 @@ def spgen(MO) }else{ L0=0; L1=MO+1; } - if(MO<=0){ + if(M0<=0){ MO=-MO; if(iand(MO,1)==1) return []; - if(MO>1){ - if(isMs()==0) return []; - Cmd="okubo "+rtostr(-MO); - MO/=2; - if(L1>0) Cmd=Cmd+"+"+rtostr(L0)+"-"+rtostr(L1); - else L1=MO+4; - Cmd=Cmd+" B"; - Id=getbyshell(Cmd); - if(Id<0) return []; - B=[]; - while((S=get_line(Id)) !=0){ - P0=str_chr(S,1,":")+1; - if(P0>1){ - P1=str_chr(S,P,"\n"); - if(P1<0) P1=str_len(S); - B=cons(sub_str(S,P0,P1-1),B); - } - } - close_file(Id); - }else{ - MO/=2; - if(L1<=1) L1=MO+4; -BB=[ -["11,11,11,11","111,111,111","1^4,1^4,22","1^6,222,33"], -["11,11,11,11,11","1^4,1^4,211","211,22,22,22","1^6,2211,33", -"2211,222,222","22211,2^4,44","2^511,444,66","1^4,22,22,31", -"2^5,3331,55","1^5,1^5,32","1^8,332,44","111,111,21,21","1^5,221,221"], -["11,11,11,11,11,11","1^4,1^4,1^4","1^4,22,22,22","111,111,111,21", -"1^6,21^4,33","21^4,222,222","221^4,2^4,44","2^41^4,444,66", -"1^5,1^5,311","1^8,3311,44","1^6,222,321","321,33,33,33", -"3321,333,333","33321,3^4,66","3^721,666,99","2^5,3322,55", -"1^6,1^6,42","222,33,33,42","1^a,442,55","1^6,33,33,51", -"222,222,33,51","1^9,333,54","2^7,554,77","1^5,2111,221", -"2^41,333,441","1^7,2221,43","211,211,22,22","2211,2211,222", -"22211,22211,44","1^4,211,22,31","2^411,3331,55","1^4,1^4,31,31", -"22,22,22,31,31","1^7,331,331","2221,2221,331","111,21,21,21,21"], -["11,11,11,11,11,11,11","111,111,111,111","1^6,1^6,33", -"1^6,222,222","222,33,33,33","1^5,1^5,221", -"1^4,211,22,22","1^4,1^4,22,31","22,22,22,22,31", -"111,111,21,21,21","21^6,2^4,44","2221^6,444,66", -"1^6,222,3111","3111,33,33,33","33111,333,333", -"333111,3^4,66","3^5111,666,99","2^5,33211,55", -"1^8,3221,44","3222,333,333","33222,3^4,66", -"3^4222,666,99","1^6,1^6,411","222,33,33,411", -"1^a,4411,55","2^4,2^4,431","431,44,44,44", -"2^6,4431,66","4431,444,444","44431,4^4,88", -"4^531,888,cc","1^a,433,55","1^7,1^7,52", -"1^c,552,66","3^4,444,552","1^8,2^4,53", -"1^8,44,44,71","3^5,555,771","21^4,2211,222", -"221^4,22211,44","2221^4,3331,55","1^6,2211,321", -"2^411,3322,55","1^7,322,331","2211,33,33,42", -"3^42,4442,77","2211,222,33,51","3^51,5551,88", -"2^611,554,77","2221,2221,322","2^41,2^41,54", -"1^5,2111,2111","222111,333,441","1^7,22111,43", -"1^5,1^5,41,41","1^9,441,441","22111,2221,331", -"1^5,221,32,41","221,221,221,41","211,211,211,22", -"2211,2211,2211","1^4,211,211,31","211,22,22,31,31", -"1^4,22,31,31,31","1^5,32,32,32","221,221,32,32","21,21,21,21,21,21"], -["11,11,11,11,11,11,11,11","1^4,1^4,22,22","1^8,2^4,44", -"1^6,2211,222","2211,33,33,33","111,111,111,21,21", -"1^5,1^5,2111","1^4,211,211,22","1^4,1^4,211,31", -"211,22,22,22,31","1^4,22,22,31,31","111,21,21,21,21,21", -"221^8,444,66","2^5,331^4,55","1^8,32111,44", -"32211,333,333","332211,3^4,66","3^42211,666,99", -"2^5,32221,55","1^7,1^7,511","1^c,5511,66", -"3^4,444,5511","541,55,55,55","5541,555,555", -"55541,5^4,aa","5^541,aaa,ff","1^8,1^8,62", -"1^a1^4,662,77","1^a,55,55,91","2^71,555,87", -"21^6,22211,44","221^6,3331,55","1^6,2211,3111", -"2^411,33211,55","1^7,3211,331","2211,33,33,411", -"3^42,44411,77","22211,2^4,431","2^511,4431,66", -"1^8,332,431","3^42,4433,77","1^8,22211,53", -"2221,2221,3211","221^5,333,441","1^7,21^5,43", -"1^b,443,65","21^5,2221,331","2^51,3332,65", -"21^4,21^4,222","221^4,221^4,44","1^6,21^4,321", -"2221^4,3322,55","21^4,33,33,42","21^4,222,33,51", -"2^51^4,554,77","2^4,3311,3311","3^411,4442,77", -"321,321,33,33","3321,3321,333","33321,33321,66", -"222,321,33,42","1^6,321,33,51","222,222,321,51", -"1^9,3321,54","1^7,322,322","3^422,5551,88", -"1^6,33,42,42","1^6,222,42,51","33,33,33,42,51", -"1^6,1^6,51,51","222,33,33,51,51","1^b,551,551", -"1^5,221,311,41","2^41,3321,441","22111,2221,322", -"2^51,443,551","222111,2^41,54","21^4,2211,2211", -"1^5,311,32,32","3331,3331,442","2211,2211,33,51", -"221,221,311,32","22111,22111,331","1^5,2111,32,41", -"2111,221,221,41","2111,221,32,32","211,211,211,211", -"211,211,22,31,31","1^4,211,31,31,31","22,22,31,31,31,31"], -["11,11,11,11,11,11,11,11,11","1^5,1^5,1^5","2^5,2^5,55", -"111,111,111,111,21","2^41,333,333","1^4,1^4,211,22", -"211,22,22,22,22","1^8,22211,44","1^4,1^4,1^4,31", -"1^4,22,22,22,31","1^7,1^7,43","1^7,2221,331", -"2221,2221,2221","1^6,21^4,222","21^4,33,33,33", -"1^6,1^6,321","222,321,33,33","1^6,33,33,42", -"222,222,33,42","1^6,222,33,51","222,222,222,51", -"33,33,33,33,51","1^6,2211,2211","111,111,21,21,21,21", -"1^5,1^5,32,41","1^5,221,221,41","1^5,221,32,32", -"221,221,221,32","1^4,211,211,211","211,211,22,22,31", -"1^4,211,22,31,31","1^4,1^4,31,31,31","22,22,22,31,31,31", -"21,21,21,21,21,21,21","21^a,444,66","1^8,31^5,44", -"321^4,333,333","3321^4,3^4,66","3^421^4,666,99", -"2^5,322111,55","32^41,3^4,66","3332^41,666,99", -"1^8,1^8,611","2^4,44,44,611","1^d,6611,77", -"4^5,66611,aa","2^6,444,651","3^4,3^4,651", -"651,66,66,66","3^6,6651,99","6651,666,666", -"66651,6^4,cc","6^551,ccc,ii","2^8,655,88", -"1^9,1^9,72","1^g,772,88","1^c,444,75", -"2^6,3^4,75","1^c,66,66,b1","3^4,444,66,b1", -"3^7,777,ba","1^7,2221,4111","2^41,333,4311", -"1^9,2^41,63","21^8,3331,55","2^411,331^4,55", -"1^7,31^4,331","2^411,32221,55","22211,2^4,422", -"2^511,4422,66","1^8,332,422","2^5,3331,541", -"22211,44,44,62","2^411,2^5,64","2^711,664,88", -"1^a,3331,64","2221,2221,31^4","21^7,333,441", -"333,333,441,81","2^6111,555,87","21^6,221^4,44", -"221^6,3322,55","2^41^6,554,77","1^6,21^4,3111", -"3111,321,33,33","33111,3321,333","333111,33321,66", -"222,3111,33,42","1^6,3111,33,51","222,222,3111,51", -"1^9,33111,54","2221^4,33211,55","1^7,3211,322", -"3^4211,5551,88","2^4,3221,3311","333221,4442,77", -"3222,3321,333","33222,33321,66","1^9,3222,54", -"21^4,33,33,411","3^411,44411,77","222,321,33,411", -"1^6,33,411,42","1^6,222,411,51","33,33,33,411,51", -"221^4,2^4,431","2^41^4,4431,66","1^8,3311,431", -"3^411,4433,77","33321,444,552","1^8,221^4,53", -"3311,44,44,53","4^42,5553,99","2^4,3311,44,71", -"3^421,555,771","4^52,7771,bb","3^611,776,aa", -"2^41,33111,441","22111,2221,3211","2^41,3222,441", -"2^61,4441,76","3331,3331,4411","22211,22211,431", -"3331,3331,433","3^41,3^41,76","1^7,1^7,61,61", -"1^d,661,661","21^5,2221,322","221^5,2^41,54", -"2^51,33311,65","21^5,22111,331","3^41,4441,661", -"1^7,331,43,61","2221,2221,43,61","2221,331,331,61", -"21^4,21^4,2211","21^4,2211,33,51","22211,3311,3311", -"1^5,311,311,32","2211,321,33,42","2211,222,321,51", -"3322,3331,442","2211,222,42,42","2^411,442,442", -"1^6,2211,42,51","2211,33,33,51,51","221,221,311,311", -"1^5,2111,311,41","222111,3321,441","22111,22111,322", -"222111,222111,54","2111,221,311,32","2111,2111,221,41", -"1^5,221,41,41,41","2221,43,43,43","1^5,32,32,41,41", -"331,331,43,43","221,221,32,41,41","221,32,32,32,41", -"211,211,211,31,31","211,22,31,31,31,31","1^4,31,31,31,31,31"]]; - B=BB[MO]; - } + MO=MO/2; + B=spbasic(-2*MO,0|str=1); + if(L1<3) L1=MO+4; if(St!=1){ for(R=[]; B!=[]; B=cdr(B)){ - RT=F?s2sp(car(B)|std=F):s2sp(car(B)); + RT= F?s2sp(car(B)|std=F): s2sp(car(B)); if(length(RT)L1) continue; R=cons(RT,R); } @@ -8563,6 +8680,198 @@ BB=[ return LL; } +def spbasic(Idx,D) +{ +/* + D<=3|Idx|+6, D<=|Idx|+2 (p>3), p<=|Idx|/2+4 + Idx=2*D^2-(D^2-\sum m_{j,\nu}^2); \sum(D-m_{j,1})>=2*D; + \sum (m_{j,1)-m_{j,\nu})*m_{j,\nu) + 0<=(2*D-\sum(D-m_{j,1})})*D=\sum_(m_{j,1}-m_{j,\mu})*m_{j,\nu} -|Idx| + (-2,0) 13ŒÂ (9+3+?) + (-4,0) 37ŒÂ (25+9+?) + (-6,0) : 8.5sec ?sec 0.05sec 69ŒÂ (46+17+?) + (-8,0) : 97 sec 1sec 0.13sec 113ŒÂ (73+29+?) <- (-2,0) + (-10,0): 4sec 0.27sec 198ŒÂ (127+50+?) +@(-12,0) 28sec 4.2sec 0.64sec 291ŒÂ (182+76+?) + (-14,0) 27sec 10.2sec 1.31sec 415ŒÂ (249+115+?) + (-16,0) 34.0sec 2.47sec 647ŒÂ (395+172+?) <- (-4,0) + (-18,0) 4.42sec 883ŒÂ (521+243+?) <- (-2,0) + (-20,0) 8.17sec 1186ŒÂ (680+345+?) +*/ + Idx=-Idx; + if((Str=getopt(str))!=1) Str=0; + if(!isint(Idx)||!isint(Idx/2)||Idx<0||!isint(D)||D<0||D==1||D>3*Idx+6) return []; + if(D==0){ + for(R=[],D=3*Idx+6;D>=2;D--) R=append(spbasic(-Idx,D|str=Str),R); + return R; + } + if(!Idx){ + R=0; + if(D==2) R="11,11,11,11"; + if(D==3) R="111,111,111"; + if(D==4) R="22,1111,1111"; + if(D==6) R="33,222,111111"; + if(!R) return []; + return [(Str==1)?R:s2sp(R)]; + } + if(D>Idx+2){ + L=3; + if(D==3*Idx+6){ + R=[[D/2,D/2],[D/3,D/3,D/3],[D/6,D/6,D/6,D/6,D/6,D/6-1,1]]; + return [(Str==1)?s2sp(R):R]; + } + if(iand(D,1)&&(D-3)/2>Idx) return []; + }else L=Idx/2+4; + V=newvect(L);SV=newvect(L); + for(S1=[],I=0;I1;T--){ + K=D%T; + if((T-K)*K<=Idx) break; + } + J=(T-K)*K;SJ=K^2+(D-K)*T; + TV=K?[K]:[]; + for(I=(D-K)/T;I>0;I--) TV=cons(T,TV); + for(I=0;I0) return []; + if(D>Idx+2 && V[0][0]+V[1][0]>=D && V[1][0]>1){ + T=V[1][0]-1;K=D%T;TV=K?[K]:[]; + for(I=(D-K)/T;I>0;I--) TV=cons(T,TV); + V[1]=V[2]=TV; + } + for(R=[];;){ + if(D>Idx+2){ + if(3*V[0][0]=D && (T=D-V[0][0]-1)>0){ + K=D%T;TV=K?[K]:[]; + for(I=(D-K)/T;I>0;I--) TV=cons(T,TV); + V[1]=V[2]=TV; + } + S2=V[0][0]+V[1][0]+V[2][0]-D; + if(V[0][0]+2*V[1][0]0){ + J=D%T; + K=J?[J]:[]; + for(J=(D-J)/T;J>0;J--) K=cons(T,K); + V[2]=K; + } + continue; + } + if(S2<0||V[2][0]<=S2){ + V[1]=V[2]=nextpart(V[1]); + continue; + }else if(S2>0){ + T=V[2][0]-S2;J=D%T; + K=J?[J]:[]; + for(J=(D-J)/T;J>0;J--) K=cons(T,K); + V[2]=K; + } + } + for(S=-2*D,IL=0;IL=0) break; + } + if((I=IL)==L){ /* reducible i.e. IL=L && S<0 */ + for(LL=L-1;LL>=0;LL--){ + if((K=car(V[LL]))+S>0){ + K+=S; + for(TV=[],TD=D;TD>=K;TD-=K) TV=cons(K,TV); + if(TD>0) V[LL]=append(TV,[TD]); + else V[LL]=TV; + break; + }else{ + S+=K-1; + V[LL]=S1; + } + } + if(LL<0) break; + continue; + } + for(S0=K=0;K<=IL;K++){ + ST=car(V[K]);J=V[K][length(V[K])-1];S0+=(ST-J)*J; + if(S0>Idx) break; + } + if(S0>Idx && car(V[K])!=1){ + ST=car(V[K]); + S0-=(ST-J)*J; + for(ST--;ST>0;ST--){ + J=D%ST; + if(S0+(ST-J)*J <= Idx) break; + } + V[K]=J?[J]:[]; + for(J=D-J;J>0;J-=ST) V[K]=cons(ST,V[K]); + for(J=K+1;JIdx && K<=IL && K!=L){ + SS0=Idx-SS+S0; + for(TV=car(V[K]);TV>1;TV--){ + U=D%TV; + if((D-U)*U<=SS0) break; + } + if(TV==car(V[K])){ + K=K-1; + V[K]=nextpart(V[K]); /* to be improves */ + }else{ + V[K]=U?[U]:[]; /* to be improved */ + for(J=D-U;J>0;J-=TV) V[K]=cons(TV,V[K]); + } + for(J=K+1;J=I||IL==2)){ + for(TR=[],K=J;K>=0;K--) TR=cons(V[K],TR); + R=cons((Str==1)?s2sp(TR):TR,R); + } + if(J>=0 && J1&& IxF-D^2+S0<0){ + for(V[J]=[],K=D-I;K>0;K--) V[J]=cons(1,V[J]); + V[J]=cons(I,V[J]); + V[J]=nextpart(V[J]); + for(I=J+1;I=0 && J(U=V[J][length(V[J])-1])+1){ + TV=reverse(V[J]); + for(S0=0,K=[];TV!=[];TV=cdr(TV),S0++){ + if((I=car(TV))1&&S0<2)){ + while(I-->0) K=cons(1,K); + }else K=cons(car(TV),K); + } + V[I=J]=K; + }else{ + if(J>=L) J=L-1; + for(I=J;I>=0&&length(V[I])==D;I--); + if(I<0) break; + } + V[I]=nextpart(V[I]); /* to be improved */ + for(J=I+1;J128){ - if(V<160 || (V>223 && V<240)) J++; + if(SJIS && (V=car(S))>128){ + if((V<160 || (V>223 && V<240))&&S!=[]) { + J++;S=cdr(S); + } } continue; } @@ -10816,6 +11127,7 @@ def my_tex_form(S) } SS = cons(S[I], SS); } + SS=str_subst(SS,"\n\\\\\n\\end{pmatrix}","\n\\end{pmatrix}"|raw=1); SS=str_subst(SS,"\\\\\n\\end{pmatrix}","\n\\end{pmatrix}"|raw=1); Subst=getopt(subst); Sub0=["{asin}","{acos}","{atan}"]; @@ -10890,7 +11202,7 @@ def my_tex_form(S) S=cons(123,S); if(F==2) SS=cdr(SS); else if(F==0) S=cons(car(SS),S); - }else if(F==2&&P-Q==3){ /* (2)^x -> 2^x*/ + }else if(F==2&&P-Q==3){ /* (2)^x -> 2^x */ SS=cdr(SS);SS=cdr(SS); S=cons(123,S);S=cons(car(SS),S);S=cons(125,S); SS=cdr(SS);SS=cdr(SS); @@ -10911,7 +11223,16 @@ def my_tex_form(S) SS=reverse(S); Top=P; } - S=asciitostr(SS); + for(F=G=0,S=[];SS!=[];SS=cdr(SS)){ /* 22^x -> 2\cdot 2^x */ + if(F==1&&G!=-1&&car(SS)==123 && length(SS)>1 && isnum(SS[1])) + S=append([116,111,100,99,92],S); + G=F; + if(car(SS)==125||car(SS)==95) F=-1; + else F=isnum(car(SS)); + S=cons(car(SS),S); + } + S=asciitostr(reverse(S)); +/* S=asciitostr(SS); */ if((K=getopt(ket))==1) S=texket(S); else if(K==2) S=texket(S|all=1); return S; @@ -11230,7 +11551,7 @@ def tocsv(L) if(type(LT)==5) LT=vtol(LT); if(type(LT)<4) LT=[LT]; for(N=0; LT!=[]; LT=cdr(LT),N++){ - if(N) str_tb(", ",Tb); + if(N) str_tb(",",Tb); if((T=car(LT))==Null) continue; if(type(T)==7){ K=str_len(T); @@ -11320,6 +11641,22 @@ def readcsv(F) return L; } +def getline(ID) +{ + if(isint(Maxlen=getopt(Max))>0) Maxlen=1024; + if(type(CR=getopt(CR))!=4) CR=[13]; + if(type(LF=getopt(LF))!=4) LF=[10]; + S=[]; + for(I=0; I<1023; I++){ + C=get_byte(ID); + if(C<0) return 0; + if(findin(C,CR)>=0) continue; + if(findin(C,LF)>=0) break; + S=cons(C,S); + } + return asciitostr(reverse(S)); +} + def showbyshell(S) { Id = getbyshell(S); @@ -11346,14 +11683,27 @@ def getbyshell(S) return open_file(F); } +def isfctr(P) +{ + if(type(P)>3) return 0; + if(type(P)==3) return (!isfctr(nm(P))||!isfctr(dn(P)))?0:1; + V=ptol(P,vars(P)|opt=0); + for(;V!=[];V=cdr(V)){ + if(type(car(V))>1||ntype(car(V))>0) return 0; + } + return 1; +} + def show(P) { T=type(P); S=P; Var=getopt(opt); + if((Raw=getopt(raw))!=1) Raw=0; if(Var=="verb"){ - dviout("{\\tt"+verb_tex_form(T)+"}\n\n"); - return; + S="{\\tt"+verb_tex_form(T)+"}\n\n"; + if(Raw) return S; + dviout(S);return; } if(type(Var)<0) Var=getopt(var); if(T==6){ @@ -11371,23 +11721,28 @@ def show(P) if(Var=="pfrac") X=var(P); else X=getopt(pfrac); if(isvar(X)){ - pfrac(P,X|dviout=1); - return; + if(Raw) return pfrac(P,X|TeX=1); + pfrac(P,X|dviout=1);return; } - Opt=cons(["dviout",1],getopt()); - if(type(Var)==2||type(Var)==4||type(Var)==7) fctrtos(P|option_list=Opt); - else{ + Opt=getopt(); + if(type(Var)!=2&&type(Var)!=4&&type(Var)!=7){ if(isdif(P)!=0) Opt=cons(["var","dif"],Opt); else Opt=cons(["br",1],Opt); - fctrtos(P|option_list=Opt); } - return; + if(!isfctr(P)){ + if(Raw) return my_tex_form(P); + else{ + dviout(P); return; + } + } + if(Raw) return fctrtos(P|option_list=cons(["TeX",3],Opt)); + fctrtos(P|option_list=cons(["dviout",1],Opt));return; }else if(T==4){ if(type(Var)==4 || type(Var)==7){ S=ltotex(P|option_list=getopt()); if(Var=="text"){ - dviout(S); - return; + if(Raw) return S; + dviout(S);return; } }else{ for(F=0,L=P;L!=[] && F!=-1;L=cdr(L)){ @@ -11415,8 +11770,8 @@ def show(P) if(F==1) S=ltotex(P|opt="spt"); else if(F==2){ M=mtranspose(lv2m(S)); - show(M|sp=1); /* GRS */ - return; + if(Raw) return show(M|sp=1,raw=1); /* GRS */ + show(M|sp=1);return; }else if(F==7) S=ltotex(P|opt="spts"); else{ for(F=0,L=P;L!=[] && F!=-1;L=cdr(L)){ @@ -11451,11 +11806,13 @@ def show(P) if(Var=="raw" || (Var !="eq" && str_chr(P,0,"\\")<0 && str_char(P,0,"^")<0 && str_char(P,0,"_")<0 && str_char(P,0,"&")<0)){ - dviout(P+"\n\n"); - return; + S=P+"\n\n"; + if(Raw) return S; + dviout(S); return; } } - dviout(S|eq=5); + if(Raw) return "\\begin{align}\\begin{split}\n &"+S+"\\end{split}\\end{align}"; + else dviout(S|eq=5); } @@ -12436,14 +12793,14 @@ def rungeKutta(F,N,Lx,Y,IY) if((Pr=getopt(prec))==1){ One=eval(exp(0)); }else{ - One=1;Pr=0; + One=deval(exp(0));Pr=0; } - if((FL=getopt(last))!=1) FL=0; + if(!isint(FL=getopt(mul))||!FL) FL=1; if(length(Lx)>2){ V=car(Lx);Lx=cdr(Lx); }else V=x; - if(Pr==0) Lx=[deval(Lx[0]),deval(Lx[1])]; - else Lx=[eval(Lx[0]),eval(Lx[1])]; + if(Pr==1) Lx=[eval(Lx[0]),eval(Lx[1])]; + else Lx=[deval(Lx[0]),deval(Lx[1])]; if(type(Y)==4){ if((Sing=getopt(single))==1||type(F)!=4) F=append(cdr(Y),[F]); @@ -12457,13 +12814,14 @@ def rungeKutta(F,N,Lx,Y,IY) } if(getopt(val)==1) V1=1; else V1=0; - H=(Lx[1]-Lx[0])/N;H2=H/2; + if(FL>0) N*=FL; + H=(Lx[1]-Lx[0])/N*One;H2=H/2; FV=findin(V,vars(F)); K=newvect(4); if(L==1){ R=[[T=Lx[0],S=IY]]; if(!H) return R; - for(;;){ + for(C=0;C0) break; + if(FL>0&&!((C+1)%FL)) R=cons([deval(T),S],R); } }else{ T=Lx[0]; R=[cons(T,V1?[car(IY)]:IY)]; S=ltov(IY); if(!H) return R; - for(;;){ + for(C=0;C0) break; + R=cons(cons(deval(T),TS),R); } } - return FL?(V1?S[0]:S):reverse(R); + L=(FL<0)?(V1?S[0]:S):reverse(R); + return L; } +def pwTaylor(F,N,Lx,Y,Ly,M) +{ + /* Pr:bigfloat, V1:last, Sf: single, Tf: autonomous, */ + if(!isint(FL=getopt(mul))||!FL) FL=1; + if(getopt(val)==1) V1=1; + else V1=0; + if(length(Lx)>2){ + V=car(Lx);Lx=cdr(Lx); + }else V=t; + if(!isvar(T=getopt(var))) V=t; + if(isint(Pr=getopt(prec))&&Pr>0){ + One=eval(exp(0)); + if(Pr>9){ + setprec(Pr); + ctrl("bigfloat",1); + } + Pr=1; + }else{ + One=deval(exp(0));Pr=0; + } + if(Pr==1) Lx=[eval(Lx[0]),eval(Lx[1])]; + else Lx=[deval(Lx[0]),deval(Lx[1])]; + Sf=(type(F)!=4)?1:0; + if(type(Y)==4){ + if(type(F)!=4) F=append(cdr(Y),[F]); + }else Y=[Y]; + if(type(Ly)!=4) Ly=[Ly]; + if(findin(V,vars(F))>=0){ + if(type(F)!=4) F=[F]; + Tf=1;F=cons(1,subst(F,V,z_z));Y=cons(z_z,Y);Ly=cons(car(Lx),Ly); + }else Tf=0; /* Tf: autonomous */ + ErF=0; + if(type(Er=getopt(err))==4){ + if(length(Er)==2) ErF=Er[1]; /* ErF&1: Raw, ErF&2: relative, ErF&4: add Sol */ + Er=car(Er); + }; + if(!isint(Er)||Er<0) Er=0; /* Šî€‰ð‚ð•Ô‚· */ + if(FL>0) N*=FL; + S=vtol(pTaylor(F,Y,M|time=V)); + FM=pmaj(F|var=x); + LS=length(S); + + if(type(Vw=getopt(view))==4){ /* Dislay on Canvas */ + Glib_math_coordinate=1; + glib_window(car(Vw)[0], car(Vw)[2],car(Vw)[1],car(Vw)[3]); + if(length(car(Vw))==6) Vr=[car(Vw)[4],car(Vw)[5]]; + else Vr=0; + if(length(Vw)>1){ + if(type(Cl=Vw[1])==4) Cl=map(os_md.trcolor,Cl); + else Cl=trcolor(Cl); + }else Cl=0; + if(length(Vw)>2){ + Mt=Vw[2]; + if(LS==1){ + if(type(Mt)>1) Mt=0; + }else{ + if(type(Mt)!=6||((Ms=size(Mt)[0])!=2&&Ms!=3)) Mt=0; + if(Ms!=3) Vr=0; + } + if(Tf&&type(Mt)==6) Mt=newbmat(2,2,[[1,0],[0,Mt]]); + }else Mt=0; + if(!Mt){ + if(LS>1+Tf){ + if(Vr){ + Mt=newmat(3,LS);Mt[2+Tf][2+Tf]=1; + } + else Mt=newmat(2,LS); + Mt[Tf][Tf]=Mt[Tf+1][Tf+1]=1; + }else Mt=1; + if(LS==1+Tf||Sf) glib_putpixel(Lx[0],Mt*Ly[Tf]|color=mcolor(Cl,0)); + else{ + YT=Mt*ltov(Ly); + glib_putpixel(YT[0],YT[1]|color=mcolor(Cl,0)); + } + } + }else Vw=0; + + T=Lx[0]; + RE=R=(Tf)?[Ly]:[cons(T,Ly)]; + H=(Lx[1]-Lx[0])/N*One; + + Ck=N+1;CB=10;Ckm=2;MM=2;C1=1; + if(Ck<5) Ck=100; + if(type(Inf=getopt(Inf))==4&&length(Inf)>1&&Inf[0]>4){ /* explosion */ + Ck=Inf[0];Ckm=Inf[1]; + if(length(Inf)>2) MM=Inf[2]; + if(!isint(MM)||MM<1) MM=2; + if(length(Inf)>3) C1=Inf[3]; + if(type(C1)!=1||C1<0) C1=1; + if(length(Inf)>4) CB=Inf[4]; + }else if(isint(Inf)&&Inf>0&&Inf<100){ + MM=Inf+1;Ck=100; + }else Inf=0; + Ckm*=Ck; + + SS=subst(S,V,H);N0=N; + if(Er>0){ + HE=H/(Er+1);SSE=subst(S,V,HE);LyE=Ly; + } + for(C=CC=CF=0;C=Ck){ /* check explosion */ + CC=0; + D0=dnorm(Ly|max=1); + if(Er&&CF){ + DE=dnorm(ladd(LyE,Ly,-1)|max=1); + if(CB*DE>D0) break; + } + for(Dy=F,TY=Y,TL=Ly;TY!=[];TY=cdr(TY),TL=cdr(TL)) + Dy=subst(Dy,car(TY),One*car(TL)); + D1=dnorm(Dy|max=1);D2=subst(FM,x,2*D0+C1);D3=D1+D2; + HH=2*(D0+C1)/Ckm; + if(HHHH) H/=2; + if(H*7/51) N*=MM; + MM=0; + } + CC=0; + } + + T+=H; + for(Dy=SS,TY=Y,TL=Ly;TY!=[];TY=cdr(TY),TL=cdr(TL)) + Dy=subst(Dy,car(TY),One*car(TL)); + Ly=Dy; + + if(Er>0){ /* estimate error */ + for(CE=0;CE<=Er;CE++){ + for(Dy=SSE,TY=Y,TL=LyE;TY!=[];TY=cdr(TY),TL=cdr(TL)) + Dy=subst(Dy,car(TY),One*car(TL)); + LyE=Dy; + } + } + if(FL<0||(C+1)%FL) continue; + if(Vw){ + if(LS==1+Tf||Sf) CR=CC/N0; + else{ + YT=Mt*ltov(Ly); + CR=(!Vr)?CC/N0:(YT[2]-Vr[0])/(Vr[1]-Vr[0]); + } + if(LS==1+Tf||Sf) glib_putpixel(deval(T),Mt*Ly[Tf]|color=mcolor(Cl,CR)); + else glib_putpixel(YT[0],YT[1]|color=mcolor(Cl,CR)); + continue; + } + TR=(V1)?[car(Ly)]:Ly; + if(!Tf) TR=cons((Inf)?eval(T):deval(T),TR); + R=cons(TR,R); + if(Er){ + TRE=(V1)?[car(LyE)]:LyE; + if(!Tf) TRE=cons((Inf)?eval(T):deval(T),TRE); + RE=cons(TRE,RE); + } + } + if(Vw) return 1; + L=(FL<0)?((V1)?car(Ly):Ly):reverse(R); + if(Er){ /* Estimate error */ + LE=(FL<0)?((V1)?car(LyE):LyE):reverse(RE); + if(FL>0){ + for(S=L,T=LE,D=[];S!=[];S=cdr(S),T=cdr(T)) D=cons(os_md.ladd(car(S),car(T),-1),D); + F=map(os_md.dnorm,reverse(D)); + if(iand(ErF,2)){ /* relative error */ + G=llget(LE,-1,[0]); + G=map(os_md.dnorm,G); + for(R=[];G!=[];G=cdr(G),F=cdr(F)){ + if(car(G)) R=cons(car(F)/car(G),R); + else R=cons(0,R); + } + F=reverse(R); + } + if(!iand(ErF,1)) F=map(os_md.nlog,F); + if(!iand(ErF,8)) F=map(deval,F); + }else if(V1){ + D=ladd(L,LE,-1);F=dnorm(D); + if(iand(ErF,2)){ + G=dnorm(cdr(L)); + if(!G) D/=G; + else D=1; + } + F=(!iand(ErF,1))?nlog(D):D; + if(!iand(ErF,8)) F=deval(F); + }else{ + D=abs(L-LE); + if(iand(ErF,2)){ + G=abs(L); + if(!G) D/=G; + else D=1; + } + F=(!iand(ErF,1))?nlog(D):D; + if(!iand(ErF,8)) F=deval(F); + } + return iand(ErF,4)?[L,F,LE]:[L,F]; + } + return L; +} + def xy2graph(F0,N,Lx,Ly,Lz,A,B) { /* (x,y,z) -> (z sin B + x cos A cos B + y sin A cos B, @@ -13084,10 +13644,15 @@ def mytan(Z) def mylog(Z) { if(type(Z=eval(Z))>1) return todf(os_md.mylog,[Z]); - if((Im=imag(Z))==0) return dlog(Z); + if(imag(Z)==0&&Z>=0) return dlog(Z); return dlog(dabs(Z))+@i*myarg(Z); } +def nlog(X) +{ + return mylog(X)/dlog(10); +} + def mypow(Z,R) { if(type(Z=eval(Z))>1||type(R=eval(R))>1) return todf(os_md.mypow,[Z,R]); @@ -13854,6 +14419,105 @@ def fcont(F,LX) return reverse(L); } +def xyplot(L,LX,LY) +{ + Vw=getopt(view); + if(type(Vw)!=1 && type(Vw)!=7 && Vw!=0) Vw=-1; + if(!LX){ + L0=llget(L,1,[0]|flat=1); + LX=[lmin(L0),lmax(L0)]; + S=LX[1]-LX[0]; + if(S>0){ + if(Vw) LX=[LX[0]-S/32,LX[1]+S/32]; + }else LX=[LX[0]-1,LX[0]+1]; + } + LX=map(deval,LX); + if(!LY){ + L0=llget(L,1,[1]|flat=1); + LY=[lmin(L0),lmax(L0)]; + S=LY[1]-LY[0]; + if(S>0){ + if(Vw) LY=[LY[0]-S/32,LY[1]+S/32]; + }else LY=[LY[0]-1,LY[0]+1]; + } + LY=map(deval,LY); + if(getopt(raw)==1) mycat([LX,LY]); + if(Vw!=-1){ + if(Vw!=1){ + if(type(Vw)==7) Vw=trcolor(Vw); + Opt=[["color",Vw]]; + }else Opt=[]; + Glib_math_coordinate=1; + glib_window(LX[0],LY[0],LX[1],LY[1]); + for(; L!=[];L=cdr(L)) + glib_putpixel(car(L)[0],car(L)[1]|option_list=Opt); + return [LX,LY]; + } + Opt=getopt();Opt0=delopt(Opt,["dviout","proc"]); + if(type(R=getopt(to))!=4) To=[12,8]; + R=[To[0]/(LX[1]-LX[0]),RY=To[1]/(LY[1]-LY[0])]; + R=[sint(R[0],4|str=0),sint(R[1],4|str=0)]; + S="% "; + if(type(C=getopt(scale))!=1&&type(C)!=4){ + Opt0=cons(["scale",R],Opt0); + S+="scale="+rtostr(R)+", "; + } + S+=rtostr(LX)+", "+rtostr(LY)+"\n"; + for(L0=[],TL=L;TL!=[];TL=cdr(TL)){ + TTL=map(deval,car(TL)); + if(TTL[0]LX[1]||TTL[1]LY[1]){ + S+=xylines(reverse(L0)|option_list=Opt0); + L0=[]; + }else{ + L0=cons(TTL,L0); + } + } + if(length(L0)>1) S+=xylines(reverse(L0)|option_list=Opt0); + AX=getopt(ax);Opt=delopt(Opt0,"opt"); + if(type(AX)==4) S+="% axis\n"+xygraph([0,0],0,LX,LX,LY|option_list=Opt); + else if((LX[0]<=0&&LX[1]>=0)||(LY[0]<=0&&LY[1]>=0)) + S+="% axis\n"+xygraph([0,0],0,LX,LX,LY|option_list=cons(["ax",[0,0]],Opt)); + if(getopt(dviout)!=1) return S; + xyproc(S|dviout=1); + return [LX,LY]; +} + +def xyaxis(A,X,Y) +{ + if(isint(Vw=getopt(view))&&Vw!=0){ + CL=getopt(opt); + if(type(CL)==7) CL=trcolor(CL); + if(type(CL)!=0) CL=0; + if(CL) Opt=[[color,CL]]; + else Opt=[]; + Glib_math_coordinate=1; + UX=(X[1]-X[0])/50;UY=(Y[1]-Y[0])/50; + glib_window(X[0],Y[0],X[1],Y[1]); + glib_line(A[0],Y[0],A[0],Y[1]|option_list=Opt); + glib_line(X[0],A[1],X[1],A[1]|otpion_list=Opt); + if(length(A)>2&&A[2]){ + I0=-floor((A[0]-X[0])/A[2]);I1=floor((X[1]-A[0])/A[2]); + for(I=I0;I<=I1;I++){ + IX=A[0]+A[2]*I; + if(iand(Vw,2)) glib_print(IX-UX,A[1]-UY/2,rtostr(IX)); + glib_line(IX,A[1],IX,A[1]+UY); + } + } + if(length(A)>3&&A[3]){ + I0=-floor((A[1]-Y[0])/A[3]);I1=floor((Y[1]-A[1])/A[3]); + for(I=I0;I<=I1;I++){ + IY=A[1]+A[3]*I; + if(iand(Vw,4)) glib_print(A[0]-UX*2,IY+UY,rtostr(IY)); + glib_line(A[0],IY,A[0]+UX,IY); + } + } + return; + } + Opt=getopt(); + Opt=cons(["ax",A],Opt); + return xygraph([0,0],0,[0,1],X,Y|option_list=Opt); +} + def xygraph(F,N,LT,LX,LY) { if((Proc=getopt(proc))!=1&&Proc!=2&&Proc!=3) Proc=0; @@ -14166,7 +14830,7 @@ def xygraph(F,N,LT,LX,LY) if(length(Ax)>3){ D=Ax[3]; if(type(D)==1 && D>0){ - I0=ceil((LY[0]-Ax[1])/D); I1=floor((LY[1]-Ax[0])/D); + I0=ceil((LY[0]-Ax[1])/D); I1=floor((LY[1]-Ax[1])/D); for(DD=[],I=I0; I<=I1; I++){ if(length(Ax)<5) DD=cons(I*D,DD); else if(I!=0){ @@ -15757,6 +16421,27 @@ def ptcopy(L,V) } } +def regress(L) +{ + E=deval(exp(0)); + for(S0=T0=0,S=L;S!=[];S=cdr(S)){ + S0+=car(S)[0]*E;T0+=car(S)[1]*E; + } + K=length(L);S0/=K;T0/=K; + for(SS=TT=0,S=L;S!=[];S=cdr(S)){ + SS+=(car(S)[0]-S0)^2*E;TT+=(car(S)[1]-T0)^2*E; + ST+=(car(S)[0]-S0)*(car(S)[1]-T0)*E; + } + if(!SS||!TT) return []; + A=ST/SS; + L=[A,A*S0-T0,ST/dsqrt(SS*TT),S0,dsqrt(SS/K),T0,dsqrt(TT/K)]; + if(isint(N=getopt(sint))){ + R=reverse(L); + for(L=[];R!=[];R=cdr(R)) L=cons(sint(car(R),N|str=0),L); + } + return L; +} + def average(L) { if(getopt(opt)=="co"){ @@ -17427,6 +18112,39 @@ def shiftop(M,S) return [QQ,P,RS]; } + +def shiftPfaff(A,B,G,X,M) +{ + if(type(G)==4){ + G0=G[1];G1=G[0]; + } + if(type(G)==6){ + G=map(red,G); + G0=llcm(G);G1=map(red,G0*G); + } + if(type(G)==3){ + G=red(G);G0=dn(G);G1=nm(G); + } + if(type(M)==4){ + M0=M[0];M1=M[1]; + }else{ + M0=M;M1=0; + } + X=vweyl(X); + D0=mydeg(G0,X[0]);D1=mydeg(G1,X[0]); + if(M1>=0){ + D=(D1-M1>D0)?D1-M1:D0; + G0=muldo(X[1]^D,G0,X);G1=muldo(X[1]^(D+M1),G1,X); + }else{ + D=(D0+M1>D1)?D0+M1:D1; + G0=muldo(X[1]^(D-M1),G0,X);G1==muldo(X[1]^D,G1,X); + } + G0=map(mc,G0,X,M0);G1=map(mc,G1,X,M0+M1); + G0=appldo(G0,A,X|Pfaff=1); + G1=sppldo(G1,B,X|Pfaff=1); + return rmul(myinv(G0),G1); +} + def conf1sp(M) { if(type(M)==7) M=s2sp(M); @@ -17753,6 +18471,7 @@ def mcvm(N) S+=N[I]; } M=newbmat(K,K,reverse(M)); + NR=N; N=S; }else{ if(type(X)==7) X=strtov(X); @@ -17767,22 +18486,73 @@ def mcvm(N) } } } - if(getopt(get)==1){ + if((Get=getopt(get))==1){ for(R=[],I=0;I0;I--) V=cons(makev(["a0",I]),V); + MI=myinv(M); + V=ltov(V)*MI; + for(R=[],I=0;I0;I--) V1=cons(os_md.makev([X[0],I]),V1); + for(V2=[],I=NR[1];I>0;I--) V2=cons(os_md.makev([X[1],I]),V2); + R=subst(R,car(V1),0,car(V2),0); + V1=subst(V1,car(V1),0); + V2=subst(V2,car(V2),0); + for(V=[],S=V1;S!=[];S=cdr(S)) for(T=V2;T!=[];T=cdr(T)) V=cons(car(T)-car(S),V); + V=reverse(V); + Mx=length(V); + for(A0=[],I=J=NR[0]-1;J>=0;I+=--J) for(K=0;K=0;I--) F0=cons(1/(x-V[I]), F0); + MV=confexp([F0,V]|sym=3); + RR=newvect(Mx); + for(K=0;K3) M=map(red,M); + else M=red(M); + R=cons(M,R); + } + R=reverse(R); + if(Sym==2) return R; + M=length(R);N=length(S[1]); + E=newmat(M,N); + for(I=0;I