=================================================================== RCS file: /home/cvs/OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v retrieving revision 1.47 retrieving revision 1.95 diff -u -p -r1.47 -r1.95 --- OpenXM/src/asir-contrib/packages/src/os_muldif.rr 2019/03/05 01:52:51 1.47 +++ OpenXM/src/asir-contrib/packages/src/os_muldif.rr 2022/07/22 23:09:09 1.95 @@ -1,12 +1,12 @@ -/* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.46 2019/03/04 02:02:04 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.94 2022/07/19 04:14:49 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 - July. 2022) * * For polynomials and differential operators with coefficients * in rational funtions (See os_muldif.pdf) @@ -21,6 +21,7 @@ module os_md; static Muldif.rr$ static TeXEq$ static TeXLim$ +static TeXPages$ static DIROUT$ static DIROUTD$ static DVIOUTL$ @@ -59,6 +60,10 @@ localf countin$ localf mycoef$ localf mydiff$ localf myediff$ +localf mypdiff$ +localf difflog$ +localf pTaylor$ +localf pwTaylor$ localf m2l$ localf m2ll$ localf mydeg$ @@ -77,15 +82,23 @@ localf ndict$ localf nextsub$ localf nextpart$ localf transpart$ +localf getCatalan$ +localf pg2tg$ +localf pg2tgb; +localf pgpart$ +localf xypg2tg; localf trpos$ -localf sprod$ +localf sprod$ +localf sadj$ localf sinv$ localf slen$ +localf sexps$ localf sord$ localf vprod$ localf dvangle$ localf dvprod$ localf dnorm$ +localf dext$ localf mulseries$ localf pluspower$ localf vtozv$ @@ -93,12 +106,16 @@ localf dupmat$ localf matrtop$ localf mytrace$ localf mydet$ +localf permanent$ localf mperm$ localf mtranspose$ localf mtoupper$ localf mydet2$ localf myrank$ -localf meigen$ +localf lext2$ +localf meigen$ +localf pf2kz$ +localf mext2$ localf transm$ localf vgen$ localf mmc$ @@ -113,9 +130,14 @@ localf myimage$ localf mymod$ localf mmod$ localf ladd$ +localf lsub$ localf lchange$ localf llsize$ localf llbase$ +localf llget$ +localf lcut$ +localf rev$ +localf qsortn$ localf lsort$ localf rsort$ localf lpair$ @@ -139,6 +161,7 @@ localf mpower$ localf mrot$ localf texlen$ localf isdif$ +localf isfctr$ localf fctrtos$ localf texlim$ localf fmult$ @@ -147,11 +170,15 @@ localf getel$ localf ptol$ localf rmul$ localf mtransbys$ +localf trcolor$ +localf mcolor$ localf drawopt$ localf execdraw$ localf execproc$ localf myswap$ localf mysubst$ +localf sort2$ +localf n2a$ localf evals$ localf myval$ localf myeval$ @@ -171,8 +198,12 @@ localf myasin$ localf myacos$ localf myatan$ localf mylog$ +localf nlog$ +localf dlog10$ localf mypow$ localf scale$ +localf catalan$ +localf iceil$ localf arg$ localf sqrt$ localf gamma$ @@ -200,13 +231,15 @@ localf mmulbys$ localf appldo$ localf appledo$ localf muldo$ +localf caldo$ localf jacobian$ localf hessian$ localf wronskian$ localf adj$ localf laplace1$ localf laplace$ -localf mce$ +localf mce$ +localf mcme$ localf mc$ localf rede$ localf ad$ @@ -216,7 +249,8 @@ localf addl$ localf cotr$ localf rcotr$ localf muledo$ -localf mulpdo$ +localf mulpdo$ +localf transppow$ localf transpdosub$ localf transpdo$ localf translpdo$ @@ -237,6 +271,10 @@ localf sftpowext$ localf polinsft$ localf pol2sft$ localf polroots$ +localf sgnstrum$ +localf polstrum$ +localf polrealroots$ +localf polradiusroot$ localf fctri$ localf binom$ localf expower$ @@ -245,7 +283,10 @@ localf seriesMc$ localf seriesTaylor$ localf mulpolyMod$ localf solveEq$ +localf res0$ +localf eqs2tex$ localf baseODE$ +localf baseODE0$ localf taylorODE$ localf evalred$ localf toeul$ @@ -260,6 +301,7 @@ localf expat$ localf polbyroot$ localf polbyvalue$ localf pcoef$ +localf pmaj$ localf prehombf$ localf prehombfold$ localf sub3e$ @@ -297,6 +339,7 @@ localf iscoef$ localf iscombox$ localf sproot$ localf spgen$ +localf spbasic$ localf chkspt$ localf cterm$ localf terms$ @@ -327,6 +370,8 @@ localf s2euc$ localf s2sjis$ localf r2ma$ localf evalma$ +localf evalcoord$ +localf readTikZ$ localf ssubgrs$ localf verb_tex_form$ localf tex_cuteq$ @@ -337,6 +382,7 @@ localf divmattex$ localf dviout0$ localf myhelp$ localf isMs$ +localf getline$ localf showbyshell$ localf readcsv$ localf tocsv$ @@ -344,6 +390,7 @@ localf getbyshell$ localf show$ localf dviout$ localf rtotex$ +localf togreek$ localf mtotex$ localf ltotex$ localf texbegin$ @@ -352,9 +399,13 @@ localf texsp$ localf getbygrs$ localf mcop$ localf shiftop$ +localf shiftPfaff; localf conf1sp$ localf confexp$ localf confspt$ +localf vConv$ +localf mcvm$ +localf s2cspb$ localf s2csp$ localf partspt$ localf pgen$ @@ -382,7 +433,8 @@ localf powsum$ localf bernoulli$ localf lft01$ localf linfrac01$ -localf nthmodp$ +localf nthmodp$ +localf issquare$ localf issquaremodp$ localf rootmodp$ localf rabin$ @@ -390,13 +442,16 @@ 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$ +localf xylabel$ localf xybox$ localf xyline$ localf xylines$ @@ -414,6 +469,8 @@ localf periodicf$ localf cmpf$ localf areabezier$ localf saveproc$ +localf xyplot$ +localf xyaxis$ localf xygraph$ localf xy2graph$ localf addIL$ @@ -424,13 +481,21 @@ localf xyarrows$ localf xyang$ localf xyoval$ localf xypoch$ +localf xycircuit$ +localf ptline$ localf ptcommon$ +localf ptinversion$ +localf ptcontain$ localf ptcopy$ localf ptaffine$ localf ptlattice$ localf ptpolygon$ localf ptwindow$ +localf pt5center$ +localf ptconvex$ localf ptbbox$ +localf darg$ +localf dwinding$ localf lninbox$ localf ptcombezier$ localf ptcombz$ @@ -446,6 +511,7 @@ localf msort$ extern Muldif.rr$ extern TeXEq$ extern TeXLim$ +extern TeXPages$ extern DIROUT$ extern DIROUTD$ extern DVIOUTL$ @@ -465,19 +531,24 @@ extern XYPrec$ extern XYcm$ extern TikZ$ extern XYLim$ +extern TeXPages$ 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="00190301"$ +extern Glib_math_coordinate$ +extern Glib_canvas_x$ +extern Glib_canvas_y$ +Muldif.rr="00220719"$ AMSTeX=1$ TeXEq=5$ TeXLim=80$ +TeXPages=20$ TikZ=0$ XYcm=0$ XYPrec=3$ @@ -644,7 +715,6 @@ def fcat(S,X) [getenv("HOME"),get_rootdir(),"/"])+"/"; if(isMs()) DIROUTD=str_subst(DIROUTD,"/","\\"|sjis=1); } - if(S==-1) return; T="fcat"; if(S>=2&&S<=9) T+=rtostr(S); T=DIROUTD+T+".txt"; @@ -782,6 +852,77 @@ 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 difflog(L) +{ + if(!isvar(X=getopt(var))) X=x; + if(type(L)!=4) return 0; + for(S=0;L!=[];L=cdr(L)){ + if(type(L0=car(L))==4) S+=L0[1]*mydiff(L0[0],X)/L0[0]; + if(type(L0)<4) S+=mydiff(L[0],X); + } + S=red(S); + if(type(F=getopt(mc))>0){ + X=vweyl(X); + S=mc(X[1]-S,X,F); + } + return red(S); +} + +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; @@ -822,8 +963,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(D 3) + if(type(P) > 3){ #ifdef USEMODULE return map(os_md.simplify,P,L,T); #else return map(simplify,P,L,T); #endif + } if(type(L[0]) == 4){ if(length(L[0]) > 1) #if USEMODULE @@ -1031,6 +1184,7 @@ def vnext(V) def ldict(N, M) { Opt = getopt(opt); + F=iand(Opt,4)/4;Opt=iand(Opt,3); R = S = []; for(I = 2; N > 0; I++){ R = cons(irem(N,I), R); @@ -1045,11 +1199,11 @@ def ldict(N, M) J++; } T[I-1] = 1; - S = cons(LL-I+1, S); + S = cons(LL-I+F+1, S); } for(I = 0; I <= LL; I++){ if(T[I] == 0){ - S = cons(LL-I, S); + S = cons(LL-I+F, S); break; } } @@ -1060,13 +1214,14 @@ def ldict(N, M) return 0; } T = []; - for(I = --M; I > LL; I--) - T = cons(I,T); + for(I = --M; I > LL;I--) + T = cons(I+F,T); S = append(S,T); if(Opt == 2 || Opt == 3) S = reverse(S); if(Opt != 1 && Opt != 3) return S; + M+=2*F; for(T = []; S != []; S = cdr(S)) T = cons(M-car(S),T); return T; @@ -1075,6 +1230,7 @@ def ldict(N, M) def ndict(L) { Opt = getopt(opt); + if(type(L)==5) L=vtol(L); R = []; if(Opt != 1 && Opt != 2) L = reverse(L); @@ -1142,6 +1298,10 @@ def transpart(L) def trpos(A,B,N) { + if(!N){ + N=(AR){ + W=newvect(L); + for(I=0;I= 0) - V[L] = S[T[L]]; - return V; + while(--L >= 0) V[L] = S[T[L]]; + return (F==2)?ndict(V):V; } + +def sadj(S,T) +{ + return sprod(sprod(S,T),sinv(S)); +} def sinv(S) { - L = length(S); + if(F=isint(S)) S=ltov(ldict(S,0)); + L = length(S); + if(L==2) return S; V = newvect(L); while(--L >= 0) V[S[L]] = L; - return V; + return (F)?ndict(V):V; } def slen(S) @@ -1176,6 +1367,18 @@ def slen(S) return V; } +def sexps(S) +{ + K=length(S);S=ltov(S); + for(R=[],I=0;I=0&&S[J]>S[J+1];J--){ + T=S[J];S[J]=S[J+1];S[J+1]=T; + R=cons(J,R); + } + } + return R; +} + def sord(W,V) { L = length(W); @@ -1209,6 +1412,7 @@ def sord(W,V) def vprod(V1,V2) { + V1=lsub(V1);V2=lsub(V2); for(R = 0, I = length(V1)-1; I >= 0; I--) R = radd(R, rmul(V1[I], V2[I])); return R; @@ -1216,24 +1420,36 @@ def vprod(V1,V2) def dnorm(V) { - if(type(V)<2) return dabs(V); + if(type(V)<2) return ctrl("bigfloat")?abs(V):dabs(V); + if((M=getopt(max))==1||M==2){ + if(type(V)==5) V=vtol(V); + for(S=0;V!=[];V=cdr(V)){ + if(M==2) S+=ctrl("bigfloat")?abs(car(V)):dabs(car(V)); + else{ + if((T=ctrl("bigfloat")?abs(car(V)):dabs(car(V)))>S) 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) { if(type(V1)<2) return V1*V2; R=0; + V1=lsub(V1); + V2=lsub(V2); if(type(V1)!=4) for(I = length(V1)-1; I >= 0; I--) R += V1[I]*V2[I]; @@ -1244,6 +1460,13 @@ def dvprod(V1,V2) return R; } +def ptline(L,R) +{ + P=L[0];Q=L[1]; + return (Q[1]-P[1])*(R[0]-P[0])-(Q[0]-P[0])*(R[1]-P[1]); +} + + def dvangle(V1,V2) { if(V2==0 && type(V1)==4 && length(V1)==3 && @@ -1275,6 +1498,741 @@ def mulseries(V1,V2) return VV; } +def catalan(K) +{ + if(isint(K)) return catalan([K,K]); + if(type(K)==4){ + if(length(K)==2){ + M=K[0];N=K[1]; + if(MN) return 0; + return fac(N)/fac(K)/fac(N-K); + } + if(K<1||N<1||K>N) return 0; + if(N==K) return 1; + if(T==1){ + if(K==1) return fac(N-1); + return catalan([1,N-1,K-1])+(N-1)*catalan([1,N-1,K]); + }else if(T==2){ + if(K==1) return 1; + return catalan([2,N-1,K-1])+ K*catalan([2,N-1,K]); + } + } + } + return 0; +} + +def sort2(L) +{ + if(L[0]<=L[1]) return L; + if(type(L)==4) return [L[1],L[0]]; + T=L[0];L[0]=L[1];L[1]=T; + return L; +} + +/* 01: 01 list + * s : 01 str + * T : tounament + * # : #lines of vertexes (vector) + * P : Polygon with tg + */ +def getCatalan(X,N) +{ + if(type(To=getopt(to))!=7) To=0; + if(type(X)==7){ /* string: s or T */ + X=strtoascii(X); + N=length(X); + if(X[0]==48){ + if(To=="s") return R; + R=calc(X,["-",48]); /* s -> 01 */ + if(To) R=getCatalan(R,0|to=To); + return R; + } + if(To=="T") return X; + if(To!="P"&&To!="#"){ /* T -> 01 */ + for(R=[];X!=[];X=cdr(X)){ + if(car(X)==41) R=cons(1,R); + else if(car(X)==42) R=cons(0,R); + } + R=cdr(reverse(R)); + if(To!="01") R=getCatalan(R,0|to=To); + return R; + } + if(N%3!=1) return 0; + M=(N+2)/3; /* T -> # */ + V=newvect(M+1); + V[0]=V[M]=-1; + for(;X!=[];X=cdr(X)){ + if(car(X)==40||car(X)==41) V[I]++; + else I++; + } + V[M]+=F; + if(To!="P") return V; + X=V; + } + if(type(X)==5){ /* vector: # -> P */ + if(To=="#") return X; + Y=newvect(length(X));K=dupmat(X); + M=length(X); + for(R=[],I=F=0;;I++){ + if(I>=M){ + if(!F) break; + F=0;I=-1;continue; + } + if(X[I]>0){ + if(I+1>=M ||K[I+1]>0) continue; + for(J=I+2;J=M||findin([I,J],R)>=0) continue; + R=cons([I,J],R); + K[I]--;K[J]--;Y[J]++; + I=J-1; + F++; + } + } + if(To&&To!="P"){ + for(V=[],J=0;J0; Y[J]--) V=cons(1,V); + } + V=reverse(V); + if(To!=0&&To!="01") V=getCatalan(V,0|to=To); + return V; + } + R=qsort(R); + return R; + } + if(!isint(F=getopt(opt))) F=0; + if(!isint(X)){ + if(type(X)==4&&type(car(X))==4){ /* ptg */ + N=length(X)+3; + V=newvect(N);R=newvect(N); + for(TX=X;TX!=[];TX=cdr(TX)){ + V[car(TX)[0]]++;R[car(TX)[1]]++; + } + if(To=="#"){ + for(I=0;I0;J--) K=K+")"; + for(J=V[I];J>0;J--) K=K+"("; + if(++IN) R+=catalan([K-N-1,K-M]); + M++; + }else N++; + } + return R; + } + if(!isint(X)||X++<0) return 0; + /* integer: */ + if(!N){ + for(Y=N=1;X>Y;N++) Y*=(4*N+2)/(N+2); + }else{ + Y=catalan(N); + if(X>Y) return 0; + } + if(F){ + X--; + if(N<3){ + if(N==2) R=X>0?"0011":"0101"; + else if(N==1) R="01"; + else R=""; + } + else for(I=0;I=V) X-=V; + else{ + J=X%catalan(N-I-1); + K=(X-J)/catalan(N-I-1); + R=(I==0)?"01":"0"+getCatalan(K,I|opt=F+1)+"1"; + if(N-I>1) R=R+getCatalan(J,N-I-1|opt=F+1); + break; + } + } + if(To=="s"||F>1) return R; + R=calc(strtoascii(R),["-",48]); + }else{ + for(R=[],M=N;M>0||N>0;){ + Z=Y*(M-N)*(M+1)/(M-N+1)/(M+N); + if(X>Z){ + N--;X-=Z;Y-=Z;R=cons(0,R); + }else{ + M--;Y=Z;R=cons(1,R); + } + } + R=reverse(R); + } + if(To=="s") R=asciitostr(calc(R,["+",48])); + else if(To=="T"||To=="#"||To=="P") R=getCatalan(R,0|to=To); + return R; +} + +def xypg2tg(K) +{ + D=3.1416/2;Or=[0,0];Op="red";Every="";M=0.5;V=0.15;W=0.2;Num=St=Pr=F=0;Line=R=[]; + if(isint(T=getopt(pg))) S=T; + if(isint(T=getopt(skip))) F=T; + if(type(T=getopt(r))==1) M=T; + else if(type(T)==4){ + M=T[0]; + if(length(T)>1) V=T[1]; + if(length(T)>2) W=T[2]; + } + if(isint(T=getopt(proc))) Pr=T; + if(type(T=getopt(org))==4) Or=T; + if(type(T=getopt(rot))==1||T==0) D=T; + if(type(T=getopt(dviout))==1) Dvi=T; + if(type(T=getopt(num))==1) Num=T; + if(type(T=getopt(every))==7) Every=T; + + if(type(car(K)[0])==4){ + if(type(T=getopt(line))==4) Line=T; + S=length(K); + Opt=delopt(getopt(),["Opt","skip","proc","dviout","num","line"]); + if(type(car(Or))!=4||length(Or)!=S){ + Or0=[0,0]; Or1=[1.5,0]; Or2=[0,1.5]; M=10; + if(car(Or)==0&&type(Or[1])==4){ + Or0=Or[1]; + Or=cdr(cdr(Or)); + } + if(length(Or)>1&&type(Or[1])==4){ + M=Or[0]; Or1=Or[1]; + } + if(length(Or)>2) Or2=Or[3]; + for(R=[],I=0;I0){ + Tb=str_tb("%%\n",Tb); + if(type(car(Line))!=4) Line=[Line]; + } + for(S="";Line!=[]; Line=cdr(Line)){ + T=car(Line); + if(length(T)>2){ + S=T[2]; + if(S!="") S="["+S+"]"; + } + Tb=str_tb("\\draw"+S+"(S"+rtostr(T[0])+")--(S"+rtostr(T[1])+");\n",Tb); + } + S=str_tb(0,Tb); + if(Dvi==1) xyproc(S|dviout=1); + else if(Dvi==-1) S=xyproc(S); + return S; + } + } + + if(type(L=getopt(V))>3){ + if(type(L)==4) L=ltov(L); + S=length(L); + }else{ + S=length(K)+3; + L=newvect(S); + } + if(Pr==1){ + if(!L[0]) + for(I=0;I2){ + if(I<=0) Tb=str_tb(";\n",Tb); + TOp="["+car(T)[2]+"]";I=-1; + }else TOp=Op; + if(!I) Tb=str_tb("\\draw "+TOp,Tb); + Tb=str_tb((I%M2)?" ":"\n",Tb); + if(!iand(F,256)) + Tb=str_tb("($(S)+"+ car(L[car(T)[0]]) +"$)--($(S)+" +car(L[car(T)[1]]) +"$)",Tb); + else + Tb=str_tb(car(L[car(T)[0]])+"--"+car(L[car(T)[1]]),Tb); + } + Tb=str_tb(";\n",Tb); + } + if(iand(F,32)) for(I=0;I [I-F,J-F] +* 1 : list of circulate +* 2 : sorted above list +* 3 : minimum in the above +* 4 : minimum mirror +* 5 : minimam extend +* 6 : extend + "std" : to normal form + "#" : #lines (10) + "-#" : inverse of "#" (11) + "del": reduction points + F=[I,J] => another diagonal (flip option) + F=[I] : the other ends of diagonal starting from I + ["ext",I] : [I-1,I] に頂点を付加 + ["del",I] : Iを潰す + ["pair",I] : Iを通る対角線の他方 + ["pairb",I] : Iとつながった頂点（辺でのつながりを含む） + ["mirror",K]：[I,J] -> [K-I,K-J] + ["flip0",[I,J]] : [I,J] の他方の対角線 + ["flip",[I,J]] : [I,J] でフリップ + ["res",I] : + */ +def pgpart(K,F) +{ + S=length(K)+3; + if(type(F)==4){ + if(length(F)==1){ + F=car(F); + for(R=[];K!=[];K=cdr(K)){ + if(car(K)[0]==F) R=cons(car(K)[1],R); + else if(car(K)[1]==F) R=cons(car(K)[0],R); + } + return R; + } + if(length(F)==2){ + if(isint(F[0])){ + F=sort2(F); + K0=pgpart(K,["pair",F[0]]);K0=cons((F[0]+1)%S,K0);K0=cons((F[0]+S-1)%S,K0); + K1=pgpart(K,["pair",F[1]]);K1=cons((F[1]+1)%S,K1);K1=cons((F[1]+S-1)%S,K1); + if(findin(F[1],K0)<0) return []; + R=lsort(K1,K2,"cap"); + if(length(R)!=2) return []; + R=sort2(R); + if(getopt(flip)==1){ + for(RR=[R];K!=[];K=cdr(K)) + RR=cons((F==car(K))?R:car(K),RR); + R=pgpart(RR,0); + } + return R; + } + if(F[0]=="ext"){ + if(F[1]=="all"){ + for(I=0,R=[];I=F)I++; + J=car(K)[1];if(J>=F)J++; + R=cons([I,J],R); + } + return pgpart(R,0); + } + if(F[0]=="res"){ + F1=F[1]; + T=sort2([(F1-1)%S,(F1+1)%S]); + for(R=[];K!=[];K=cdr(K)){ + if(car(K)==T) continue; + if((I=car(K)[0])>F1)I--; + if((J=car(K)[1])>F1)J--; + R=cons([I,J],R); + } + if(length(R)!=S-4) return 0; + return pgpart(R,0); + } + if(F[0]=="pair"||F[0]=="pairb"){ + for(R=[];K!=[];K=cdr(K)){ + if(car(K)[0]==F[1]) R=cons(car(K)[1],R); + if(car(K)[1]==F[1]) R=cons(car(K)[0],R); + } + if(F[0]=="pairb"){ + R=cons((F[1]+1)%S,R);R=cons((F[1]-1)%S,R); + } + return qsort(R); + } + if(F[0]=="flip"||F[0]=="flip0"){ + S=sort2(F[1]); + I=pgpart(K,["pairb",S[0]]);J=pgpart(K,["pairb",S[1]]); + I=lsort(I,J,"cap"); + if(length(I)!=2) return 0; + if(F[0]=="flip0") return I; + for(R=[];K!=[];K=cdr(K)) + R=cons((car(K)==S)?I:car(K),R); + return qsort(R); + } + if(F[0]=="mirror"){ + for(R=[];K!=[];K=cdr(K)) + R=cons([(F[1]-car(K)[0])%S,(F[1]-car(K)[1])%S],R); + return pgpart(R,0); + } + } + } + if(F=="check"){ + K=pgpart(K,0); + for(;K!=[];){ + L=pgpart(K,"res"); + if(L==[]) return 0; + for(L=reverse(L);L!=[];L=cdr(L)){ + R=pgpart(K,["res",car(L)]); + if(R==0) return 0; + K=R; + } + } + return 1; + } + if(F=="std") F=0; + if(type(F)==7){ + S0=[7,8,12,0,13];S1=["#","-#","res","std","cat"]; + I=findin(F,S1); + if(I>=0) F=S0[I]; + } + if(isint(F) && F<=0){ + for(R=[];K!=[];K=cdr(K)){ + I=(car(K)[0]-F)%S; + J=(car(K)[1]-F)%S; + R=cons(sort2([I,J]),R); + } + return qsort(R); + } + if(F>0&&F<4){ + for(R=[],I=0;I1) R=lsort(R,[],1); + if(F==3) R=R[0]; + return R; + } + if(F==4){ + for(R=[];K!=[];K=cdr(K)){ + I=S-car(K)[0]-1; + J=S-car(K)[1]-1; + R=cons([J,I],R); + } + return pgpart(R,3); + } + if(F==5){ + K=pgpart(K,1); + for(R=[];K!=[];K=cdr(K)){ + TK=cons([0,S-1],car(K)); + R=cons(pgpart(TK,3),R); + } + return lsort(R,[],1); + } + if(F==6){ + K=cons([0,S-1],K); + return lsort(pgpart(K,2),[],1); + } + if(F==7||F=="#"){ +/* + for(R=newvect(S);K!=[];K=cdr(K)){ + R[car(K)[0]]++; + R[car(K)[1]]++; + } + return vtol(R); +*/ + return vtol(getCatalan(K,0|to="#")); + } + if(F=="-#"||F==8){ + if(type(K)==4) K=ltov(K); + return getCatalan(K,0|to="P"); + } + if(F==10||F==11){ + S=length(K); + K=ltov(K);L=newvect(S); + for(R=[],T=S-3;T>0;T--){ + for(I=0;I0;J--) if(K[T1=(I+J)%S]) break; + if(T1==T0||T0==I||T1==I) return []; + K[T0]--;K[T1]--;L[I]--; + R=cons([T1,T0],R); + break; + } + if(I==S) return []; + } + if(F==11) return reverse(pgpart(R,8)); + return pgpart(R,0); + } + if(F==12||F=="res"){ + K=pgpart(K,7); + for(I=0,R=[];K!=[];I++,K=cdr(K)) if(!car(K)) R=cons(I,R); + return reverse(R); + } + if(F==13||F==14||F=="0"||F=="("||(F=="T"&&type(K)==7)){ + ST=(F==13||F=="0")?48:40; + S=length(K)+3; + J=newvect(S);I=newvect(S);RR=newvect(S); + for(;K!=[];K=cdr(K)){ + I[car(K)[0]]++; + J[car(K)[1]]++; + } + J[S-1]++; + for(R=[],K=S-1;K>1;K--){ + for(T=J[K];T>0;T--) R=cons(ST+1,R); + for(T=I[K-2];T>0;T--) R=cons(ST,R); + } + R=cons(ST,R); + if(F!="T") return asciitostr(R); + F=="TT"; + } + if(F==9){ + for(R=[];K!=[];K=cdr(K)){ + I=S-car(K)[0]-1; + J=S-car(K)[1]-1; + R=cons([J,I],R); + } + T=pgpart(R,3); + if(imod(S,1))return T; + for(R=[];K!=[];K=cdr(K)){ + I=(-car(K)[0])%S; + J=(-car(K)[1])%S; + R=cons([J,I],R); + } + R=pgpart(R,3); + return T=0;N--){ + T=(Bis==2)?getCatalan(N,K-2|to="P",opt=2):getCatalan(N,K-2|to="P",opt=1); + if(Zig&&length(pgpart(T,"res"))!=2) continue; + if(Al==0||Al==2){ + for(I=1;T!=0&&I3){ + while(K-- > 3) R=pg2tgb(R|verb=F,red=M,all=Al); + return R; + }else if(K<-3){ + for(RR=[],K=-K-3;K>0;K--) RR=cons(R=pg2tgb(R|verb=F,red=M,all=Al),RR); + return reverse(RR); + } + return []; + } + if(K==[]) return (Al==1)?[[[0,2]],[[1,3]]]:[[[0,2]]]; + S=length(car(K))+3; + for(R=[],I=N=0;K!=[];K=cdr(K),I++){ + TR=pgpart(car(K),(Al==1)?6:5); + if(!Al){ + TR=append(pgpart(pgpart(car(K),4),5),TR); + for(T=TR,TR=[];T!=[];T=cdr(T)) if(pgpart(car(T),4) >= car(T)) TR=cons(car(T),TR); + /* 4 => 9 */ + TR=reverse(TR); + } + N+=length(TR); + R=append(TR,R); + if(N>M){ + R=lsort(R,[],1); + M=length(R); + if(F) mycat([M,N]); + N=0; + } + } + R=lsort(R,[],1); + if(F) mycat([length(R),N]); + return R; +} + +def n2a(T) +{ + Opt=[40,41];M=61; + if(type(U=getopt(opt))==7){ + Opt=strtoascii(U); + } + if(!isint(S=getopt(s))) S=0; + if(isint(N=getopt(m))&&N>8&&N<62) M=N; + if(T>M){ + TR=[Opt[1]]; + TR=append(strtoascii(rtostr(T)),TR); + TR=cons(Opt[0],TR); + if(S==1) TR=asciitostr(TR); + return TR; + } + if(T<10) T+=48; + else if(T<36) T+=87; + else if(T<62) T+=29; + if(S) T=[T]; + if(S==1) T=asciitostr(T); + return T; +} + def scale(L) { T=F=0;LS=1; @@ -1614,6 +2572,30 @@ def mydet(M) } } +def permanent(M) +{ + SS=size(M); + if((S=SS[0]) != SS[1] || S==0) return 0; + if((Red=getopt(red))!=1){ + MM = matrtop(M); + for(Dn = 1, I = 0; I < S; I++) + Dn *= MM[1][I]; + return (!Dn)?0:red(permanent(MM[0]|red=1)/Dn); + } + if(S<3){ + if(S==1) return M[0][0]; + else return M[0][0]*M[1][1]+M[0][1]*M[1][0]; + } + LL=m2ll(M); + for(V=I=0;I=0) continue; + N0=M0-diagm(S,[E0]); + for(TR1=R1;TR1!=[];TR1=cdr(TR1)){ + E1=car(TR1)[1]; + if(findin(zz,vars(E1))>=0) continue; + N=newbmat(2,1,[[N0],[M1-diagm(S,[E1])]]); + L=mykernel(N|opt=1); + if(length(L)>0) V=cons([[E0,E1], L],V); + } + } + } + if(type(M)==6){ + M=mtranspose(M); + S=size(M)[0]; + R=meigen(M|mult=1); + for(TR=R;TR!=[];TR=cdr(TR)){ + E=car(TR)[1]; + if(findin(zz,vars(E))>=0) continue; + N=M-diagm(S,[E]); + V=cons([E,mykernel(N)],V); + } + } + V=reverse(V); + if(getopt(TeX)==1||getopt(dviout)==1){ + Sp0="&:"; + if(type(Sp=getopt(sep))!=7){ + if(type(Sp)==4){ + Sp0=Sp[0];Sp=Sp[1]; + }else Sp="\\\\\n"; + } + for(S="",TV=V;TV!=[];TV=cdr(TV)){ + S+=my_tex_form(car(TV)[0])+Sp0+mtotex(mtranspose(lv2m(car(TV)[1]))); + if(length(TV)>1) S+=Sp; + } + if(getopt(dviout)==1) dviout(texbegin("align*",S)); + return S; + } + return V; + } F = getopt(mult); if(type(M)==4 || type(M)==5){ II=length(M); @@ -2106,8 +3135,77 @@ def meigen(M) S = size(M)[0]; P = mydet2(mgen(S,0,[zz],0)-M); return (F==1)?getroot(P,zz|mult=1):getroot(P,zz); +} + +def lext2(L) +{ + if(length(L)==2){ + for(S=0,I=1;IL[0]) break; + } + return [I-1,L[0]+L[1]-S]; + } + if(L[0]==L[1]) return [0,0]; + if(L[0]=6 && Mt!=0)||(L==3&&Mt==1)){ + for(SS=2,I=3; I GRS */ - G=s2sp(M|std=1); + G=M; L=length(G); for(V=[],I=L-2;I>=0;I--) V=cons(makev([I+10]),V); V=cons(makev([L+9]),V); @@ -2265,17 +3398,21 @@ def mmc(M,X) if(Mt!=1) Mt=0; if(R[2]!=2 || R[3]!=0 || !(R=getbygrs(G,1|mat=1))) return 0; MZ=newmat(1,1); - SS=length(G); - if(Mt==1) SS=SS*(SS-1)/2; + SS=length(G)-1; + if(Mt==1) SS=SS*(SS+1)/2; for(M=[],I=0;ISS){ /* addition */ - for(I=0;I=SS){ /* addition */ + for(I=0;I=0; I--){ if(I==J){ @@ -2316,17 +3458,17 @@ def mmc(M,X) if(J==0) M1=MM[0]; else M1=radd(M1,MM[J]); } - /* middle convolution */ + /* convolution of KZ */ for(P=0,Q=1;J=0; I--){ for(RR=[],K=SS-1;K>=0;K--){ MT=MZ; if(I==K){ MT=N[J]; - if(I==P) MT-=N[Q]; - else if(I==Q) MT-=N[P]; - }else if(I==P && K==Q) MT=N[Q]; - else if(I==Q && K==P) MT=N[P]; + if(I==P) MT+=N[Q]; + else if(I==Q) MT+=N[P]; + }else if(I==P && K==Q) MT=-N[Q]; + else if(I==Q && K==P) MT=-N[P]; RR=cons(MT,RR); } R=cons(RR,R); @@ -2335,7 +3477,9 @@ def mmc(M,X) if(++Q==SS){ P++;Q=P+1; } - } + } + if(getopt(homog)==1) MM[L-1]-=diagm(S*SS,[X]); + /* middle convolution */ for(R=[],I=SS-1; I>=0; I--){ for(RR=[N[I]],J=0; J 1 */ if(J>0) M[0][J]= red(M[0][J]/P); - if(Tr) GR[0][J]=red(GR[0][J]/P); + if(Tr) GC[0][J]=red(GC[0][J]/P); } if(S0>1 && S1>1) N=newmat(S0-1,S1-1); else N=0; @@ -2912,24 +4060,30 @@ def mdsimplify(L) return L; } +#if 1 def m2mc(M,X) -{ +{ if(type(M)<2){ mycat([ "m2mc(m,t) or m2mc(m,[t,s])\t Calculation of Pfaff system of two variables\n", " m : list of 5 residue mat. or GRS/spc for rigid 4 singular points\n", -" t : [a0,ay,a1,c], swap, GRS, GRSC, sp, irreducible, pair, pairs, Pfaff, All\n", +" t : [a0,ay,a1,c], swap, GRS, GRSC, extend (eigen), sp, irreducible, pair, pairs, Pfaff, All\n", " s : TeX, dviout, GRSC\n", " option : swap, small, simplify, operator, int\n", " Ex: m2mc(\"21,21,21,21\",\"All\")\n" ]); return 0; } - if(type(M)==7) M=s2sp(M); + if(type(M)==7) M=s2sp(M); + if(type(M)==4&&type(M[0])==6&&length(M)==10){ + N=newvect(5); + N[0]=M[1];N[1]=M[0];N[2]=M[2];N[3]=M[4];N[4]=M[5]; + M=vtol(N); + } if(type(X)==7) X=[X]; Simp=getopt(simplify); if(Simp!=0 && type(Simp)!=1) Simp=2; - Small=(getopt(small)==1)?1:0; + Small=(getopt(small)==1)?1:0; if(type(M[0])==4){ if(type(M[0][0])==1){ /* spectral type */ XX=getopt(dep); @@ -2966,7 +4120,7 @@ def m2mc(M,X) if(type(X)==4 && type(X[0])==7) return m2mc(N,X|keep=Keep,small=Small); return N; - } + } if(type(X)==4 && type(X[0])==7){ Keep=(getopt(keep)==1)?1:0; if(X[0]=="All"){ @@ -3001,9 +4155,61 @@ def m2mc(M,X) if(length(X)>1){ if(X[1]=="dviout") Show=2; if(X[1]=="TeX") Show=1; - } - if(X[0]=="GRS"||X[0]=="GRSC"||X[0]=="sp"){ + } + if(X[0]=="GRS"||X[0]=="GRSC"||X[0]=="sp"||X[0]=="extend"){ Y=radd(-M[0],-M[1]-M[2]); + if(X[0]=="extend"){ + R=[M[1],M[0],M[2],Y, M[3],M[4],radd(-M[1],-M[3]-M[4]), + radd(Y,-M[3]-M[4]),radd(M[1],M[2]+M[4]), radd(M[0],M[1]+M[3])]; + if(length(X)>1){ + if(X[1]=="eigen"){ + L=["x,y","x,0","x,1","x,\\infty","y,0","y,1","y,\\infty", + "0,1","0,\\infty","1,\\infty"]; + U="\\\\\n";S=""; + for(TR=R,TL=L;TR!=[];TL=cdr(TL),TR=cdr(TR)){ + if(length(TL)==1) U="\n"; + S+="A_{"+car(TL)+"}&\\rightarrow\ " + +meigen(car(TR)|vec=1,TeX=1,sep=[":","\\quad"])+U; + } + return S; + } + if(type(X[1])==4){ + TL=[x,y,0,1,2]; + for(T=TL,TO=[];T!=[];T=cdr(TL)){ + J=findin(TL[J],X[1][0]); + if(J<0) return 0; + TO=cons(X[1][1][J],TO); + } + TO=reverse(TO); + for(T=[],I=0;I<4;I++) + for(J=I+1;J<5;J++) T=cons([TL[I],TL[J]],T); + T=reverse(T); + for(R=[],I=0;I<4;I++){ + for(J=I+1;J<5;J++){ + K=findin([TO[I],T0[J]],T); + if(K<0) K=findin([TO[J],T0[I]],T); + R=cons(R[K],R); + } + } + return reverse(R); + } + } + if(Show){ + TL=["x","y","0","1","\\infty"]; + for(S="",TR=R,I=0,J=1;;TR=cdr(TR)){ + S+="A_{"+TL[I]+","+TL[J]+"}&="+mtotex(car(TR)); + if(length(X)>2&&X[2]=="spt") + S+="&&"+ltotex(meigen(car(TR)|mult=1)|vec=1,TeX=1,opt="spt"); + if(J++==4){ + if(I++==3) break; + J=I+1; + } + S+="\\\\\n"; + } + return S; + } + return R; + } if(X[0]!="GRSC"){ L=meigen([M[0],M[1],M[2],M[3],M[4],Y,radd(-M[1],-M[3]-M[4]),radd(Y,-M[3]-M[4])]|mult=1); if(X[0]=="sp"){ @@ -3016,7 +4222,7 @@ def m2mc(M,X) }else{ L=meigen([M[0],M[1],M[2],M[3],M[4],Y,radd(-M[1],-M[3]-M[4]),radd(Y,-M[3]-M[4]), radd(M[0],M[1]+M[3]),radd(M[1],M[2]+M[4])]|mult=1); - S="x=0&x=y&x=1&y=0&y=1&x=\\infty&y=\\infty&x=y=\\infty&x=y=0&x=y=1\\\\\n"; + S="x:0&x:y&x:1&y:0&y:1&x:\\infty&y:\\infty&0:1&1:\\infty&0:\\infty\\\\\n"; } T=ltotex(L|opt="GRS",pre=S,small=Small); if(Show==2) dviout(T|eq=0,keep=Keep); @@ -3073,8 +4279,10 @@ def m2mc(M,X) if(getopt(swap)==1) return m2mc(m2mc(m2mc(M,"swap"),X),"swap"); N=newvect(5); - for(I=0;I<5;I++) - N[I]=M[I]; + for(I=0;I<5;I++){ + if(type(T=M[I])<4) T=diagm(1,[T]); + N[I]=T; + } S=size(N[0])[0]; if(type(X)==4){ for(I=0;I<3;I++){ @@ -3090,18 +4298,213 @@ def m2mc(M,X) MM[0] = newbmat(3,3, [[N[0]+ME,N[1],N[2]], [MZ], [MZ]]); MM[1] = newbmat(3,3, [[MZ], [N[0],N[1]+ME,N[2]], [MZ]]); MM[2] = newbmat(3,3, [[MZ], [MZ], [N[0],N[1],N[2]+ME]]); - MM[3] = newbmat(3,3, [[N[3]+N[1],-N[1]], [-N[0],radd(N[0],N[3])], [MZ,MZ,N[3]]]); + MM[3] = newbmat(3,3, [[N[3]+N[1],-N[1]], [-N[0],radd(N[0],N[3])], [MZ,MZ,N[3]]]); MM[4] = newbmat(3,3, [[N[4]], [MZ,N[4]+N[2],-N[2]], [MZ,-N[1],radd(N[4],N[1])]]); M0 = newbmat(3,3, [[N[0]], [MZ,N[1]], [MZ,MZ,N[2]]]); M1 = radd(MM[0],MM[1]+MM[2]); KE = append(mykernel(M0|opt=1),mykernel(M1|opt=1)); if(length(KE) == 0) return MM; KK = mtoupper(lv2m(KE),0); - for(I=0;I<5;I++) + for(I=0;I<5;I++) MM[I] = mmod(MM[I],KK); + if(Simp!=0) MM = mdsimplify(MM|type=Simp); + return MM; +} +#else +def m2mc(M,X) +{ + if(type(M)<2){ + mycat([ +"m2mc(m,t) or m2mc(m,[t,s])\t Calculation of Pfaff system of two variables\n", +" m : list of 5 residue mat. or GRS/spc for rigid 4 singular points\n", +" t : [a0,ay,a1,c], swap, GRS, GRSC, sp, irreducible, pair, pairs, Pfaff, All\n", +" s : TeX, dviout, GRSC\n", +" option : swap, small, simplify, operator, int\n", +" Ex: m2mc(\"21,21,21,21\",\"All\")\n" +]); + return 0; + } + if(type(M)==7) M=s2sp(M); + if(type(X)==7) X=[X]; + Simp=getopt(simplify); + if(Simp!=0 && type(Simp)!=1) Simp=2; + Small=(getopt(small)==1)?1:0; + if(type(M[0])==4){ + if(type(M[0][0])==1){ /* spectral type */ + XX=getopt(dep); + if(type(XX)!=4 || type(XX[0])>1) XX=[1,length(M[0])]; + M=sp2grs(M,[d,a,b,c],[XX[0],XX[1],-2]|mat=1); + if(XX[0]>1 && XX[1]<2) XX=[XX[0],2]; + if(getopt(int)!=0){ + T=M[XX[0]-1][XX[1]-1][1]; + for(V=vars(T);V!=[];V=cdr(V)){ + F=coef(T,1,car(V)); + if(type(F)==1 && dn(F)>1) + M = subst(M,car(V),dn(F)*car(V)); + } + } + V=vars(M); + if(findin(d1,V)>=0 && findin(d2,V)<0 && findin(d3,V)<0) + M=subst(M,d1,d); + } + RC=chkspt(M|mat=1); + if(RC[2] != 2 || RC[3] != 0){ /* rigidity idx and Fuchs cond */ + erno(0);return 0; + } + R=getbygrs(M,1|mat=1); + if(getopt(anal)==1) return R; /* called by mc2grs() */ + Z=newmat(1,1,[[0]]); + N=[Z,Z,Z,Z,Z,Z]; + for(RR=R; RR!=[]; RR=cdr(RR)){ + RT=car(RR)[0]; + if(type(RT)==4){ + if(RT[0]!=0) N=m2mc(N,RT[0]|simplify=Simp); + N=m2mc(N,[RT[1],RT[2],RT[3]]|simplify=Simp); + } + } + if(type(X)==4 && type(X[0])==7) + return m2mc(N,X|keep=Keep,small=Small); + return N; + } + if(type(X)==4 && type(X[0])==7){ + Keep=(getopt(keep)==1)?1:0; + if(X[0]=="All"){ + dviout("Riemann scheme"|keep=1); + m2mc(M,[(findin("GRSC",X)>=0)?"GRSC":"GRS","dviout"]|keep=1); + dviout("Spectral types : "|keep=1); + m2mc(M,["sp","dviout"]|keep=1); + dviout("\\\\\nBy the decompositions"|keep=1); + R=m2mc(M,["pairs","dviout"]|keep=1); + for(R0=R1=[],I=1; R!=[]; I++, R=cdr(R)){ + for(S=0,RR=car(R)[1][0];RR!=[]; RR=cdr(RR)) S+=RR[0]; + if(S==0) R0=cons(I,R0); + else if(S<0) R1=cons(I,R1); + } + S="irreducibility\\ $"+((length(R0)==0)?"\\Leftrightarrow":"\\Leftarrow") + +"\\ \\emptyset=\\mathbb Z\\cap$"; + dviout(S|keep=1); + m2mc(M,["irreducible","dviout"]|keep=1); + if(R0!=[]) + dviout(ltotex(reverse(R0))|eq=0,keep=1, + title="The following conditions may not be necessary for the irreducibility."); + if(R1!=[]) + dviout(ltotex(reverse(R1))|eq=0,keep=1,title="The following conditions can be omitted."); + if(getopt(operator)!=0){ + dviout("The equation in a Pfaff form is"|keep=1); + m2mc(M,["Pfaff","dviout"]|keep=Keep,small=Small); + } + else if(Keep!=1) dviout(" "); + return M; + } + Show=0; + if(length(X)>1){ + if(X[1]=="dviout") Show=2; + if(X[1]=="TeX") Show=1; + } + if(X[0]=="GRS"||X[0]=="GRSC"||X[0]=="sp"||X[0]=="extend"){ + Y=radd(-M[0],-M[1]-M[2]); + if(X[0]=="extend") + return [M[1],M[0],M[2],Y, M[3],M[4],radd(-M[1],-M[3]-M[4]), + radd(Y,-M[3]-M[4]),radd(M[1],M[2]+M[4]), radd(M[0],M[1]+M[3])]; + if(X[0]!="GRSC"){ + L=meigen([M[0],M[1],M[2],M[3],M[4],Y,radd(-M[1],-M[3]-M[4]),radd(Y,-M[3]-M[4])]|mult=1); + if(X[0]=="sp"){ + L=chkspt(L|opt="sp"); + V=[L[1],L[0],L[2],L[5]]; W=[L[1],L[3],L[4],L[6]]; + if(Show==2) dviout(s2sp(V)+" : "+s2sp(W)|keep=Keep); + return [V,W]; + } + S="x=0&x=y&x=1&y=0&y=1&x=\\infty&y=\\infty&x=y=\\infty\\\\\n"; + }else{ + L=meigen([M[0],M[1],M[2],M[3],M[4],Y,radd(-M[1],-M[3]-M[4]),radd(Y,-M[3]-M[4]), + radd(M[0],M[1]+M[3]),radd(M[1],M[2]+M[4])]|mult=1); + S="x=0&x=y&x=1&y=0&y=1&x=\\infty&y=\\infty&x=y=\\infty&x=y=0&x=y=1\\\\\n"; + } + T=ltotex(L|opt="GRS",pre=S,small=Small); + if(Show==2) dviout(T|eq=0,keep=Keep); + if(Show==1) L=T; + return L; + } + if(X[0]=="Pfaff"){ + S=ltotex(M|opt=["Pfaff",u,x,x-y,x-1,y,y-1],small=Small); + if(Show==2) dviout(S|eq=0,keep=Keep); + return S; + } + if(X[0]=="irreducible"){ + L=meigen([M[0],M[1],M[2],radd(-M[0],-M[1]-M[2])]|mult=1); + S=getbygrs(L,10|mat=1); + if(Show==2) dviout(ltotex(S)|eq=0,keep=Keep); + return S; + } + if(X[0]=="pairs"||X[0]=="pair"){ + L=meigen([M[0],M[1],M[2],radd(-M[0],-M[1]-M[2])]|mult=1); + S=chkspt(L|opt=0); + V=(Show==2)?1:0; + S=sproot(L,X[0]|dviout=V,keep=Keep); + return S; + } + if(X[0]=="swap"){ + Swap=getopt(swap); + if(type(Swap)<1 || Swap==1) + return newvect(6,[M[3],M[1],M[4],M[0],M[2],M[5]]); + if(Swap==2) + return newvect(5,[radd(M[0],M[1]+M[3]),M[4],M[2],radd(-M[1],-M[3]-M[4]),M[1]]); + if(type(Swap)==4 && length(Swap)==3){ + MX=radd(-M[0],-M[1]-M[2]); MY=radd(-M[3],-M[1]-M[4]); + if(Swap[0]==1){ + MX0=M[2];MY0=M[4]; + } + else if(Swap[0]==2){ + MX0=MX;MY0=MY; + }else{ + MX0=M[0];MY0=M[3]; + } + if(Swap[1]==1){ + MX1=M[2];MY1=M[4]; + } + else if(Swap[1]==2){ + MX1=MX;MY1=MY; + }else{ + MX1=M[0];MY1=M[3]; + } + return newvect(5,MX0,M[1],MX1,MY0,MY1); + } + } + return 0; + } + if(getopt(swap)==1) + return m2mc(m2mc(m2mc(M,"swap"),X),"swap"); + N=newvect(6); + for(I=0;I<6;I++) + N[I]=M[I]; + S=size(N[0])[0]; + if(type(X)==4){ + for(I=0;I<3;I++){ + if(X[I] != 0) + N[I] = radd(N[I],X[I]); + } + if(length(X)==3) return N; + X=X[3]; + } + MZ = newmat(S,S); + ME = mgen(S,0,[X],0); + MM = newvect(6); + MM[0] = newbmat(3,3, [[N[0]+ME,N[1],N[2]], [MZ], [MZ]]); /* A01 */ + MM[1] = newbmat(3,3, [[MZ], [N[0],N[1]+ME,N[2]], [MZ]]); /* A02 */ + MM[2] = newbmat(3,3, [[MZ], [MZ], [N[0],N[1],N[2]+ME]]); /* A03 */ + MM[3] = newbmat(3,3, [[N[3]+N[1],-N[1]], [-N[0],radd(N[0],N[3])], [MZ,MZ,N[3]]]); /* A12 */ + MM[4] = newbmat(3,3, [[N[4]], [MZ,N[4]+N[2],-N[2]], [MZ,-N[1],radd(N[4],N[1])]]); /* A23 */ + MM[5] = newbmat(3,3, [[MZ,N[5]+N[2],-N[2]], [N[5]], [MZ,-N[0],radd(N[5],N[0])]]); /* A13 */ + M0 = newbmat(3,3, [[N[0]], [MZ,N[1]], [MZ,MZ,N[2]]]); + M1 = radd(MM[0],MM[1]+MM[2]); + KE = append(mykernel(M0|opt=1),mykernel(M1|opt=1)); + if(length(KE) == 0) return MM; + KK = mtoupper(lv2m(KE),0); + for(I=0;I<6;I++) MM[I] = mmod(MM[I],KK); if(Simp!=0) MM = mdsimplify(MM|type=Simp); return MM; } +#endif def easierpol(P,X) { @@ -3336,11 +4739,19 @@ def llbase(VV,L) T = length(L); for(I = 0; I < S; I++) V[I] = nm(red(V[I])); - LV = 0; + LV = 0; for(J = 0; J < T; J++){ X = var(L[J]); N = deg(L[J],X); - for(I = LV; I < S; I++){ + for(I = LV; I < S; I++){ if((C2=coef(V[I],N,X)) != 0){ + if(type(C2)==1){ + for(K=I+1;Kabs(C1)){ + I=K;C2=C1; + } + } + } if(I > LV){ Temp = V[I]; V[I] = V[LV]; @@ -3350,15 +4761,19 @@ def llbase(VV,L) if(I == LV || (C1 = coef(V[I],N,X)) == 0) continue; Gcd = gcd(C1,C2); - V[I] = V[I]*tdiv(C2,Gcd)-V[LV]*tdiv(C1,Gcd); + V[I] = V[I]*tdiv(C2,Gcd)-V[LV]*tdiv(C1,Gcd); } LV++; } - } + } } return V; } +def rev(A,B){return A>B?-1:(A=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); @@ -3413,9 +4875,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); } @@ -3771,24 +5232,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) @@ -3865,9 +5351,19 @@ 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); + if(type(M)==4){ + if(length(M)==1) + N=1-(M=car(M)); + else{ + N=M[0];M=M[1]; + } + }else N=1; + return vtol(N*X+M*Y); } def mrot(X) @@ -4151,21 +5647,21 @@ def texsp(P) def fctrtos(P) { /* extern TeXLim; */ - if(!chkfun("write_to_tb", "names.rr")) return 0; TeX = getopt(TeX); if(TeX != 1 && TeX != 2 && TeX != 3) TeX = 0; - if((Dvi=getopt(dviout)==1) && TeX<2) TeX=3; + if((Dvi=getopt(dviout)==1) && TeX<2) TeX=3; if(TeX>0){ Lim=getopt(lim); if(Lim!=0 && TeX>1 && (type(Lim)!=1||Lim<30)) Lim=TeXLim; else if(type(Lim)!=1) Lim=0; CR=(TeX==2)?"\\\\\n":"\\\\\n&"; - if(TeX==1 || Lim==0) CR=""; - else if((Pages=getopt(pages))==1) CR="\\allowdisplaybreaks"+CR; + CR2="\\allowdisplaybreaks"+CR; + if(TeX==1 || Lim==0) CR=CR2=""; + else if((Pages=getopt(pages))==1) CR2=CR; if(!chkfun("print_tex_form", "names.rr")) return 0; Small=getopt(small); @@ -4242,7 +5738,10 @@ def fctrtos(P) } VV=reverse(VV);VD=reverse(VD); Rev=(getopt(rev)==1)?1:0; - Dic=(getopt(dic)==1)?1:0; + Rdic=0; + if((Dic=getopt(dic))==2){ + Dic=Rdic=1; + }else if(Dic!=1) Dic=0; TT=terms(P,VV|rev=Rev,dic=Dic); if(TeX==0){ Pre="("; Post=")"; @@ -4250,7 +5749,7 @@ def fctrtos(P) Pre="{"; Post="}"; } Out = string_to_tb(""); - for(L=C=0,Tm=TT;Tm!=[];C++,Tm=cdr(Tm)){ + for(L=C=CC=0,Tm=TT;Tm!=[];C++,Tm=cdr(Tm)){ for(I=0,PC=P,T=cdr(car(Tm)),PW="";T!=[];T=cdr(T),I++){ PC=mycoef(PC,D=car(T),VV[I]); if(PC==0) continue; @@ -4262,7 +5761,8 @@ def fctrtos(P) else PT="^"+rtostr(D); } if(Dif>0) PW+=(Dif==1)?"d":"\\partial "; - PW+=VD[I]+PT; + if(Rdic) PW=VD[I]+PT+PW; + else PW+=VD[I]+PT; } } D=car(Tm)[0]; @@ -4271,10 +5771,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 = ""; @@ -4296,16 +5798,28 @@ def fctrtos(P) } } if(Lim>0){ + CC++; LL=texlen(OC)+texlen(PW); if(LL+L>=Lim){ if(L>0) str_tb(CR,Out); if(LL>Lim){ - if(TOC==7) OC=texlim(OC,Lim|cut=CR); - PW+=CR; L=0; + if(TOC==7) OC=texlim(OC,Lim|cut=[CR,CR2]); + if(length(Tm)!=1) PW+=CR; + L=0; }else L=LL; }else L+=LL; - }else if(length(Tm)!=1) PW += CR; /* not final term */ - if(TeX) OC=texsp(OC); + }else if(length(Tm)!=1){ + CC++; + PW += CR; /* not final term */ + } + if(CC>TeXPages) CR=CR2; + 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); @@ -4335,14 +5849,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); @@ -4350,6 +5864,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) @@ -4360,19 +5875,23 @@ def fctrtos(P) if(nmono(P[J][0])>1|| (!isvar(P[J][0])||vtype(P[J][0]))&&str_len(SS)>1) SS="("+SS+")"; write_to_tb(SS,S); - str_tb(["^", (TeX>1)?rtotex(P[J][1]):monotos(P[J][1])],S); + str_tb(["^", (TeX>=1)?rtotex(P[J][1]):monotos(P[J][1])],S); }else{ - if(nmono(P[J][0])>1) SS="("+SS+")"; + if(nmono(P[J][0])>1&&length(P)>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); } } S = str_tb(0,S); - if((Lim>0 || TP!=2) && CR!="") S=texlim(S,Lim|cut=CR); + if((Lim>0 || TP!=2) && CR!="") S=texlim(S,Lim|cut=[CR,CR2]); } if(TeX>0){ if(Small==1) S=str_subst(S,"\\frac{","\\tfrac{"); if(Dvi==1){ - dviout(strip(S,"(",")")|eq=(Pages==1)?6:0); S=1; + dviout(strip(S,"(",")")|eq=(Pages==1||Pages==2)?6:0); S=1; } } return S; @@ -4396,7 +5915,12 @@ def texlim(S,Lim) mycat(["Set TeXLim =",Lim]); return 1; } - if(type(Out=getopt(cut))!=7) Out="\\\\\n&"; + if(type(Out=getopt(cut))!=7){ + if(type(Out)!=4) Out=Out2="\\\\\n&"; + else{ + Out2=Out[1];Out=Out[0]; + } + } if(type(Del=getopt(del))!=7) Del=Out; if(Lim<30) Lim=TeXLim; S=ltov(strtoascii(S)); @@ -4427,6 +5951,7 @@ def texlim(S,Lim) SS=str_tb(0,0); L=cons(length(S),L); L=reverse(L); + if(length(L)>TeXPages) Out=Out2; for(I=0; L!=[]; I=J,L=cdr(L)){ str_tb((I==0)?"":Out,SS); J=car(L); @@ -4574,6 +6099,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; @@ -4605,6 +6159,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; @@ -4857,6 +6429,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]); @@ -4906,7 +6484,10 @@ def execdraw(L,P) } } if(MM) V=ptaffine(MM,V|option_list=Org); - if(length(T)>3) V=append(V,T[3]); + if(length(T)>3){ + if(type(T2=T[3])==7) T2=[T2]; + V=append(V,T2); + } str_tb(xyput(V),Out); }else if(T[0]==3){ F++; @@ -4936,9 +6517,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]); } @@ -5105,6 +6692,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,I,DX))) continue; + P=red(P-TP*DX^I+TP*muldo(DX^(I-1),F,L)); + } + return P; + } if(type(F) <= 3){ if(type(L) == 4 && type(L[0]) == 4) return applpdo(P,F,L); @@ -5145,6 +6741,26 @@ def appledo(P,F,L) #endif } +def caldo(P,L) +{ + for(R=0;P!=[];P=cdr(P)){ + TP=car(P); + if(type(TP)<4){ + R=red(R+TP);continue; + } + for(S=1;TP!=[];TP=cdr(TP)){ + S0=car(TP); + if(type(S0)==4){ + TP0=S0; + for(S0=1,K=TP0[1];K>0;K--) S0=muldo(S0,TP0[0],L); + } + S=muldo(S,S0,L); + } + R=red(R+S); + } + return R; +} + def muldo(P,Q,L) { if(type(Lim=getopt(lim))!=1) Lim=100; @@ -5269,13 +6885,38 @@ def mce(P,L,V,R) { L = vweyl(L); X = L[0]; DX = L[1]; - P = sftexp(laplace1(P,L),L,V,R); + P=red(P); + if(findin(DX,dn(P))>=0) return 0; + PP=fctr(nm(P)); + for(P=1;PP!=[];PP=cdr(PP)){ + TP=car(PP); + if(findin(DX,vars(TP[0]))>=0) P*=TP[0]^TP[1]; + } + 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 mcme(P,L,V,R) +{ + for(LL=[];L!=[];L=cdr(L)) LL=cons(car(L),LL); + LL=reverse(LL); + if(V==0) L=LL; + else L=delopt(LL,V|inv=1); + P=rede(P,LL); + for(Q=laplace(P,L),E=0,TL=L;TL!=[];TL=cdr(TL)){ + E+=TL[0]*TL[1]; + Q=toeul(Q,car(TL),0); + } + N=length(L); + for(R=[],TL=L;TL!=[];TL=cdr(TL)) + R=cons([-E/car(TL),-TL[0]*TL[1]+(R-1)/N],R); + Q=transpdo(Q,L,reverse(R)|ex=1); + return rede(Q,LL); } def rede(P,L) @@ -5449,9 +7090,30 @@ def mulpdo(P,Q,L); } } #endif + +def transppow(LL,M) +{ + for(L=[];LL!=[];LL=cdr(LL)) + L=cons(vweyl(car(LL)),L); + L=reverse(L);N=length(L); + if(type(M)==4) M=lv2m(M); + MM=myinv(M); + for(K=[],I=0;I3) return +#ifdef USEMODULE + mtransbys(os_md.transpdosub,P,[LL,K]); +#else + mtransbys(transpdosub,P,[LL,K]); +#endif Len = length(K)-1; if(Len < 0 || P == 0) return P; @@ -5476,13 +7138,17 @@ def transpdosub(P,LL,K) } def transpdo(P,LL,K) -{ - if(type(K[0]) < 4) - K = [K]; +{ + if(type(K)==4&&type(K[0]==4)){ + for(TK=K;;TK=cdr(TK)) if(!isint(car(TK))) break; + if(TK==[]) K=transppow(LL,K); + } + if(type(K)==6) K=transppow(LL,K); Len = length(K)-1; K1=K2=[]; if(type(LL)!=4) LL=[LL]; if(type(LL[0])!=4) LL=[LL]; + if(type(car(K)) < 4 && length(LL)!=length(K)) K = [K]; if(getopt(ex)==1){ for(LT=LL, KT=K; KT!=[]; LT=cdr(LT), KT=cdr(KT)){ L = vweyl(LL[J]); @@ -5491,10 +7157,25 @@ def transpdo(P,LL,K) } K2=append(K1,K2); }else{ + if(length(LL)==length(K) && type(car(K))!=4){ + for(DV=V=TL=[],J=length(LL)-1;J>=0;J--){ + TL=cons(vweyl(LL[J]),TL); + V=cons(car(TL)[0],V); + DV=cons(car(TL)[1],DV); + } + LL=TL; + if(type(RK=solveEq(K,V|inv=1))!=4) return TK; + if(!isint(Inv=getopt(inv))) Inv=0; + if(iand(Inv,1)){J=K;K=RK;RK=J;} + M=jacobian(RK,V|mat=1); + M=mulsubst(M,[V,K]|lpair=1); + RK=vtol(M*ltov(DV)); + if(Inv>1) return RK; + K=lpair(K,RK); + } for(J = length(K)-1; J >= 0; J--){ L = vweyl(LL[J]); - if(L[0] != K[J][0]) - K1 = cons([L[0],K[J][0]],K1); + if(L[0]!= K[J][0]) K1=cons([L[0],K[J][0]],K1); K2 = cons(K[J][1],K2); } P = mulsubst(P, K1); @@ -5878,7 +7559,7 @@ def divdo(P,Q,L) } P -= muldo(SR*(DX)^(J-I),Q,L); S += SR*(DX)^(J-I); - } + } return [S,P,M]; } @@ -6241,10 +7922,12 @@ def baseODE(L) if(type(F=getopt(f))!=1) F=0; if(isint(In=getopt(in))!=1) In=0; if(type(Ord=getopt(ord))!=1&&Ord!=0) Ord=2; + Pages=getopt(pages); + if(Pages!=1&&Pages!=2) Pages=0; if(Ord>3){ Ord-=4; Hgr=1; }else Hgr=0; - if(type(car(L))==4&&type(L[1])==7){ + if(type(car(L0=L))==4&&type(L[1])==7){ Tt=L[1];L=car(L); } M=N=length(L); SV=SVORG; @@ -6264,13 +7947,48 @@ def baseODE(L) for(Var=[],I=M-1;I>=0;I--) Var=cons(makev([SV[I]]),Var); } if(type(To=getopt(to))<2||type(To)>4) To=0; - else if(!isvar(To)){ - if(type(To)!=4) To=cons(red(To),cdr(Var)); + if(Ord<0){ /* cancell y1, z1,... by baseODE0() */ + if(Ord==-1) Ord=2; + if(type(To)==4||!isvar(To)){ + L=L0=baseODE(L0|to=To,f=-3)[1]; + To=0; + } + R=baseODE0(L|option_list= + delopt(getopt(),[["var",Var],["ord",Ord]]|inv=1)); + if(TeX){ + if(type(R)==4&&length(R)>1&&type(R[1])==4) R=R[1]; + if(type(To)==2 && !isvar(To)){ + S0=baseODE(L0|TeX=1,f=-1,to=To); + V=baseODE0(L|step=-1,to=To); + }else{ + S0=baseODE(L0|TeX=1,f=-1); + V=baseODE0(L|step=-1,to=To); + } + T=eqs2tex(R,[V,2,Pages]); + S=((F==1)?(Tt+"\n"):S0)+texbegin("align*",T); + if(TeX==2) dviout(S); + return S; + } + return R; + } + if(To&&!isvar(To)){ + if(type(To)!=4){ + To=red(To); + for(K=0;K0&&I<2&&J<2) break; + } + if(K==length(Var)) return -9; + J=To; + for(To=[],I=length(Var)-1;I>=0;I--) + if(I!=K) To=cons(Var[I],To); + To=cons(J,To); + } if(type(To)==4){ if(type(car(To))==4){ R=1;To=car(To); }else R=0; - if(type(IL=solveEq(To,Var|inv=1))!=4) return -7; + if(type(IL=solveEq(To,Var|inv=1))!=4) return IL; if(R==1){ R=To;To=IL;IL=R; } @@ -6295,8 +8013,8 @@ def baseODE(L) } } } - if(F==-3) return [Var,L]; - for(I=0;I=0 && IT=0&&!chkfun("gr",0)){ mycat("load(\"gr\"); /* <- do! */\n"); @@ -6366,7 +8105,7 @@ def baseODE(L) } if(F==-2) return [VV,L]; if(F<0) return [V,L]; - LL=(Hgr==1)?hgr(L,V,Ord):gr(L,V,Ord); + LL=(Hgr==1)?hgr(L,V,Ord):gr(L,V,Ord); if(F==2) return [V,L,LL]; if(Ord==2) P=LL[0]; else{ @@ -6376,19 +8115,279 @@ def baseODE(L) LL=tolex_tl(LL,V,Ord,V,2);P=LL[0]; } } - V0=makev([car(SV),M]); - CP=mycoef(P,mydeg(P,V0),V0); - if(cmpsimple(-CP,CP)<0) P=-P; if(TeX){ - for(V0=[makev([car(SV)])],I=1;I<=M;I++) V0=cons(makev([car(SV),I]),V0); - T="&\\!\\!\\!"+fctrtos(P|var=VV,dic=1,TeX=3); - S=((F==1)?(Tt+"\n"):S0)+texbegin("align*",texbegin("split",T)); + for(V0=[],I=1;I<=M;I++) V0=cons(makev([car(SV),I]),V0); + T=eqs2tex(P,[V0,2,Pages]); + if(!Vl||Vl==1){ + for(I=1,K=0;I1 && isvar(L[1])) L=[L]; + R=car(L);L=cdr(L);Sgn=1; + }else R=[]; + if(type(R)==4&&car(R)==0){ + Sgn=0;R=cdr(R); + } + if(L!=[]){ + Dic=car(L);L=cdr(L); + } + if(L!=[]){ + Pages=car(L);L=cdr(L); + } + if(L!=[]) Cont=car(L); + if(type(P)==4){ + for(S="";P!=[];P=cdr(P)){ + S+=eqs2tex(car(P),[R,Dic,Pages,Cont]); + if(!Cont) Cont=1; + } +/* S=str_subst(S,"\\\\&,\\\\",",\\\\&"); */ + if(getopt(dviout)==1) dviout(S|eq=6); + return S; + } + if(type(R)==2) R=[R]; + if(Sgn){ + for(;R!=[];R=cdr(R)) + if((Deg=mydeg(P,car(R)))>0) break; + if(Deg>0){ + CP=mycoef(P,Deg,car(R)); + if(cmpsimple(-CP,CP)<0) P=-P; + } + } + S="&\\!\\!\\!"; + if(Cont) + S=(Pages?",\\allowdisplaybreaks":",")+"\\\\\n"+S; + S+=fctrtos(P|var=R,dic=Dic,TeX=3,pages=Pages); + if(getopt(dviout)==1) dviout(S|eq=6); + return S; +} + +/* Opt: var, opt, dbg */ +def res0(P,Q,X) +{ + if(!isvar(X)){ + if(!isvar(P)) return -1; + Y=P;P=Q;Q=X;X=Y; + } + if(isvar(Var=getopt(var))) Var=[Var]; + else if(type(Var)!=4) Var=0; + if(type(W=getopt(w))!=4) W=[]; + if(!isint(Opt=getopt(opt))&&type(Opt)!=4) Opt=0; + if(type(Dbg=getopt(dbg))==4){ + Fct=Dbg[1];Dbg=Dbg[0]; + } + if(!isint(Dbg)) Dbg=0; + P=nm(P);Q=nm(Q); + Fctr=isfctr(P)*isfctr(Q); + DP=deg(P,X);DQ=deg(Q,X); + if(DP==DQ&&nmono(coef(P,DP,X))0){ + if(DP=2) mycat([DP,"(",nmono(P), nmono(coef(P,DP,X)),") :", + DQ, "(",nmono(Q),nmono(coef(Q,DQ,X)), ")"]); + else mycat0([DP,":",DQ,","],0); + } + TQ=coef(Q,DQ,X);TP=coef(P,DP,X); + if(Fctr){ + T=gcd(TP,TQ);M=red(TQ/T); + if(Var&&M!=car(W)&&type(TV=vars(M))==4&&lsort(TV,Var,2)!=[]) W=cons(M,W); + P=M*(P-coef(P,DP,X)*X^DP)-red(TP/T)*X^(DP-DQ)*(Q-coef(Q,DQ,X)*X^DQ); + if(Var){ +#if 1 + if(Dbg>2) mycat0(">",0); + for(S=SS=fctr(P),P=1,C=0;S!=[];S=cdr(S)){ + TV=vars(S0=car(S)[0]); + if(type(TV)==4&&lsort(TV,Var,2)!=[]){ + for(TW=W;TW!=[];TW=cdr(TW)){ + if(gcd(car(TW),S0)!=1){ + S0=1;break; + } + } + if(Dbg>1){ + if(S0==1) mycat(["Reduced by :",nmono(car(TW))]); + else if(C++>0){ + mycat(["Product :", nmono(P), nmono(S0)]); + if(Dbg==3){ + if(!Fct||Fct==[]){ + if(C>1) P=1; + }else{ + if(car(Fct)==C){ + C=10000;Fct=cdr(Fct);P=1; + }else S0=1; + } + }else if(Dbg==4) return [SS,Q,DP,DQ,W]; + } + } + P*=S0; + } + } +#else + for(TW=W;TW!=[];TW=cdr(TW)){ + if((C=gcd(P,car(TW)))!=1){ + P=red(P/C); + if(Dbg>=2&&nmono(Q)>1) mycat(["Reduce :",nmono(C)]); + } + } +#endif + } + }else{ + if(type(TQ)==1){ + Q/=TQ; + P=P-TP*X^(DP-DQ)*Q; + }else P=TQ*P-TP*X^(DP-DQ); + if(deg(P,X)==DP) P-=coef(P,DP,X)*X^DP; + } + DP=deg(P,X); + if(Opt==-2||(type(Opt)==4&&Opt[0]==DP&&Opt[1]==DQ)) return [P,Q,DP,DQ,W]; + } + if(Dbg){ + if(Dbg>1) mycat([DP,"(",nmono(P), nmono(coef(P,DP,X)),") :", + DQ, "(",nmono(Q), nmono(coef(Q,DQ,X)), ")"]); + else mycat0([DP,":",DQ," "],0); + } + if(Opt==1) Q=[P,Q,DP,DQ,W]; + return (DQ==0)?Q:0; +} + +/* Opt : f, var, ord, ord, step, f, to */ +def baseODE0(L) +{ + if(!isint(Ord=getopt(ord))) Ord=-1; + if(Ord==-1) Ord=2; + if(Ord0&&Ord>0) Ord=-1; + N=length(L); + if(type(To=getopt(to))==4&&length(To)==N){ + V=cdr(To);To=car(To); + } + if(!isvar(To)) To=V=0; + if(type(SV=Var=getopt(var))!=4){ + SV=SVORG; + if(N>10){ + R=[]; + for(K=N-1;K>9;K++) R=cons(SV[floor(K/10)-1]+SV[K%10],R); + SV=append(SV,R); + } + for(Var=[],I=N-1;I>=0;I--) Var=cons(makev([SV[I]]),Var); + } + if((J=findin(To,Var))>0){ + TV=TL=[]; + for(I=N-1;I>=0;I--){ + if(I!=J){ + TV=cons(Var[I],TV);TL=cons(L[I],TL); + } + } + Var=cons(Var[J],TV);L=cons(L[J],TL); + } + if(!To) To=car(SV); + Q=car(L); + V0=makev([To,1]); + R=[V0-Q];V0=[V0]; + for(I=2;I<=N;I++){ + P=diff(t,Q); + if(type(P)==3) P=red(P); + for(TV=Var,TL=L;TV!=[];TV=cdr(TV),TL=cdr(TL)){ + P+=diff(Q,car(TV))*car(TL); + if(type(P)==3) P=red(P); + } + Q=P; + TV=makev([To,I]); + R=cons(nm(TV-Q),R); + V0=cons(TV,V0); + } + if(Step==-1) return V0; + if(!V) V=cdr(Var); + if(Ord<0){ + for(C=1,R0=[];V!=[];V=cdr(V),C++){ + TR=R=reverse(R); + if(length(R)>1){ /* reduce common factor */ + P=car(TR);TR=cdr(TR); + for(;TR!=[]&&P!=1;TR=cdr(TR)) + P=gcd(P,car(TR)); + if(P!=1){ + for(TR=[];R!=[];R=cdr(R)) TR=cons(red(car(R)/P),TR); + R=reverse(TR); + } + } + TR=[]; + TV=car(V); + if(length(V)==1) V0=[car(V0)]; + if(C==Step) return [append(V,V0),R]; + while(R!=[]&&findin(TV,vars(car(R)))<0){ + TR=cons(car(R),TR); + R=cdr(R); + } + R0=(F==2)?append(R,R0):cons(car(R),R0); + if(R!=[]){ + for(W=[],P=car(R),R=cdr(R); R!=[]; R=cdr(R)){ + if(Dbg) mycat0(["\nStep ",C,"-",length(R)," ",TV, + (type(Dbg)==4||Dbg>=2)?"\n":" "],0); + if(findin(TV,vars(car(R)))<0){ + TR=cons(car(R),TR); + continue; + } + if(Ord>-3){ + if(Dstep&&Dstep[0]==C&&Dstep[1]==length(R)) + return res0(P,car(R),TV|var=V0,opt=cdr(cdr(Dstep)),dbg=Dbg); + else TQ=res0(P,car(R),TV|var=V0,opt=1,dbg=Dbg,w=W); + if(Dbg==4&&type(car(TQ))==4) return TQ; + if(Ord==-2) P=car(TQ); + W=TQ[4];TQ=TQ[1]; + }else{ + TQ=res(TV,P,car(R)); + Q=fctr(TQ); /* irreducible one */ + for(TQ=1;Q!=[];Q=cdr(Q)) + if(lsort(V0,vars(car(Q)[0]),2)!=[]) TQ*=car(Q)[0]; + } + TR=cons(TQ,TR); + } + } + R=TR; + } + if(Dbg==1) mycat([]); + return (F==1)?car(R):(F==2?append(R,R0):cons(car(R),R0)); + } + V=append(V,[makev([To,N])]); + if(Step==1) return [R,V]; + R=gr(R,V,Ord); + return (F==1)?car(R):R; /* hgr(R,V,Ord); */ +} + + def taylorODE(D){ Dif=(getopt(dif)==1)?1:0; if(D==0) return Dif?f:f_00; @@ -6479,37 +8478,40 @@ def taylorODE(D){ } return R; } + +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); + if(II==I) S=II-J; + else if(P!=0 && II-J>S) S=II-J; + } + F *= X^S; + } + 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); + if(getopt(raw)!=1){ + 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; + } + for(R = 0 ; I >= 0; I--) + R += (red(mycoef(F,I,DX)/X^I))*DX^I; + return pol2sft(R,DX); +} + -def toeul(F,L,V) -{ - L = vweyl(L); - X = L[0]; DX = L[1]; - I = mydeg(F,DX); - if(V == "infty"){ - for(II=I; II>=0; II--){ - J = mydeg(P=mycoef(F,I,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--) - 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--) - R += (red(mycoef(F,I,DX)/X^I))*DX^I; - return pol2sft(R,DX); -} - /* def topoldif(P,F,L) { @@ -6541,8 +8543,11 @@ 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){ + R=nm(R); + while(mycoef(R,0,X) == 0) + R = tdiv(R,X); + } if(V != "infty" && V != 0) R = mysubst(R,[X,X-V]); return R; @@ -6551,11 +8556,10 @@ 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=getopt()),[DX,DX+N]); + return fromeul(P,L,V|option_list=getopt()); } - def fractrans(P,L,N0,N1,N2) { L = vweyl(L); @@ -6789,7 +8793,7 @@ def getroot(F,X) C2=mycoef(P,2,X);C1=mycoef(P,1,X);C0=mycoef(P,0,X); C=sqrt2rat(C1^2-4*C0*C2); C0=[]; - if(type(C)==0&&ntype(C)==0&&pari(issquare,-C)) C0=sqrt(C); + if(type(C)==0&&ntype(C)==0&&issquare(-C)) C0=sqrt(C); else if(Cpx>1) C0=sqrtrat(C); if(C0==[]&&Cpx>2) C0=C^(1/2); if(C0!=[]){ @@ -6947,6 +8951,10 @@ def expat(F,L,V) def polbyroot(P,X) { + if(isvar(V=getopt(var))&&length(P)>1&&isint(car(P))){ + for(Q=[],I=car(P);I<=P[1];I++) Q=cons(makev([V,I]),Q); + P=Q; + } R = 1; while(length(P)){ R *= X-car(P); @@ -7135,6 +9143,30 @@ 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); + Y=getopt(var); + Abs=(Y==1)?1:0; + if(!(K=length(V))) return Y==1?1: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|var=Abs):(Y==1?1:abs(Q)))*X^D; + } + if(isvar(Y)) for(;V!=[];V=cdr(V)) R=subst(R,car(V),Y); + return R; +} + def prehombf(P,Q) { if((Mem=getopt(mem))!=1 && Mem!=-1) @@ -7534,9 +9566,10 @@ def stoe(M,L,N) L = vweyl(L); Size = size(M); S = Size[0]; - NN = 0; + NN = -1; if(type(N) == 4){ NN=N[0]; N=N[1]; + if(N==NN) return 1; }else if(N < 0){ NN=-N; N=0; } @@ -7546,7 +9579,7 @@ def stoe(M,L,N) MN = dupmat(M); MD = newmat(S,S); DD = D[0]; - DD[N] = 1; DD[S] = 1; + DD[N]=1; DD[S] = 1; for(Lcm = I = 1; ; ){ DD = D[I]; MM = MN[N]; @@ -7559,9 +9592,9 @@ def stoe(M,L,N) DD[J] = red(DD[J]*Lcm); if(I++ >= S) break; - if(I==S && NN>0){ + if(I==S && NN>=0){ DD = D[I]; - DD[0]=-z_zz; DD[NN]=1; + DD[S]=z_zz; DD[NN]=1; break; } Mm = dupmat(MN*M); @@ -7579,7 +9612,7 @@ def stoe(M,L,N) if(mydeg(P[I][0],L[1]) > 0) R *= P[I][0]^P[I][1]; } - if(NN > 0) + if(NN >= 0) R = -red(coef(R,0,z_zz)/coef(R,1,z_zz)); return R; } @@ -7783,7 +9816,6 @@ def okuboetos(P,L) Phi[J+1] = Phi[J]*(X-L[J]); for(ATT = AT[N], J = 0; J < N; J++) ATT[J] = mycoef(P,J,DX); - for(K = 1; K <= N; K++){ for(J = N; J >= K; J--){ Aj = A[J-1]; @@ -7804,7 +9836,6 @@ def okuboetos(P,L) ATj[J-K-1] = red(ATj[J-K-1]-DAT); } } - ATT = newmat(N,N); for(J = 0; J < N; J++){ for(K = 0; K < N; K++){ @@ -7985,34 +10016,36 @@ def sgn(X) def calc(X,L) { - if(type(X)<4||type(X)==7){ - if(type(L)==4||type(L)==7){ - V=L[1]; - if(type(X)!=7){ - if((L0=L[0])=="+") X+=V; - else if(L0=="-") X-=V; - else if(L0=="*") X*=V; - else if(L0=="/") X/=V; - else if(L0=="^") X^=V; - } - if((L0=L[0])==">") X=(X>V); - else if(L0=="<") X=(X=") X=(X>=V); - else if(L0=="<=") X=(X<=V); - else if(L0=="!=") X=(X!=V); - }else if(type(L)==7&&type(X)<4){ - if(L=="neg") X=-X; - else if(L=="abs") X=abs(X); - else if(L=="neg") X=-X; - else if(L=="sqr") X*=X; - else if(L=="inv") X=1/X; - else if(L=="sgn"){ - if(X>0)X=1; - else if(X<0) X=-1; - } + if((T=type(X))==4||T==5) return map(os_md.calc,X,L); + if(type(L)==4){ + V=L[1]; + if((L0=L[0])==">") X=(X>V); + else if(L0=="<") X=(X=") X=(X>=V); + else if(L0=="<=") X=(X<=V); + else if(L0=="!=") X=(X!=V); + else if(type(X)==6 || type(X)<4){ + if((L0=L[0])=="+") X+=V; + else if(L0=="-") X-=V; + else if(L0=="*") X*=V; + else if(L0=="/") X/=V; + else if(L0=="^") X^=V; } + return X; } + if(type(L)!=7||T>7||T==4||T==5) return X; + if(L=="neg") X=-X; + else if(L=="sqr") X*=X; + else if(L=="inv"){ + if(T==6) X=myinv(X); + else X=1/X; + }else if(T==6) return X; + if(L=="abs") X=abs(X); + else if(L=="sgn"){ + if(X>0) X=1; + else if(X<0) X=-1; + } return X; } @@ -8272,156 +10305,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); } @@ -8523,6 +10415,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;JJ[1]) J=[J[1],J[0]]; L=lsort(L,J,1); - if(length(L)!=1) return 0; + if(length(L)!=1) return 0; return [mc2grs(G,["get0",I]),mc2grs(G,["get0",[I[0],J[0]],[I[1],J[1]]]), mc2grs(G,["get0",[I[0],J[1]],[I[1],J[0]]]),mc2grs(G,["get0",[I[0],I[1],L[0]]])]; } @@ -9408,29 +11493,36 @@ def mc2grs(G,P) } for(LL=[],T=L;T!=[];T=cdr(T)) LL=cons(mc2grs(G,["get0",[I,car(T)]]),LL); - LL=reverse(LL); - for(R=[],Q=mc2grs(G,["get0",I]);Q!=[];Q=cdr(Q)){ + LL=reverse(LL); +/* mycat(["I",I,"\nL",L,"\nLL",LL]); +mycat([I,mc2grs(G,["get0",I])]); */ + for(R=[],Q=mc2grs(G,["get0",I]);Q!=[];Q=cdr(Q)){ /* Q : simultaneous spct */ +/* mycat(["Q",Q,car(Q)[1]]); */ for(T=[],J=2;J>=0;J--){ - V=anal2sp(LL[J],["get1",(I[0]128){ - 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; } @@ -10578,6 +12693,132 @@ def evalma(S) return S; } +def evalcoord(L) +{ + if(type(L)==7) L=strtoascii(L); + I=str_str(L,"("); + if(I>=0) J=str_pair(L,I+1,"(",")"); + if(I<0 || J32&&(C<40||C>58)){F=0;break;} + } + S0=str_cut(L,I+1,J-1); + for(;J>=0;J--) L=cdr(L); + while(L!=[]&&car(L)<33) L=cdr(L); + if(F){ + S="["+S0+"]"; + return [eval_str(S),L]; + }else return [[S0],L]; +} + +def readTikZ(L) +{ + if(type(L)!=4) L=strtoascii(L); + R=[]; + CMD=["draw","fill","filldraw","shade","shadedraw","clip","pattern","node","begin"]; + while(L!=0&&L!=[]){ + while(L!=[]&&car(L)<33) L=cdr(L); + if(L==[]) break; + if(car(L)==34){ /* % */ + while(L!=[]&&car(L)!=10) L=cdr(L); + continue; + } + if(car(L)!=92) {L=0;break;} /* \ */ + for(DF=0;DF<9;DF++) if(str_str(L,CMD[DF]|top=1,end=1)==1) break; + if(DF<7){ + S=T=0; + I=str_str(L,"(");J=str_str(L,"["); + if(J>0&&I>J){ + K=str_str(L,"]"); + S=str_cut(L,J+1,K-1); + } + F0=F=0;C=[]; + while(L!=0&&L!=[]){ + V=evalcoord(L); + L=V[1]; + if(L==[]) break; + if(F0){ + if (!F) C=cons(0,C); + else if(F0!=3) C=cons(1,C); + } + C=cons(V[0],C); + F0=F;F=0; + if(L[0]==34){ /* % */ + while(L!=[]&&car(L)!=10) L=cdr(L); + continue; + } + if(!str_str(L,"..")){ /* .. */ + L=cdr(L);L=cdr(L); + F=1; + }else if(!str_str(L,"--")){ /* -- */ + L=cdr(L);L=cdr(L); + F=2; + } + while(L!=[]&&car(L)<33) L=cdr(L); + if(L==[]){L=0; break;} + if(!str_str(L,"cycle")){ + if(F==2) C=cons(1,C); + C=cons(-1,C); + F0=F=0; + continue; + } + if(!str_str(L,"and")||!str_str(L,"control")) + F=3; /* control, and */ + else if(car(L)==59){ /* ; */ + L=cdr(L); + break; + }else if(isalpha(car(L))){ + T=[]; + while(car(L)!=40 && car(L)!=59){ /* ( ; */ + T=cons(car(L),T); + if((L=cdr(L))==[]){L=0;break;} + } + T=asciitostr(reverse(T)); + if(car(L)==59){ /* ; */ + L=cdr(L); + break; + } + F0=0;continue; + }else if(F!=1&&F!=2){ + L=0;break; + } + } + if(T){ + if(length(C)==1||length(C)==2) S=(!S)?["",T]:[S,T]; + else{ + L=0;break; + } + } + S=(!S)? []:[["opt",S]]; + if(DF) S=S=cons(["cmd",CMD[DF]],S); + if(T&&length(C)) R=cons((length(C)==1)?[3,S,C[0],DF]:[3,S,C[1],C[0]],R); + else R=cons([1,S,reverse(C)],R); + }else{ /* \node */ + U=0; + I=str_str(L,"(");J=str_str(L,"["); + if(J>0&&I>J){ + K=str_str(L,"]"); + U=str_cut(L,J+1,K-1); + } + V=evalcoord(L); + C=V[0];L=V[1]; + J=str_str(L,"{");K=str_pair(L,J+1,"{","}"); + S=str_cut(L,J+1,K-1); + if(U) S=[U,S]; + R=cons([2,[],C,[S]],R); + for(;K>=0;K--) L=cdr(L); + K=str_str(L,";"); + for(;K>=0;K--) L=cdr(L); + }; + } + if(!L){ + mycat("Can't understand!"); + return -1; + } + return reverse(R); +} + def i2hex(N) { Opt=getopt(); @@ -10776,6 +13017,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}"]; @@ -10850,7 +13092,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); @@ -10871,7 +13113,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; @@ -11028,7 +13279,7 @@ def str_subst(S, L0, L1) def dviout0(L) { - Cmd=["TikZ","TeXLim","TeXEq","DVIOUT","XYPrec","XYcm","XYLim","Canvas"]; + Cmd=["TikZ","TeXLim","TeXEq","DVIOUT","XYPrec","XYcm","XYLim","Canvas","TeXPages"]; if(type(Opt=getopt(opt))==7){ if((F=findin(Opt,Cmd)) < 0) return -1; if(L==-1){ @@ -11041,7 +13292,8 @@ def dviout0(L) if(F==4) V=XYPrec; else if(F==5) V=XYcm; else if(F==6) V=XYLim; - else V=Canvas; + else if(F==7) V=Canvas; + else if(F==8) V=TeXPages; } return V; } @@ -11059,6 +13311,7 @@ def dviout0(L) else if(F==4) XYPrec=L; else if(F==5) XYcm=L; else if(F==6) XYLim=L; + else if(F==8) TeXPages=L; } mycat0([Cmd[F],"=",L],1); return 1; @@ -11093,6 +13346,7 @@ def dviout0(L) mycat0(["DVIOUTL=\"", DVIOUTL,"\""],1); mycat(["Canvas =", Canvas]); mycat(["TeXLim =", TeXLim]); + mycat(["TeXPages =", TeXPages]); mycat(["TeXEq =", TeXEq]); mycat(["AMSTeX =", AMSTeX]); mycat(["TikZ =", TikZ]); @@ -11190,7 +13444,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); @@ -11280,6 +13534,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); @@ -11306,14 +13576,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){ @@ -11331,23 +13614,51 @@ 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(["pages",2],cons(["dviout",1],Opt)));return; }else if(T==4){ + F=0;N=length(getopt()); + if(Raw) N--; + if(N==1){ + if(type(Var=getopt(var))>1){ + if(isvar(Var)) Var=[0,Var]; + else if(type(Var)==4&&Var[0]!=0) Var=cons(0,Var); + else Var=0; + }else if(type(Var=getopt(eqs))!=4) Var=0; + }else if(N==0) Var=[]; + else Var=0; + if(type(Var)==4){ + for(F=0,L=P;L!=[];L=cdr(L)){ /* */ + if(type(car(L))==2) F+=nmono(car(L)); + else{ + F=0;break; + } + } + } + if(F>50){ + S=texbegin("align*",eqs2tex(P,Var)); + if(Raw) return S; + dviout(S);return; + } 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)){ @@ -11375,8 +13686,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)){ @@ -11408,14 +13719,23 @@ def show(P) } } }else if(T==7){ - 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; + if(Var=="raw") S=P+"\n\n"; + else if(Var != "eq" &&str_str(P,"\\begin"|end=128)<0){ + if((TikZ&&str_str(P,"\\draw"|end=128)>=0)||(!TikZ&&str_str(P,"\\ar@"|end=128)>=0)) + S=xyproc(P); + }else if(Var !="eq"){ + if(str_str(P,"\\begin{align")>=0 || str_str(P,"\\[")>=0 + || str_str(P,"\\begin{equation")>=0 + || (str_char(P,0,"^")<0 && str_char(P,0,"_")<0 && str_char(P,0,"&")<0)) + S=P+"\n\n"; } + if(P!=S){ + 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); } @@ -11615,6 +13935,25 @@ def rtotex(P) return (str_len(S) == 1)?S:"{"+S+"}"; } +def togreek(P,T) +{ + R0=[a,b,c,d,e,i,k,l,m,n,o,p,r,s,t,u,x,z]; + R1=[alpha,beta,gamma,delta,epsilon,iota,kappa,lambda, + mu,nu,omega,pi,rho,sigma,theta,tau,xi,zeta]; + if(T==0||T==[]) T=[a,b,c]; + for(S=[],TR=T;TR!=[];TR=cdr(TR)){ + if(type(TR[0])!=4){ + if((I=findin(car(TR),R0))>=0) S=cons([car(TR),R1[I]],S); + }else if((I=findin(car(TR)[0],R0))>=0){ + for(U=car(TR)[1];U!=[];U=cdr(U)) + S=cons([makev([R0[I],car(U)]),makev([R1[I],car(U)])],S); + } + } + if(getopt(raw)==1) return S; + if(getopt(inv)==1) return mysubst(P,S|inv=1); + else return mysubst(P,S); +} + def mtotex(M) { /* extern TexLim; */ @@ -11652,10 +13991,15 @@ def mtotex(M) for(I=0; I1)?fctrtos(P|TeX=2,lim=0,var=Var):fctrtos(P|TeX=2,lim=0); - if(type(P)==1 && str_str(SS[I][J],"\\frac{-"|end=0)==0) - SS[I][J]="-\\frac{"+str_cut(SS[I][J],7,100000); + if(P!=0 || Null == 0 || (Null==2 && I==J)){ + if(type(V=getopt(pfrac))==2) + SS[I][J]=pfrac(P,V|TeX=1); + else{ + SS[I][J]=(type(Var)>1)? + fctrtos(P|TeX=2,lim=0,var=Var):fctrtos(P|TeX=2,lim=0); + if(type(P)==1 && str_str(SS[I][J],"\\frac{-"|end=0)==0) + SS[I][J]="-\\frac{"+str_cut(SS[I][J],7,100000); + } } }else if(type(P)==6){ ST= mtotex(P|small=1,len=1); @@ -11845,6 +14189,86 @@ def frac2n(N) #endif } +/* Option : opt */ +def ptconvex(L) +{ + if(!(isint(Opt=getopt(opt)))) Opt=0; + L0=car(L);X=L0[0];Y=L0[1]; + for(TL=cdr(L);TL!=[];TL=cdr(TL)){ /* find the most left pt L0 */ + if(X0?Y^2:-Y^2)/S,S],L0); + R=cons(L0,R); + } + L=qsort(R); + if(Opt==2) return L; + + for(R=[],TL=L;TL!=[];TL=cdr(TL)){ + if(Opt==4){ + L0=car(TL); + V=car(L0); + L0=append(cdr(cdr(L0)),[V]); + }else L0=cdr(cdr(car(TL))); + R=cons(L0,R); + } + L=reverse(R); + if(Opt==1) return L; + R=[cons(V0=-8,L0=car(L))]; + for(TL=cdr(L);TL!=[];TL=cdr(TL)){ + V=darg(L0,L1=car(TL)); + if(V<-4) continue; + while(V2){ + if((V-=4)>4) return -4; + }else if(V<=-2) V+=4; + return V; + } + X=Q[0]-P[0];Y=Q[1]-P[1]; + if(!(S=X^2+Y^2)) return -8; + V=Y^2/S; + if(Y<0) V=-V; + return X<=0?2-V:V; +} + +def dwinding(P,Q) +{ + V=V0=V1=darg(P,Q0=car(Q)); + Q=cons(Q0,reverse(Q)); + for(Q=cdr(Q);Q!=[];Q=cdr(Q)){ + if((V2=darg(P,car(Q)))<-4) return 1/3; + V1=V2-V1; + if(V1==2||V1==-2) return 1/2; + if(V1<-2) V1+=4; + else if(V1>2) V1-=4; + V+=V1; + V1=V2; + } + return floor((V0-V+1/2)/4); +} + def xyproc(F) { if(type(Opt=getopt(opt))!=7) Opt=""; @@ -11922,6 +14346,9 @@ def xypos(P) def xyput(P) { + if(type(T=car(P))==4||type(car(P)==5)){ + P=cdr(P);P=cons(T[1],P);P=cons(T[0],P); + } if((type(Sc=getopt(scale))==1 && Sc!=1) || type(Sc)==4){ if(type(Sc)==1) Sc=[Sc,Sc]; Sx=Sc[0];Sy=Sc[1]; @@ -11934,6 +14361,31 @@ def xyput(P) return "\\"+xypos(P)+";\n"; } +def xylabel(P,S) +{ + if(type(P[0])==4){ + if(type(S)==4){ + if(type(Put=getopt(put))!=7) Put=""; + if(type(Opt=getopt(opt))!=7) Opt=0; + for(R="";P!=[];P=cdr(P),S=cdr(S)){ + T=car(S); + if(Opt) T=[Opt,T]; + R+=xyput([car(P),Put,T]|option_list=getopt()); + } + return R; + } + if(type(S)==7){ + B=getopt(base); + if(!isint(B)) B=0; + for(SS=[],I=length(P)-1;I>=0;I--) SS=cons(S+rtostr(I+B),SS); + return xylabel(P,SS); + } + } + if(P[0]==0||type(P[0])==1){ + if(type(S)==7) return xylabel([P],[S]|option_list=getopt()); + } +} + def xyline(P,Q) { if(!TikZ) return "{"+xypos(P)+" \\ar@{-} "+xypos(Q)+"};\n"; @@ -12396,14 +14848,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]); @@ -12417,13 +14869,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)); + } + } + if(isint(W=getopt(wait))) sleep(W*100); + }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, @@ -13044,10 +15700,29 @@ 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 dlog10(X) +{ + if(X==0||imag(Z)!=0) return []; + if(X<0) X=-X; + Neg=1; + if(X<1){X=1/X;Neg=-1;} + C=10^(2000); + for(V=0;X>10^(2000);X/=C, V+=2000); + C=10^(200); + for(;X>10^(200);X/=C, V+=200); + V+=nlog(X); + return Neg*V; +} + def mypow(Z,R) { if(type(Z=eval(Z))>1||type(R=eval(R))>1) return todf(os_md.mypow,[Z,R]); @@ -13064,7 +15739,7 @@ def mypow(Z,R) def myarg(Z) { - if(type(Z=map(eval,Z))==4){ + if(type(Z=map(eval,Z))==4||type(Z)==5){ if(length(Z)!=2) return todf(os_md.myarg,[Z]); Re=Z[0];Im=Z[1]; }else if(type(Z)>1){ @@ -13109,17 +15784,29 @@ def arg(Z) if(vars(Z=map(eval,Z))!=[]) return todf(os_md.arg,[Z]); return (type(Z)==4)?pari(arg,Z[0],Z[1]):arg(sqrt,Z); } + +def issquare(X) +{ + if(X==0) return 1; + if(type(X)==1&&ntype(X)==0){ + for(N=nm(X),I=0;I<2;N=dn(X),I++){ + if(isqrt(N)^2!=N) return 0; + } + return 1; + } +} -def sqrt(Z){ +def sqrt(Z) +{ if(vars(Z=map(eval,Z))!=[]) return todf(os_md.sqrt,[Z]); R=(type(Z)==4)?Z[1]:Z; if(ntype(R)==0){ if(R==0) return 0; if(R>0){ - if(pari(issquare,R)) return pari(isqrt,nm(R))/pari(isqrt,dn(R)); + if(issquare(R)) return pari(isqrt,nm(R))/pari(isqrt,dn(R)); }else{ R=-R; - if(pari(issquare,R)) return pari(isqrt,nm(R))/pari(isqrt,dn(R))*@i; + if(issquare(R)) return pari(isqrt,nm(R))/pari(isqrt,dn(R))*@i; } } return (type(Z)==4)?pari(sqrt,Z[0],Z[1]):pari(sqrt,Z); @@ -13814,6 +16501,126 @@ 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),LXm=lmax(L0)]; + S=SX=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),LYm=lmax(L0)]; + S=SY=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); + if((AX=getopt(ax))==1||AX==2){ + if(LY[0]<0&&LY[1]>0){ + glib_line(LX[0],0,LX[1],0); + if(AX==2&&LXm>0){ + E=floor(dlog(LXm)/dlog(10)); + V=floor(LXm*10^(-E)+1/128)*10^E; + glib_line(V,0,V,SY/64); + glib_print(V,-SY/128,rtostr(V)); + } + } + if(LX[0]<0&&LX[1]>0){ + glib_line(0,LY[0],0,LY[1]); + if(AX==2&&LYm>0){ + E=floor(dlog(LYm)/dlog(10)+1/64); + V=floor(LYm*10^(-E)+1/128)*10^E; + glib_line(0,V,SX/64,V); + glib_print(SX/96,V,rtostr(V)); + } + + } + } + 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; @@ -14126,7 +16933,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){ @@ -14364,6 +17171,60 @@ def polroots(L,V) return reverse(SS); } +def lsub(P) +{ + if((T=type(P[0]))==4){ + Q=reverse(P[1]);P=reverse(P[0]); + for(R=[];P!=[];P=cdr(P),Q=cdr(Q)) R=cons(car(Q)-car(P),R); + return R; + }else if(T==5){ + L=length(P[0]);Q=P[1];P=P[0]; + R=newvect(L); + for(V=[],L--;L>=0;L--) R[L]=Q[L]-P[L]; + return R; + } + return P; +} + +def dext(P,Q) +{ + P=lsub(P);Q=lsub(Q); + return P[0]*Q[1]-P[1]*Q[0]; +} + +def ptinversion(P) +{ + if(type(P)==4&&type(P[1])==4){ + for(R=[];P!=[];P=cdr(P)) + R=cons(ptinversion(car(P)|option_list=getopt()),R); + return reverse(R); + } + if(type(V=getopt(org))!=0) V=[0,0]; + if(P==[0,0]) return 0; + if(type(P[0])==4){ + R=P[1];P=P[0]; + } + X=P[0]-V[0];Y=P[1]-V[1];N=X^2+Y^2; + if(type(T=getopt(sc))==1||T==0){ + S=(T<0)?(-T^2):T^2; + }else S=-1; + if(!R){ + if(!N) return 0; + return [X/N+V[0],S*Y/N+V[1]]; + } + N-=R^2; + if(!N){ + if(X+R!=0) X0=X+R; + else X0=X-R; + S=[]; + S=cons(ptinversion([X0,Y]|option_list=getopt()),S); + if(Y+R!=0) Y0=Y+R; + else Y0=Y-R; + return cons(ptinversion([X,Y0]|option_list=getopt()),S); + } + return [[X/N+V[0],S*Y/N+V[1]],T^2*R/N]; +} + def ptcommon(X,Y) { if(length(X)!=2 || length(Y)!=2) return 0; @@ -14389,7 +17250,7 @@ def ptcommon(X,Y) } XX=X[1][0]-X[0][0];YY=X[1][1]-X[0][1]; Arg=(length(Y)<2)?0:Y[1]; - Arg=deval(Arg); + Arg=feval(Arg); if(Arg!=0){ S=dcos(Arg)*XX-dsin(Arg)*YY; YY=dsin(Arg)*XX+dcos(Arg)*YY; @@ -14406,19 +17267,29 @@ def ptcommon(X,Y) T=[Y[0][0]+(Y[1][0]-Y[0][0])*y_-S[0], Y[0][1]+(Y[1][1]-Y[0][1])*y_-S[1]]; R=lsol(T,[x_,y_]); - if(type(R[0])==4&&type(R[1])==4&&R[0][0]==x_&&R[1][0]==y_){ - if(!In || (R[0][1]>=0&&R[0][1]<=1&&R[1][1]>=0&&R[1][1]<=1) ) - return subst(S,x_,R[0][1],y_,R[1][1]); + if(type(R[0])==4&&type(R[1])==4&&R[0][0]==x_&&R[1][0]==y_){ + /* unique sol of parameters */ + if(In && (R[0][1]<0||R[0][1]>1||R[1][1]<0||R[1][1]>1) ) return 0; + return subst(S,x_,R[0][1],y_,R[1][1]); } - if((type(R[0])>0&&type(R[0])<4)||(type(R[1])>0&&type(R[1])<4)) return 0; - if(!In) return 1; - I=(X[0][0]==X[1][0]&&Y[0][0]==Y[1][0]&&X[0][0]==Y[0][0])?1:0; - if(X[0][I]<=X[1][I]){ - X0=X[0][I];X1=X[1][I]; - }else{ - X1=X[0][I];X0=X[1][I]; + if((type(R[0])>0&&type(R[0])<4)||(type(R[1])>0&&type(R[1])<4)) return 0; /* no solution */ + F=0; + if(X[0]==X[1]) F=1; + else if(Y[0]==Y[1]) F=2; + if(!In){ + if(!F) return 1; + else if(F==1) return X[0]; + else if(F==2) return Y[0]; } - return ((Y[0][I]X1&&Y[1][I]>X1))?0:1; + X0=X[0];X1=X[1]; + if(X0>X1){R=X0;X0=X1;X1=R;} + Y0=Y[0];Y1=Y[1]; + if(Y0>Y1){R=Y0;Y0=Y1;Y1=R;} + if(X0Y1) X1=Y1; + if(X0>X1) return 0; + if(X01 && R[0]==R[1]) R=cdr(R); + return R; + } + if(dext([L[0],L[1]],[L[0],L[2]])<0) L=[L[0],L[2],L[1]]; + L=cons(L[2],L); + for(I=F=1;I<4;I++,L=cdr(L)){ + if((V=dext([L[0],L[1]],[L[0],P])) < 0) return 0; + if(!V) F++; + } + return F; +} + def tobezier(L) { if((Div=getopt(div))==1||Div==2){ @@ -14750,21 +17654,39 @@ def areabezier(V) else R=cmpf([V[0],V[2]]); V=[R,V[1],[0,1]]; } - if(type((Int=getopt(int)))==1 && type(V[0])<4 && (V1=V[1])>=0){ + if(type((Int=getopt(int)))==1 && type(V[0])<4 && (V1=V[1])>=0){ + if(Int==3||Int==4){ + R=newvect(V1); + Opt=delopt(getopt(),[["int",7-2*Int]]|inv=1); + for(I=0;I=I;J--) R[J]=(P*R[J]-R[J-1])/(P-1); + } + return R[V1-1]; + } if(Int==2&&iand(V1,1)) V1++; if(!V1) V1=32; - Opt=cons(["raw",1],getopt()); - W=xygraph(V[0],V[1],V[2],Mx,Mx|option_list=Opt); - SS=W[0][1]; + Opt=cons(["raw",1],getopt()); + VV=eval(V[2][1]-V[2][0]); + if(Int==-1) + V=[V[0],V[1]-1,[V[2][0]+VV/V[1]/2,V[2][1]-VV/V[1]/2]]; + W=xygraph(V[0],V[1],V[2],Mx,Mx|option_list=Opt); + SS=W[0][1]; for(S0=S1=0,I=0,L=W;L!=[] && I<=V1;I++, L=cdr(L)){ if(iand(I,1)) S1+=car(L)[1]; - else S0+=car(L)[1]; + else S0+=car(L)[1]; + if(I==3&&(getopt(Acc)==1||(isint(Prec=getopt(prec))&&Prec>10))){ + S0=tobig(S0);S1=tobig(S1); + } if (I==V1) SS+=car(L)[1]; - } - VV=deval(V[2][1]-V[2][0]); + } + if(getopt(Acc)==1) VV=tobig(VV); if(Int==2) return (2*S0+4*S1-SS)*VV/(3*V1); - else + else if(Int==-1){ + return V1==1? S0:(S0+S1)/V1; + }else return (2*S0+2*S1-SS)*VV/(2*V1); } Opt=cons(["opt",0],getopt()); @@ -14811,6 +17733,213 @@ def ptbezier(V,L) return [subst(B,t,L[1]),subst(BB,t,L[1])]; } +/* +def isroot(P,Q,I) +{ + if(subst(P,X,X0=I[0])*subst(P,X,I[1])<=0) return 1; + XM=(I[1]+I[0])/2;W=XM-X0; + if(W<0) W=-W; + X=var(P); + if(!Q) Q=diff(P,X); + Q=subst(Q,X,X+I2);D=deg(Q,X); + for(M=0,P=1,I=deg(Q,X);I<=D;I++){ + V=coef(Q,I,X); + M+=(V<0?-V:V)*P; + P*=W; + } + V=subst(P,X,X0); + if(V<0) V=-V; + return (V-M<=0) 2:0; +} +*/ + +def sgnstrum(L,V) +{ + X=var(car(L)); + if(X==0) X=var(L[1]); + for(F=N=0;L!=[];L=cdr(L)){ + P=car(L); + if(type(V)==7){ + C=coef(P,D=deg(P,X),X); + if(V=="-"&&iand(D,1)) C=-C; + }else C=subst(P,X,V); + if(!C) continue; + if(C*F<0) N++; + F=C; + } + return N; +} + +def polstrum(P) +{ + X=vars(P0=P); + if(!length(X)) return []; + X=car(X); + if(isfctr(P)){ + D=gcd(P,Q=diff(P,X)); + P=sdiv(P,D); + if(getopt(mul)==1&&type(getopt(num))<0) + return append(polstrum(D|mul=1),[P]); + } + D=deg(P,X); + P=P/coef(P,deg(P,X),X); + Q=diff(P,X)/D; + for(L=[Q,P];D>0;){ + R=urem(P,Q); + if((D=deg(R,X))<0) break; + C=coef(R,D,X); + if(C>0) C=-C; + R/=C; + L=cons(R,L); + P=Q;Q=R; + } + if(type(N=getopt(num))>0){ + if(getopt(mul)!=1){ + if(type(N)==1) N=["-","+"]; + return sgnstrum(L,N[0])-sgnstrum(L,N[1]); + } + if(!isfctr(P0)) return -1; + R=polstrum(P0|mul=1); + for(C=0;R!=[];R=cdr(R)) C+=polstrum(car(R)|num=N); + return C; + } + return reverse(L); +} + +def iceil(X) +{ + S=(X>0)?1:-1; + X*=S; + if(X>1) X=ceil(X); + else if(X>1/2) X=1; + else if(X) X=1/floor(1/X); + return S*X; +} + +def polradiusroot(P) +{ + X=var(P);D=deg(P,X); + if(D<1) return -1; + C=coef(P,D,X); + P/=-C; + Int=getopt(int); + if(getopt(comp)==1){ + for(ND=0,TD=0;TD0){ + N2++; + if(!iand(D-TD,1)) N1++; + }else if(iand(D-TD,1)) N1++; + } + for(V1=V2=0,TD=0;TD0){ + TV=eval((C*N2)^(1/(D-TD))); + if(V20){ + TV=eval((C*N1)^(1/(D-TD))); + if(V11&&MC<10001) MC1=MC; + else MC1=MC=32; + if(type(I=getopt(in))!=4){ + I=polradiusroot(P); + W=(I[1]-I[0])/1024; + I=[I[0]-W,I[1]+W]; + } + if(type(L=type(getopt(strum)))!=4) L=polstrum(P); + N0=sgnstrum(L,I[0]);N1=sgnstrum(L,I[1]); + P=car(L);X=var(P); + if(N0<=N1) return []; /* [L,I,N0,N1]; */ + LT=[[0,I[0],I[1],N0,N1]];R=[]; + Z=eval(exp(0)); + while(LT!=[]){ + T=car(LT);LT=cdr(LT); + C=T[0];X0=T[1];X1=T[2];N0=T[3];N1=T[4]; + if(N0<=N1)continue; + if(N0==N1+1){ + V0=subst(P,X,X0); + V1=subst(P,X,X1); + while(C++0&&V2>0)||(V0<0&&V2<0)) X0=X2; + else X1=X2; + } + R=cons([X0,X1,1],R); + continue; + } + while(++CN2){ + if(N2>N1) LT=cons([C,X2,X1,N2,N1],LT); + X1=X2; + N1=N2; + if(N0==N1+1){ + LT=cons([C,X0,X1,N0,N1],LT); + C=MC+1; + } + }else{ + X0=X2; + N0=N2; + } + } + if(C!=MC+2) R=cons([X0,X1,N0-N1],R); + } + if(isint(Nt=getopt(nt)) && Nt>0){ + if(Nt>256) Nt=256; + Q=diff(P,X); + for(S=[],TR=R;TR!=[];TR=cdr(TR)){ + if(car(TR)[2]>1) continue; + V0=subst(P,X,car(TR)[0]); + V1=subst(P,X,car(TR)[1]); + if(abs(V0)0;Tn--){ + X1=X0-V0/subst(Q,X,X0); + V1=subst(P,X,X1); + if(abs(V1)>=abs(V0)) break; + X0=X1;V0=V1; + } + S=cons(X0,S); + } + for(TR=R;TR!=[];TR=cdr(TR)) + if(car(TR)[2]>1) S=cons(car(TR),S); + return reverse(S); + } + return reverse(cons(P,R)); +} + +/* +def ptcombezier0(P,Q) +{ + PB=subst(tobezier(P|div=1),t,s); + QB=tobezier(Q|Div=1); + Z=res(PB[0]-QB[0],PB[1]-QB[1],s); + D=pmaj(diff(Z,t)|val=t); +} +*/ + def ptcombezier(P,Q,T) { if(type(T)<2){ @@ -15019,7 +18148,7 @@ def lbezier(L) else{ if(R!=[]&&F!=0) R=cons(0,R); R=cons(G=car(LT),R); - if(In==3) In==2; + if(In==3) In=2; } for(LT=cdr(LT);LT!=[];LT=cdr(LT)) R=cons(car(LT),R); @@ -15038,7 +18167,7 @@ def lbezier(L) } RT=cons(T,RT); }else if(T==0){ - if(RT==[]) R=cons(reverse(RT),R); + if(RT!=[]) R=cons(reverse(RT),R); RT=[];F=0; }else if(T==1){ if(RT!=[]){ @@ -15060,6 +18189,7 @@ def lbezier(L) def xybezier(L) { + if(type(L)==4&&type(car(L))==4&&type(car(L)[0])==4) L=lbezier(L|inv=1); if(L==0 || (LS=length(L))==0) return ""; Out=str_tb(0,0); if(type(VF=getopt(verb))==4){ @@ -15227,49 +18357,76 @@ def xybox(L) def xyang(S,P,Q,R) { - Opt=getopt(); + Opt=delopt(getopt(),"ar"); + if(type(S)>2) S=dnorm([S,P]); if(type(Prec=getopt(prec))!=1) Prec=0; if(type(Q)>2){ - if(R==1||R==-1){ /* 直角 */ - P1=ptcommon([Q,P],[-S,0]); - S*=R; - P2=ptcommon([P,P1],[S,@pi/2]); - P3=ptcommon([P1,P2],[S,@pi/2]); - return xylines([P1,P2,P3]|option_list=Opt); - }else if((AR=abs(R))==0||AR==2||AR==3||AR==4){ /* 矢印 */ - Ang=myarg([Q[0]-P[0],Q[1]-P[1]]); - if(R<0) Ang+=3.14159; - ANG=[0.7854,0.5236,1.0472]; - X=(AR==0)?1.5708:ANG[AR-2]; - U=[P[0]+S*dcos(Ang+X),P[1]+S*dsin(Ang+X)]; - V=[P[0]+S*dcos(Ang-X),P[1]+S*dsin(Ang-X)]; /* 矢先 */ - V=(X==0)?[U,V]:[U,P,V]; - if(getopt(ar)==1) V=append([Q,P,0],V); /* 心棒 */ - return xylines(V|option_list=Opt); - }else if(AR>4&&AR<9){ - Ang=myarg([Q[0]-P[0],Q[1]-P[1]]); - ANG=[0.7854,0.5236,0.3927,0.2618]; - X=ANG[AR-5]; - U=[P[0]+S*dcos(Ang+X),P[1]+S*dsin(Ang+X)]; - V=[P[0]+S*dcos(Ang-X),P[1]+S*dsin(Ang-X)]; - W=ptcommon([P,U],[P,Q]|in=-2); - W1=[(U[0]+P[0]+W[0])/3,(U[1]+P[1]+W[1])/3]; - W2=[(V[0]+P[0]+W[0])/3,(V[1]+P[1]+W[1])/3]; - L=[U,W1,P,1,W2,V]; - if(getopt(ar)==1) L=append([Q,P,0],L); - if(type(Sc=getopt(scale))>0){ - if(type(Sc)==1) Sc=[Sc,Sc]; - L=ptaffine(diagm(2,Sc),L); + if(isint(S)&&S<0&&S>-8){ + if((S=-S)==6||S==7){ + H=ptcommon([Q,R],[P,0]); + if(S==6) return xyang(H,P,0,0|option_list=getopt()); /* 円 */ + return xylines([P,H]|option_list=getopt()); /* 垂線 */ } - Opt=getopt(opt); - if(type(Opt)>0) OL=[["opt",Opt]]; - else OL=[]; - if(getopt(proc)==1) return append([2,OL],L); - S=xybezier(L|optilon_list=OL); - if(getopt(dviout)!=1) return S; - dviout(S); - return 1; + O=pt5center(P,Q,R); + if(S==2) H=P; /* 外心 */ + else{ + if(S>2) S++; /* 内心，傍心 */ + H=ptcommon([P,Q],[O[S],0]); + } + return xyang(H,O[S],0,0|option_list=getopt()); } + if(type(Ar=getopt(ar))!=1) Ar=0; + if(isint(R)){ + if(R==1||R==-1){ /* 直角 */ + P1=ptcommon([Q,P],[-S,0]); + S*=R; + P2=ptcommon([P,P1],[S,@pi/2]); + P3=ptcommon([P1,P2],[S,@pi/2]); + return xylines([P1,P2,P3]|option_list=Opt); + }else if((AR=abs(R))==0||AR==2||AR==3||AR==4||AR>=10){ /* 矢印 */ + Ang=myarg([Q[0]-P[0],Q[1]-P[1]]); + if(R<0) Ang+=3.14159; + if(AR>10) X=deval(@pi/180*AR); + else{ + ANG=[0.7854,0.5236,1.0472]; + X=(AR==0)?1.5708:ANG[AR-2]; + } + U=[P[0]+S*dcos(Ang+X),P[1]+S*dsin(Ang+X)]; + V=[P[0]+S*dcos(Ang-X),P[1]+S*dsin(Ang-X)]; /* 矢先 */ + L=(X==0)?[U,V]:[U,P,V]; + if(X&&iand(Ar,2)){ + L=append([V],L); + if((X=ptcommon([P,Q],[U,V]|in=1))!=0) P=X; + } + if(iand(Ar,1)) + L=append([Q,P,0],L); /* 心棒 */ + return xylines(L|option_list=Opt); + }else if(AR>4&&AR<9){ + Ang=myarg([Q[0]-P[0],Q[1]-P[1]]); + ANG=[0.7854,0.5236,0.3927,0.2618]; + X=ANG[AR-5]; + U=[P[0]+S*dcos(Ang+X),P[1]+S*dsin(Ang+X)]; + V=[P[0]+S*dcos(Ang-X),P[1]+S*dsin(Ang-X)]; + W=ptcommon([P,U],[P,Q]|in=-2); + W1=[(U[0]+P[0]+W[0])/3,(U[1]+P[1]+W[1])/3]; + W2=[(V[0]+P[0]+W[0])/3,(V[1]+P[1]+W[1])/3]; + L=iand(Ar,2)?[V,U,1,W1,P,1,W2,V]:[U,W1,P,1,W2,V]; + if(iand(Ar,1)){ + if(iand(Ar,2)) P=ptcommon([P,Q],[U,V]); + L=append([Q,P,0],L); + }; + if(type(Sc=getopt(scale))>0){ + if(type(Sc)==1) Sc=[Sc,Sc]; + L=ptaffine(diagm(2,Sc),L); + } + Opt=delopt(Opt,"proc"); + if(getopt(proc)==1) return append([2,Opt],L); + S=xybezier(L|option_list=Opt); + if(getopt(dviout)!=1) return S; + dviout(xyproc(S)); + return 1; + } + } } if(type(Q)<3){ X=deval(Q); Y=deval(R); @@ -15431,6 +18588,178 @@ def xypoch(W,H,R1,R2) return S; } +def xycircuit(P,S) +{ + if(type(Sc=getopt(scale))!=1) Sc=1; + if(type(Opt0=getopt(opt))!=7) Opt0=""; + if(type(At=getopt(at))!=1) At=(S=="E"||S=="EE")?1:1/2; + Rev=(getopt(rev)==1)?-1:1; + if(type(P)==4&&type(car(P))==4&&P[0][0]==P[1][0]) Rev=-Rev; + W=R=B2=B3=0;Opt=Opt2=Opt3=""; + if(S=="L"||S=="VL"||S=="LT"){ + G=[1/8*x-2/5*cos(x)+2/5,1/2*sin(x)+1/2]; + B=xygraph(G,-21,[0,7*@pi],[-1,10],[-2,2]|scale=0.3/1.06466,opt=0); + B=append(B,[1,[1,0]]); + B=append([[0,0],car(B),1],cdr(B)); + W=1;Opt="thick"; + if(S=="VL"){ + B2=xyang(0.2,[0.5+0.4*Rev,0.45],[0.5-0.435*Rev,-0.3],3|ar=3,opt=0); + Opt2="thick,fill"; + }else if(S=="LT"){ + B2=[[0.5+0.4*Rev,0.45],[0.5-0.435*Rev,-0.3],0,[0.45+0.4*Rev,0.394],[0.55+0.4*Rev,0.506]]; + Opt2="thick"; + } + }else if(S=="C"||S=="VC"||S=="C+"||S=="C-"||S=="CT"){ + B=[[0,-0.2],[0,0.2],0,[0.15,-0.2],[0.15,0.2]]; + W=0.15;Opt="very thick"; + if(S=="VC"){ + B2=xyang(0.2,[1/3+0.075,0.3*Rev],[-1/3+0.075,-0.3*Rev],3|ar=3,opt=0); + Opt2="thick,fill"; + }else if(S=="CT"){ + B2=[[1/3+0.075,0.3*Rev],[-1/3+0.075,-0.3*Rev],0,[1/3+0.125,0.244*Rev], + [1/3+0.025,0.356*Rev]]; + Opt2="thick"; + }else if(S=="C+") + B2=[[0,0.05],[0.15,-0.05],0,[0,0.15],[0.15,0.05],0,[0,-0.05],[0.15,-0.15], + 0,[0.29,0.04*Rev],[0.29,0.24*Rev],0,[0.19,0.14*Rev],[0.39,0.14*Rev]]; + else if(S=="C-") + B2=[[0,0.05],[0.15,-0.05],0,[0,0.15],[0.15,0.05],0,[0,-0.05],[0.15,-0.15]]; + }else if(S=="R"||S=="VR"||S=="VR3"||S=="RT"){ + for(I=0,B=[[0,0]];I<12;I++) + if(iand(I,1)) B=cons([I,(-1)^((I+1)/2)],B); + B=reverse(cons([12,0],B)); + B=xylines(B|scale=[1/18,0.15],opt=0); + W=2/3;Opt="thick"; + if(S=="VR"){ + B2=xyang(0.2,[2/3,0.3*Rev],[0,-0.3*Rev],3|ar=3,opt=0); + Opt2="thick,fill"; + }else if(S=="RT"){ + B2=[[2/3,0.3*Rev],[0,-0.3*Rev],0,[0.717,0.244*Rev],[0.617,0.357*Rev]]; + Opt2="thick"; + }else if(S=="RN3"){ + B2=xyang(0.2,[1/3,0.2*Rev],[1/3,0.5*Rev],3|ar=3,opt=0); + Opt2="thick,fill"; + } + }else if(S=="RN"||S=="VRN"||S=="RN3"||S=="NRT"){ + B=xylines([[0,0.1],[2/3,0.1],[2/3,-0.1],[0,-0.1],[0,0.1]]|opt=0); + W=2/3;Opt="thick"; + if(S=="VRN"){ + B2=xyang(0.2,[2/3,0.3*Rev],[0,-0.3*Rev],3|ar=3,opt=0); + Opt2="thick,fill"; + }else if(S=="RN3"){ + B2=xyang(0.2,[1/3,0.2*Rev],[1/3,0.5*Rev],3|ar=3,opt=0); + Opt2="thick,fill"; + }else if(S=="NRT"){ + B2=[[2/3,0.3*Rev],[0,-0.3*Rev],0,[0.717,0.244*Rev],[0.617,0.357*Rev]]; + Opt2="thick"; + } + }else if(S=="circle"){ + W=1; + B=xyang(0.5,[0.5,0],0,0|opt=0); + }else if(S=="gap"){ + W=0.3; + B=xyang(0.15,[0.15,0],0,3.1416|opt=0); + }else if(S=="E"){ + W=0.1; + B=[[0,0.2],[0,-0.2],0,[0,0.05],[0.1,-0.05],0,[0,0.15],[0.1,0.05],0,[0,-0.05],[0.1,-0.15]]; + }else if(S=="EE"){ + W=0.15; + B=[[0,0.2],[0,-0.2],0,[0.075,0.13],[0.075,-0.13],0,[0.15,-0.06],[0.15,0.06]]; + }else if(S=="Cell"){ + W=0.1; + B=[[0,-0.2],[0,0.2]]; + B2=[[0.1,-0.1],[0.1,0.1]];Opt2="very thick"; + }else if(S=="Cell2"){ + W=0.3; + B=[[0,-0.2],[0,0.2],0,[0.2,-0.2],[0.2,0.2]]; + B2=[[0.1,-0.1],[0.1,0.1],0,[0.3,-0.1],[0.3,0.1]];Opt2="very thick"; + }else if(S=="Cells"){ + W=0.6; + B=[[0,-0.2],[0,0.2],0,[0.5,-0.2],[0.5,0.2],0,[0.1,0],[0.18,0],0, + [0.24,0],[0.34,0],0,[0.40,0],[0.5,0]]; + B2=[[0.1,-0.1],[0.1,0.1],0,[0.6,-0.1],[0.6,0.1]];Opt2="very thick"; + }else if (S=="Sw"){ + W=0.5; + B=xyang(0.05,[0.05,0],0,0|opt=0); + B0=ptaffine(1,B|shift=[0.4,0]); + B=ptaffine("union",[B,B0]); + B=ptaffine("union",[B,[[0.0908,0.025*Rev],[0.45,0.17*Rev]]]); + }else if(S=="D"){ + W=0.3;Opt="thick"; + B=[[0,0],[0.3,0.173],0,[0.3,0.173],[0.3,-0.173],0,[0.3,-0.173],[0,0],0, + [0,0.173],[0,-0.173]]; + }else if(S=="NPN"||S=="PNP"||S=="NPN0"||S=="PNP0"){ + W=0.6; + C=[[0.6,0],[0.37,0.23],[0,0],[0.23,0.23]]; + if(Rev==-1) C=[C[2],C[3],C[0],C[1]]; + if(S=="PNP"||S=="PNP0") C=[C[1],C[0],C[2],C[3]]; + B=[[0,0],[0.23,0.23],0,[0.6,0],[0.37,0.23],0,[0.3,0.23],[0.3,0.6]]; + B=ptaffine("union",[xyang(0.15,C[0],C[1],18|ar=1,opt=0),B]); + if(S=="PNP"||S=="NPN") B=ptaffine("union",[xyang(0.3354,[0.3,0.15],0,0|opt=0),B]); + B2=[[0.07,0.23],[0.53,0.23]]; + Opt2="very thick"; + }else if(S=="JN"||S=="JP"){ + W=0.6; + B=[[0,0],[0.2,0],1,[0.2,0.23],0,[0.6,0],[0.4,0],1,[0.4,0.23],0,[0.3,0.23],[0.3,0.6]]; + C=[[0.3,0.23],[0.3,0.4854]]; + if(S=="JP") C=reverse(C); + B=ptaffine("union",[B,xyang(0.15,C[0],C[1],18|opt=0)]); + B=ptaffine("union",[B,xyang(0.3354,[0.3,0.15],0,0|opt=0)]); + B2=[[0.07,0.23],[0.53,0.23]]; + Opt2="very thick"; + }else if(S=="") R=(Opt0=="")?xyline(P[0],P[1]):xyline(P[0],P[1]|opt=Opt0); + else if(S=="arrow") R=xyang(0.2*Sc,P[1],P[0],3|ar=1,opt=Opt0); + else if(type(S)==4&&type(car(S))==7){ + if(type(car(P))!=4) P=[P]; + for(R="";P!=[];P=cdr(P)) R+=xyput([car(P)[0],car(P)[1],car(S)]); + } + if(W){ + R=""; + if(type(P)==4){ + if(type(car(P))==4){ + T=ptcommon([[0,0],[1,0]],P|in=2); + L=dnorm(P); + W*=Sc; + L1=L*At-W/2;L2=L*(1-At)-W/2; + if(L1>0){ + P1=[P[0][0]+L1*dcos(T),P[0][1]+L1*dsin(T)]; + R+=xyline(P[0],P1); + } + if(L2>0){ + P2=[P[1][0]-L2*dcos(T),P[1][1]-L2*dsin(T)]; + R+=xyline(P2,P[1]); + } + B=ptaffine(Sc,B|shift=P1,arg=T); + if(B2) B2=ptaffine(Sc,B2|shift=P1,arg=T); + if(B3) B3=ptaffine(Sc,B3|shift=P1,arg=T); + }else{ + B=ptaffine(Sc,B|shift=P1); + if(B2) B2=ptaffine(Sc,B2|shift=P1); + if(B3) B3=ptaffine(Sc,B3|shift=P1); + } + }else{ + B=ptaffine(Sc,B); + if(B2) B2=ptaffine(Sc,B2); + if(B3) B3=ptaffine(Sc,B3); + } + if(Opt=="") Opt=Opt0; + else if(Opt0!="") Opt=Opt+","+Opt0; + R+=(Opt=="")?xybezier(B):xybezier(B|opt=Opt); + if(B2){ + if(Opt2=="") Opt2=Opt0; + else if(Opt0!="") Opt2=Opt2+","+Opt0; + R+=(Opt2=="")?xybezier(B2):xybezier(B2|opt=Opt2); + } + if(B3){ + if(Opt3=="") Opt3=Opt0; + else if(Opt0!="") Opt3=Opt3+","+Opt0; + R+=(Opt3=="")?xybezier(B3):xybezier(B3|opt=Opt3); + } + } + return R; +} + + def ptaffine(M,L) { if(type(L)!=4&&type(L)!=5){ @@ -15573,8 +18902,16 @@ def ptlattice(M,N,X,Y) for(L=[],I=M-1;I>=0;I--){ for(P0=P1=0,J=N-1;J>=0;J--){ P=Org+I*X+J*Y; - for(C=Cond; C!=[]; C=cdr(C)) - if(subst(car(C),x,P[0],y,P[1])<0) break; + for(C=Cond; C!=[]; C=cdr(C)){ + TC=car(C); + if(type(TC)<4) + if(subst(TC,x,P[0],y,P[1])<0) break; + else{ + for(;TC!=[];TC=cdr(TC)) + if(subst(car(TC),x,P[0],y,P[1])>=0) break; + if(TC==[]) break; + } + } if(C!=[]) continue; if(Line) F[I][J]=1; else L=cons(vtol(S*P),L); @@ -15647,6 +18984,75 @@ def ptwindow(L,X,Y) return reverse(R); } +def pt5center(P,Q,R) +{ +/* P=[1,[0,0]];Q=[[0,0],[1,0]];R=[[0,0],[0,1]]; */ + if(length(P)==2&&type(P[0])==4){ /* circle */ + if(type(Q)==4&&type(Q[1])==4){ /* line */ + A=myarg(lsub(Q));B=myarg(lsub(R));X0=ptcommon(Q,R); + M=mrot(-A);N=mrot(A);X=M*ltov(X0);O=M*ltov(P[0]); + if(!(L=B-A)) return 0; + Pi=deval(@pi);for(;L<0;L+=Pi);for(;L>Pi;L-=Pi); + XX=X[0]+y*deval(cos(L/2))/deval(sin(L/2)); + XY=X[1]+y; + if(getopt(neg)==1){ + XX=subst(XX,y,-y);XY=subst(XY,y,-y); + } +/* mycat([[P[0],O],XX,XY]); */ + V=(XX-O[0])^2+(XY-O[1])^2; +/* mycat(V-(y+P[0])^2); */ + S=polroots(V-(y+P[1])^2,y); + S=append(polroots(V-(y-P[1])^2,y),S); + S=qsort(S);V=ltov([XX,XY]); +/* mycat([S,V,M,N,M*N]); */ + for(R0=[],ST=S;ST!=[];ST=cdr(ST)) R0=cons([vtol(N*subst(V,y,car(ST))), car(ST)],R0); +/* mycat(R0); */ + for(R=[],F=1;R0!=[];R0=cdr(R0)){ + if(car(R0)[1]>=0) R=cons(car(R0),R); + else{ + if(F){ + F=0; R=reverse(R); + } + R=cons(car(R0),R); + } + } +/* mycat(R); */ + if(!F) R=reverse(R); + return R; + } + } + L=newvect(7); + L[2]=ptcommon([P,Q],[P,R]|in=-1); + Q1=ptcommon([P,R],[Q,0]);R1=ptcommon([P,Q],[R,0]); + L[3]=ptcommon([Q,Q1],[R,R1]); + P=ltov(P);Q=ltov(Q);R=ltov(R); + A=dnorm([Q,R]);B=dnorm([P,R]);C=dnorm([P,Q]); + L[0]=vtol((P+Q+R)/3); + L[1]=vtol((A*P+B*Q+C*R)/(A+B+C)); + L[4]=vtol((-A*P+B*Q+C*R)/(-A+B+C)); + L[5]=vtol((A*P-B*Q+C*R)/(A-B+C)); + L[6]=vtol((A*P+B*Q-C*R)/(A+B-C)); + if((V=getopt(opt))==0){ + H1=ptcommon([Q,R],[1,1]|in=1); + H2=ptcommon([R,P],[1,1]|in=1); + H3=ptcommon([P,Q],[1,1]|in=1); + return [L(0),H1,H2,H3]; + }else if(V==1||V==4||V==5||V==6){ + H1=ptcommon([Q,R],[L[1],0]); + H2=ptcommon([R,P],[L[1],0]); + H3=ptcommon([P,Q],[L[1],0]); + return [[L[1],dnorm(L[1],H1)],H1,H2,H3]; + }else if(V==2){ + return [L[2],dnorm([L[2],P])]; + }else if(V==3){ + H1=ptcommon([Q,R],[P,0]); + H2=ptcommon([R,P],[Q,0]); + H3=ptcommon([P,Q],[R,0]); + return [L[3],H1,H2,H3]; + } + return vtol(L); +} + def lninbox(L,W) { if(L[0]==L[1]) return 0; @@ -15717,6 +19123,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"){ @@ -16422,16 +19849,16 @@ def getbygrs(M, TT) if(MT[S] > MT[0]) S = 0; } } - M = reverse(R); + M = reverse(R); R = getopt(var); if(type(R)<1){ for(R = [], I = J-1; I >= 0; I--) R = cons(asciitostr([97+I]), R); - } + } Sft=(Sft>=0)?1:0; if(General < 0) - Sft=-Sft-1; - M = sp2grs(M,R,Sft|mat=Mat); + Sft=-Sft-1; + M = sp2grs(M,R,Sft|mat=Mat); } for(M0=[],MM=M;MM!=[];MM=cdr(MM)){ /* change "?" -> z_z */ for(M1=[],Mm=car(MM);Mm!=[];Mm=cdr(Mm)){ @@ -17156,19 +20583,19 @@ def getbygrs(M, TT) } } }else /* Rigid case */ - P = dx^(AL[0][1][0][0][0]); + P = dx^(AL[0][1][0][0][0]); /* additions and middle convolutions */ for(ALT = AL; ALT != []; ALT = cdr(ALT)){ AI = car(ALT)[0]; if(type(AI) != 4) continue; V = ltov(AI); - if(V[0] != 0) P = mc(P,x,V[0]); + if(V[0] != 0) P = mc(P,x,V[0]); for(I = 1; I < NP; I++){ if(V[I] != 0) - P = sftexp(P,x,L[I],-V[I]); - } + P = sftexp(P,x,L[I],-V[I]); + } } - P = (Simp>=0)? simplify(P,Fuc,4|var=[dx]):simplify(P,Fuc,4); + P = (Simp>=0)? simplify(P,Fuc,4|var=[dx]):simplify(P,Fuc,4); if(TeX >= 0){ Val = 1; if(mydeg(P,dx) > 2 && AMSTeX == 1 && TeXEq > 3) @@ -17387,6 +20814,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); @@ -17476,13 +20936,59 @@ def conf1sp(M) return P; } -/* ((1)(1)) ((1)) 111|11|21 [[ [2,[ [1,[1]],[1,[1]] ]], [1,[[1,[1]]]] ]] */ +def s2cspb(L) +{ + Sub=(getopt(sub)); + if(Sub!=1&&Sub!=2&&Sub!=-1) Sub=0; + if(type(L)==7){ + if(Sub>0){ + I=str_char(L,0,"("); + if(I<0) return car(s2sp(L)); + for(R=[];;){ + J=str_char(L,I,"("); + if(J<0) return reverse(R); + K=str_pair(L,J+1,"(",")"); + if(K<0) return reverse(R); + R=cons(s2cspb(str_cut(L,J+1,K-1)|sub=1),R); + I=K; + } + }else{ + K=str_len(L); + for(R=[],I=0;I0){ + if(N==-1) S=s2cspb(S); for(D=0,S=reverse(S);S!=[];S=cdr(S),D++){ if(D) U=","+U; T=str_subst(rtostr(car(S)),","," "); @@ -17522,27 +21028,15 @@ def s2csp(S) } return U; } - S=strtoascii(S); if(type(N=getopt(n))>0){ - S=ltov(S); - L=length(S); - R=""; - for(I=J=N=0, V=[];J47&&S[J]<58) N=N*10+S[J]-48; - else{ - if(N>0){ - V=cons(N,V); - N=0; - } - if(S[J]==41){ /* ) */ - - }else if(S[J]==44){ /* , */ - - } - } + if(N==-1) return s2cspb(S|sub=-1); + else{ + S=s2cspb(S); + if(N==1) S=s2csp(S); + return S; } } + S=strtoascii(S); for(P=TS=[],I=D=0; S!=[]; S=cdr(S)){ if((C=car(S))==44){ /* , */ P=cons(D,P);D=0; @@ -17574,18 +21068,56 @@ def s2csp(S) return reverse(R); } +/* +def confspt(S) +{ + if(!isint(F=getopt(sub))) F=0; + N=length(S); + P=newmat(N,N); + for(I=0;I=0) TP=cons([I],TP); + for(F=1;F;){ + for(T=TP,F=0,S=length(car(TP));T!=[];T=cdr(T)){ + if(length(T0=car(T))length(T)) return []; if(type(Op=getopt(opt))!=1) Op=0; - else{ - VS=ltov(S); - L=length(S)-1; - VT=ltov(qsort(T)); - } + VS=ltov(S); + L=length(S)-1; + VT=ltov(qsort(T)); if(length(S)==length(T)){ - if(S==T||qsort(S)==qsort(T)) R=S; + if((R=S)==T|| (R=qsort(S))==qsort(T)){ + for(S=[];R!=[];R=cdr(R)) S=cons([car(R),[car(R)]],S); + return S; + } else return []; }else if(getopt(sort)==1){ S0=S1=[]; @@ -17680,12 +21212,124 @@ def confspt(S,T) } #endif +def vConv(K,I,J) +{ + if(type(X=getopt(var))!=7) X="a"; + if(getopt(e)==2) return subst(vConv(K,I+1,J+1),makev([X,1]),0); + if(J>K){L=J;J=K;K=L;} + if(K>I||J<1||K+J0;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;I0){ FE=cons([C*(X+S)^2-C1^(1/2)*(X+S)+C0,car(FT)[1]],FE); @@ -18193,7 +21850,7 @@ def pfrac(F,X) Q += P/(FC[I][0]^K); Q = red(Q); } - } + } L=reverse(L); Q = nm(red(red(Q*FD)-FN)); Q = ptol(Q,X); @@ -18210,33 +21867,31 @@ def pfrac(F,X) for(;RT!=[];RT=cdr(RT)){ RTT=car(RT); R=mtransbys(os_md.substblock,R,[RTT^(1/2),(RTT^(1/2))^2,RTT]); - } - TeX=getopt(TeX); - if((Dvi=getopt(dviout))==1||TeX==1){ - V=strtov("0"); - for(S=L=0,RR=R;RR!=[];RR=cdr(RR),L++){ - RT=car(RR); + } + if(Dvi==1||TeX==1){ + V=strtov("0"); + for(S=L=0,RR=R;RR!=[];RR=cdr(RR),L++){ + RT=car(RR); S+=(RT[0]/RT[1]^RT[2])*V^L; - } + } if(TeX!=1) fctrtos(S|var=[V,""],dviout=1); - else return fctrtos(S|var=[V,""],TeX=3); + else return fctrtos(S|var=[V,""],TeX=3,lim=0); } return reverse(R); } def cfrac(X,N) { - F=[floor(X)]; - if(N<0){ - Max=N=-N; - } - X-=F[0]; - if(Max!=1) - M=mat([F[0],1],[1,0]); + if(!ntype(X)&&N==0) N=2*dn(X)+1; + if(N<0) Max=N=-N; + Ng=(getopt(neg)==1)?1:0; + F=[Ng?ceil(X):floor(X)]; + X=Ng?F[0]-X:X-F[0]; + if(Max!=1) M=mat([F[0],1],[1,0]); for(;N>0 && X!=0;N--){ X=1/X; - F=cons(Y=floor(X),F); - X-=Y; + F=cons(Y=Ng?ceil(X):floor(X),F); + X=Ng?Y-X:X-Y; if(Max){ M0=M[0][0];M1=M[1][0]; M=M*mat([Y,1],[1,0]); @@ -18381,12 +22036,79 @@ def sqrt2rat(X) def cfrac2n(X) { + if(((Q=getopt(q))==1||Q==-1)&&isall(os_md.isint,X)){ /* q-def */ + F=car(X); + R=[red((1-q^F)/(1-q))];X=cdr(X); + for(SQ=1/q;X!=[];SQ=1/SQ,X=cdr(X)){ + G=car(X); + V=(Q==1)?[(1/SQ)^F,(1-SQ^G)/(1-SQ)]:[-q^(F-1),(1-q^G)/(1-q)]; + R=cons(red(V),R); + F=G; + } + return cfrac2n(reverse(R)|ex=1); + } + if((Ex=getopt(ex))!=1) Ex=0; +/* + if(!isint(Fn=getopt(dn))) Dn=0; + else{ + RD=[];V=0; + } +*/ + if(isvar(car(X))&&length(X)==4&&isint(X[1])){ /* [x,n,g(x),f(x)] */ + A=newmat(2,2);B=mgen(2,"diag",[1],0); + for(I=1;I<=X[1];I++){ + A[1][1]=0;A[0][1]=1; + A[1][0]=Ex?myfeval(X[2],[X[0],I]):subst(X[2],X[0],I); + A[0][0]=Ex?myfeval(X[3],[X[0],I]):subst(X[3],X[0],I); + if(vars(A)!=[]) A=red(A); + B=B*A; + if(vars(B)!=[]) B=red(B); + } + if(getopt(var)==1) return [B[1][0],B[0][0]]; + B=B[1][0]/B[0][0]; + if(vars(B)!=[]) B=red(B); + return B; + } + if(Ex||(type(car(X))==4&&length(car(X))==2)){ + if(type(car(X))!=4){ + N=car(X);X=cdr(X); + } + if(getopt(reg)==1){ /* normal */ + for(R=[N], F=1;X!=[];X=cdr(X)) { + if(type(car(X))==4){ + F=car(X)[0]/F; + R=cons(car(X)[1]/F,R); + }else{ + F=1/F; + R=cons(car(X)/F,R); + } + } + return reverse(R); + } + A=newmat(2,2);B=mgen(2,"diag",[1],0); + for(I=0,TX=X;TX!=[];TX=cdr(TX),I++){ + A[1][1]=0;A[0][1]=1; + if(type(car(TX))!=4){ + A[1][0]=1;A[0][0]=car(TX); + }else{ + A[1][0]=car(TX)[0];A[0][0]=car(TX)[1]; + } + if(vars(A)!=[]) A=red(A); + B=B*A; + if(vars(B)!=[]) B=red(B); + } + if(getopt(var)==1) return [N,B[1][0],B[0][0]]; + B=N+B[1][0]/B[0][0]; + if(vars(B)!=[]) B=red(B); + return B; + } if(type(L=getopt(loop))==1&&L>0) C=x; else{ C=0;L=0; } - if(L>1){ + Sg=getopt(neg)==1?-1:1; + if(L>1){ /* circulate */ for(Y=[];L>1;L--){ Y=cons(car(X),Y); X=cdr(X); @@ -18394,21 +22116,21 @@ def cfrac2n(X) if(X!=[]){ P=cfrac2n(X|loop=1); for(V=P,Y=reverse(Y);Y!=[];Y=cdr(Y)) - V=sqrt2rat(car(Y)+1/V); + V=sqrt2rat(car(Y)+Sg/V); return V; }else{ C=0;X=reverse(Y); } } for(V=C,X=reverse(X);X!=[];X=cdr(X)){ - if(V!=0) V=1/V; + if(V!=0) V=Sg/V; V+=car(X); } if(C!=0){ V=red(V);P=dn(V)*x-nm(V); S=getroot(P,x|cpx=2); T=map(eval,S); - V=(T[0]>0)?S[0]:S[1]; + V=T[0]>0?S[0]:S[1]; } return V; } @@ -18588,7 +22310,7 @@ def sp2grs(M,A,L) else{ if(LM == 1) MN = cons([V, (J==0)?0:makev([AA])], MN); - else if(I == 1 && Mat == 0) + else if(I == 1 && getopt(con)==1 && Mat == 0) MN = cons([V, (J==length(MI)-1)?0:makev([AA,J+Sft])], MN); else MN = cons([V, (J==0)?0:makev([AA,J])], MN); @@ -19818,129 +23540,120 @@ def bernoulli(N) } /* linfrac01([x,y]) */ -/* linfrac01(newvect(10,[0,1,2,3,4,5,6,7,8,9]) */ -/* 0:x=0, 1:x=y, 2:x=1, 3:y=0, 4:y=1, 5:x=\infty, 6:y=\infty, 7:x=y=0, 8:x=y=1, 9:x=y=\infty - 10:y_2=0, 11:y_2=x, 12:y_2=y, 13: y_2=1, 14: y_2=\infty - 15:y_3=0, 16:y_3=x, 17:y_3=y, 18: y_3=y_2, 19: y_3=1, 20:y_3=\infty - X[0],X[11],X[2],X[10],X[13],X[5],X[14],X[7],X[8],X[9], - X[3],X[1],X[12],X[4],X[6] +/* (x_0,x_1,x_2,x_3,...,x_{q+3})=(x,0,1,y_1,...,y_q,\infty) T=0 (x_2,x_1,x_3,x_4,...) T=-j (x_1,x_2,..,x_{j-1},x_{j+1},x_j,x_{j+2},...) T=1 (1-x_1,1-x_2,1-x_3,1-x_4,...) T=2 (1/x_1,1/x_2,1/x_3,1/x_4,...) - T=3 (x_1,x_1/x_2,x_1/x_3,x_1/x_4,...) + T=3 (1/x_1,x_2/x_1,x_3/x_1,x_4/x_1,...) + ... */ - def lft01(X,T) -{ - MX=getopt(); +{ + S=0; if(type(X)==4){ + if(type(car(X))==4){ + S=X[1];X=car(X); + } K=length(X); if(K>=1) D=1; } - if(type(X)==5){ - K=length(X); - for(J=5, F=K-10; F>0; F-=J++); - if(F==0) D=2; - } - if(D==0) return 0; - if(T==0){ /* x <-> y */ - if(D==1){ - R=cdr(X); R=cdr(R); - R=cons(X[0],R); - return cons(X[1],R); + if(D==0) return 0; + if(type(T)==4&&type(car(T))==4&&length(T)==2){ + U=newvect(K+3); + for(I=0;I=0;I--) U=cons(I,U); + if(length(T)==2) T=mperm(U,[T],0); + L=sexps(T); + for(R=[X,S];L!=[];L=cdr(L)){ + if(!(I=car(L))) I=4; + /* else if(I==1) I=1; */ + else if(I==2) I=5; + else if(I==K+1) I=6; + else if(I>2) I=2-I; + R=lft01(R,I); } - R=newvect(K,[X[3],X[1],X[4],X[0],X[2],X[6],X[5]]); - for(I=7;I=0;I--) S=cons(I,S); + else S=0; + if(T<=0){ /* y_i <-> y_{i+1}, y_0=x=x_0, y_i=x_{i+2} */ + R=mperm(X,[[-T,1-T]],0); + if(S){ + if(!T) S=mperm(S,[[0,3]],0); + else S=mperm(S,[[2-T,3-T]],0); /* : J J=3,...,K; */ + R=[R,S]; } - R=newvect(K,[X[2],X[1],X[0],X[4],X[3],X[5],X[6],X[8],X[7],X[9]]); - for(I=11;I1;X=cdr(X)) R=cons(red(car(X)*(1-T)/(car(X)-T)),R); + R=cons(1-T,R); + if(S) S=mperm(S,[[K+1,K+2]],0); + }else if(T==7){ /* x_2=1 <-> x_{K+2}=infty */ + for(R=[];X!=[];X=cdr(X)) R=cons(red(car(X)/(car(X)-1)),R); + if(S) S=mperm(S,[[2,K+2]],0); + }else if(T==8) { /* x_0=x <-> x_{K+2}=infty */ + T=car(X); + for(R=[1-car(X)], X=cdr(X); X!=[]; X=cdr(X)) R=cons(red((T-1)*car(X)/(T-car(X))),R); + if(S) S=mperm(S,[[0,K+2]],0); + }else return 0; + R=reverse(R); + return S?[R,S]:R; } def linfrac01(X) { - if(type(X)==4) K=length(X)-2; - else if(type(X)==5){ - L=length(X); - for(K=0,I=10,J=5; I=0;I--) U=cons(I,U); + X=[car(X),U]; + }else U=0; } if(K>3 && getopt(over)!=1) return(-1); II=(K==-1)?3:4; @@ -19955,7 +23668,7 @@ def linfrac01(X) } } } - return L; + return reverse(L); }