=================================================================== RCS file: /home/cvs/OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v retrieving revision 1.35 retrieving revision 1.51 diff -u -p -r1.35 -r1.51 --- OpenXM/src/asir-contrib/packages/src/os_muldif.rr 2018/09/26 04:49:44 1.35 +++ OpenXM/src/asir-contrib/packages/src/os_muldif.rr 2019/06/28 01:54:31 1.51 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.34 2018/09/25 00:13:52 takayama Exp $ */ +/* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.50 2019/06/27 02:53:26 takayama Exp $ */ /* The latest version will be at ftp://akagi.ms.u-tokyo.ac.jp/pub/math/muldif scp os_muldif.[dp]* ${USER}@lemon.math.kobe-u.ac.jp:/home/web/OpenXM/Current/doc/other-docs */ @@ -6,7 +6,7 @@ /* #undef USEMODULE */ /* os_muldif.rr (Library for Risa/Asir) - * Toshio Oshima (Nov. 2007 - Sep. 2018) + * Toshio Oshima (Nov. 2007 - Feb. 2019) * * For polynomials and differential operators with coefficients * in rational funtions (See os_muldif.pdf) @@ -43,6 +43,7 @@ static Canvas$ static ID_PLOT$ static Rand$ static LQS$ +static SVORG$ localf spType2$ localf erno$ localf chkfun$ @@ -116,6 +117,7 @@ localf lchange$ localf llsize$ localf llbase$ localf lsort$ +localf rsort$ localf lpair$ localf lmax$ localf lmin$ @@ -242,6 +244,8 @@ localf seriesHG$ localf seriesMc$ localf seriesTaylor$ localf mulpolyMod$ +localf solveEq$ +localf baseODE$ localf taylorODE$ localf evalred$ localf toeul$ @@ -302,6 +306,7 @@ localf cutgrs$ localf mcgrs$ localf mc2grs$ localf mcmgrs$ +localf spslm$ localf anal2sp$ localf delopt$ localf str_char$ @@ -349,6 +354,10 @@ localf mcop$ localf shiftop$ localf conf1sp$ localf confexp$ +localf confspt$ +localf mcvm$ +localf s2csp$ +localf partspt$ localf pgen$ localf diagm$ localf mgen$ @@ -461,11 +470,12 @@ extern Canvas$ extern ID_PLOT$ extern Rand$ extern LQS$ +extern SV=SVORG$ #endif static S_Fc,S_Dc,S_Ic,S_Ec,S_EC,S_Lc$ static S_FDot$ extern AMSTeX$ -Muldif.rr="00180925"$ +Muldif.rr="00190620"$ AMSTeX=1$ TeXEq=5$ TeXLim=80$ @@ -484,6 +494,7 @@ LCOPT=["red","green","blue","yellow","cyan","magenta", COLOPT=[0xff,0xff00,0xff0000,0xffff,0xffff00,0xff00ff,0,0xffffff,0xc0c0c0]$ LPOPT=["above","below","left","right"]$ LFOPT=["very thin","thin","dotted","dashed"]$ +SVORG=["x","y","z","w","u","v","p","q","r","s"]$ Canvas=[400,400]$ LQS=[[1,0]]$ @@ -623,6 +634,7 @@ def mycat(L) Do = 1; } if(CR) print(""); + else print("",2); } def fcat(S,X) @@ -661,6 +673,7 @@ def mycat0(L,T) Do = 1; } if(T) print(""); + else print("",2); } def findin(M,L) @@ -888,6 +901,7 @@ def mulsubst(F,L) if(N == 0) return F; if(type(L[0])!=4) L=[L]; + if(getopt(lpair)==1||(type(L[0])==4&&length(L[0])>2)) L=lpair(L[0],L[1]); if(getopt(inv)==1){ for(R=[];L!=[];L=cdr(L)) R=cons([car(L)[1],car(L)[0]],R); L=reverse(R); @@ -1685,6 +1699,7 @@ def mtoupper(MM, F) if(type(St = getopt(step))!=1) St=0; Opt = getopt(opt); if(type(Opt)!=1) Opt=0; + if(type(Main=getopt(main))!=1) Main=0; TeX=getopt(dviout); if(type(Tab=getopt(tab))!=1 && Tab!=0) Tab=2; Line="\\text{line}"; @@ -1715,6 +1730,13 @@ def mtoupper(MM, F) Top+=(TeX)?"\\ ":" "; } PC=IF=1; + if(Opt>3){ + for(P=[1],K=0;K4 && length(Var=vars(nm(M[J][K])))==1){ J0=J;Jv=mydeg(nm(M[J0][K]),car(Var)); for(I=JJ;I1) continue; if(mydeg(MIK,T[0])0){ QF=1;Q0*=T; continue; } @@ -1855,10 +1881,10 @@ def mtoupper(MM, F) if(TeX){ Lout=cons(["\\hspace{",Tab*(St-3)-1,"mm}\\text{If }", Var,"=",T,","] ,Lout); - Lout=append(mtoupper(M0,F|step=St+1,opt=Opt,dviout=-2,tab=Tab),Lout); + Lout=append(mtoupper(M0,F|step=St+1,opt=Opt,dviout=-2,tab=Tab,main=Main),Lout); }else{ mycat([str_times(" ",St-1)+"If",Var,"=",T,","]); - mtoupper(M0,F|step=St+1,opt=Opt); + mtoupper(M0,F|step=St+1,opt=Opt,main=Main); } } } @@ -1875,11 +1901,13 @@ def mtoupper(MM, F) KRC=-red((T[2]*dn(M[J0][K]))/(T[1]*dn(M[I][K]))); for(II=K;II GRS */ G=s2sp(M|std=1); L=length(G); @@ -3332,6 +3364,16 @@ def llbase(VV,L) return V; } +def rsort(L,T,K) +{ + for(R=[];L!=[];L=cdr(L)) + R=cons((type(car(L))==4)?rsort(car(L),T-1,K):car(L),R); + if(T>0||iand(T,iand(K,2)/2)) return reverse(R); + R=qsort(R); + return (iand(K,1))? reverse(R):R; +} + + def lsort(L1,L2,T) { C1=getopt(c1);C2=getopt(c2); @@ -5415,6 +5457,12 @@ def mulpdo(P,Q,L); def transpdosub(P,LL,K) { + if(type(P)>3) 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; @@ -5440,12 +5488,11 @@ def transpdosub(P,LL,K) def transpdo(P,LL,K) { - if(type(K[0]) < 4) - K = [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]); @@ -5454,10 +5501,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); @@ -5510,7 +5572,8 @@ def texbegin(T,S) { if(type(Opt=getopt(opt))==7) Opt="["+Opt+"]\n"; else Opt="\n"; - return "\\begin{"+T+"}"+Opt+S+"%\n\\end{"+T+"}\n"; + U=(str_chr(S,str_len(S)-1,"\n")<0)?"%\n":""; + return "\\begin{"+T+"}"+Opt+S+U+"\\end{"+T+"}\n"; } def mygcd(P,Q,L) @@ -6165,6 +6228,203 @@ def mulpolyMod(P,Q,X,N) return R; } +def solveEq(L,V) +{ + Inv=0;K=length(V); + H=(getopt(h)==1)?1:0; + if(getopt(inv)==1){ + if(K!=length(L)) return -5; + Inv=1; + VN=makenewv(vars(L)|num=K); + for(TL=[],I=K-1;I>=0;I--) TL=cons(VN[I]-L[I],TL); + S=solveEq(TL,V|h=H); + if(type(S)!=4) return S; + return mysubst(S,[VN,V]|lpair=1); + } + for(TL=[];L!=[];L=cdr(L)) TL=cons(nm(red(car(L))),TL); + S=gr(TL,reverse(V),2); + if(length(S)!=K) return -1; + for(R=[],I=F=0;I3){ + Ord-=4; Hgr=1; + }else Hgr=0; + if(type(car(L))==4&&type(L[1])==7){ + Tt=L[1];L=car(L); + } + M=N=length(L); SV=SVORG; + if(type(Var=getopt(var))==4&&(In>0||length(Var)==N)){ + SV=Var; + M=length(SV); + if(type(car(SV))==2){ + for(R=[];SV!=[];SV=cdr(SV)) R=cons(rtostr(car(SV)),R); + SV=reverse(R); + } + }else{ + if(N>10){ + R=[]; + for(K=M-1;K>9;K++) R=cons(SV[floor(K/10)-1]+SV[K%10],R); + SV=append(SV,R); + } + 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=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 IL; + if(R==1){ + R=To;To=IL;IL=R; + } + L=mulsubst(L,[Var,IL]|lpair=1); + if(!In){ /* X_i'=\sum_j(\p_{x_j}X_i)*x_j' */ + for(TL=[],I=M-1;I>=0;I--){ + P=To[I];Q=mydiff(P,t); + for(J=0;J=0;I--){ + P=IL[I];Q=mydiff(P,t); + for(J=0;J=0;I--) L=cons(num(TL[I]),L); + } + } + } + if(F==-3) return [Var,L]; + for(I=0;I=0 && IT=0&&!chkfun("gr",0)){ + mycat("load(\"gr\"); /* <- do! */\n"); + F=-1; + } + if(F==-2) return [VV,L]; + if(F<0) return [V,L]; + 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{ + P=LL[length(LL)-1]; + for(RV=reverse(V), I=0;I1&&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); @@ -8533,13 +8797,13 @@ def mcgrs(G, R) { NP = length(G); Mat = (getopt(mat)==1)?0:1; - if(Mat==1 && type(SM=getopt(slm))==4){ + if(Mat==0 && type(SM=getopt(slm))==4){ SM0=SM[0];SM1=anal2sp(SM[1],["*",-1]); if(findin(0,SM0)>=0){ for(SM=[],I=length(G)-1;I>0;I--) if(findin(I,SM0)<0) SM=cons(I,SM); SM=[SM,SM1]; - G=mcgrs(G,R|slm=SM); + G=mcgrs(G,R|mat=1,slm=SM); return [G[0],anal2sp(G[1],["*",-1])]; } }else SM0=0; @@ -8548,14 +8812,15 @@ def mcgrs(G, R) L = length(G)-1; RT = car(R); if(type(RT) == 4){ - if(length(RT)==L&&RT[0]!=0){ - R=cons(cdr(RT),R); + if(length(RT)==L+1&&RT[0]!=0){ + R=cons(cdr(RT),cdr(R)); R=cons(RT[0],R); + R=cons(0,R); continue; } /* addition */ RT = reverse(RT); S = ADS = 0; - for(G = reverse(G); G != []; G = cdr(G), L--){ - AD = car(RT); RT = cdr(RT); + for(G = reverse(G); G != []; G = cdr(G), L--, RT=cdr(RT)){ + AD = car(RT); if(L > 0){ S += AD; if(SM && findin(L,SM0)>=0) ADS+=AD; @@ -8582,6 +8847,7 @@ def mcgrs(G, R) OD += car(GT)[0]; if(car(GT)[1] == RTT && car(GT)[0] > M){ S += car(GT)[0]-M; + M=car(GT)[0]; VP[I] = J; } } @@ -8590,7 +8856,7 @@ def mcgrs(G, R) for(GN = []; L >= 0; L--){ GT = GV[L]; RTT = (L==0)?(-RT):RT; - GTN = [(VP[L] >= 0 || S == 0)?[]:[-S,(L==0)?(Mat-RT):0]]; + GTN = (VP[L]>=0 || S == 0)?[]:[[-S,(L==0)?(Mat-RT):0]]; for(J = 0; GT != []; GT = cdr(GT), J++){ if(J != VP[L]){ GTN = cons([car(GT)[0],car(GT)[1]+RTT], GTN); @@ -8601,31 +8867,70 @@ def mcgrs(G, R) print("Not realizable"); return; } - GTN = cons([K,(L==0)?(Mat-RT):0], GTN); + if(K>0) GTN = cons([K,(L==0)?(Mat-RT):0], GTN); } - if(VP[L]<0) GTN=cons([-S,(L==0)?(Mat-RT):0],GTN); GN = cons(reverse(GTN), GN); } - if(SM0){ - for(M=0,L=length(G)-1;L>0;L--) - if(findin(L,SM0)>=0&&VP[L]>=0) M+=GV[L][VP[L]][0]; - Mx=sp2anal(SM1,["max",1,0]); - for(SM2=[],J=0;SM1!=[];J++,SM1=cdr(SM1)){ - if(J!=Mx[0]) SM2=cons([car(SM1)[0],car(SM1)[1]+RT],SM2); - else if((V=car(SM1)[0]-M)!=0) SM2=cons([V,RT],SM2); + if(SM0&&RT!=0){ + for(M0=M1=-OD,L=length(G)-1;L>=0;L--){ + if(findin(L,SM0)>=0){ + M0+=OD; + if(VP[L]>=0) M0-=GV[L][VP[L]][0]; + }else{ + M1+=OD; + if(VP[L]>=0) M1-=GV[L][VP[L]][0]; + } } - Mx=sp2anal(SM1,["max",1,0]); - for(J=0;SM2!=[];J++,SM2=cdr(SM2)){ - if(J!=Mx[0]) SM1=cons(car(SM2),SM1); - else if((V=car(SM1)[0]-S+M)!=0) SM1=cons([V,0],SM2); + SM2=[]; + if((Mx1=anal2sp(SM1,["max",1,-RT])[0])<0){ + if(M1>0) SM2=cons([M1,0],SM2); + }else M1+=car(SM1[Mx1]); + if((Mx0=anal2sp(SM1,["max",1,0])[0])<0){ + if(M0>0) SM2=cons([M0,RT],SM2); + }else M0+=car(SM1[Mx0]); + for(J=0;SM1!=[];J++,SM1=cdr(SM1)){ + if(J==Mx0){ + if(M0>0) SM2=cons([M0,-RT],SM2); + }else if(J==Mx1){ + if(M1>0) SM2=cons([M1,0],SM2); + }else SM2=cons([car(SM1)[0],car(SM1)[1]+RT],SM2); } - if(Mx[0]<0) SM1=cons([-S,0],SM1); + SM1=reverse(SM2); } G = cutgrs(GN); } - return SM0?[G[0],SM1]:G; + return SM0?[G,SM1]:G; } +def spslm(M,TT) +{ + R=getbygrs(M,1|mat=1); + if(type(R)!=4||type(R[0])!=4||type(S=R[0][1])!=4){ + errno(0);return0; + } + if(S[1]!=[[1,0]]){ + print("Not rigid!");return0; + } + if((F=S[0][0][1])!=0){ + for(V=vars(F);V!=[];V=cdr(V)){ + if(mydeg(F,car(V))==1){ + T=lsol(F,car(V)); + break; + } + } + if(V==[]){ + print("Violate Fuchs condition!"); + return0; + } + } + for(P=[];R!=[];R=cdr(R)) + P=cons(car(R)[0],P); + if(F!=0){ + S=mysubst(S,[car(V),T]);P=mysubst(P,[car(V),T]); + } + return mcgrs(S,P|mat=1,slm=[TT,[[1,0]]]); +} + /* F=0 : unify F=["add",S] : @@ -9219,7 +9524,7 @@ def mc2grs(G,P) T=(K==6)?["reduction"]:[]; S=cons(append([x0,x1,x2,x3,x4,"idx"],T),S); M=ltotex(S|opt="tab",hline=[0,1,z], - vline=(K==6)?[0,1,z-2,z-1,z]:[0,1,z-2,z-1,z], + vline=(K==6)?[0,1,z-2,z-1,z]:[0,1,z-1,z], left=["","$x_0$","$x_1$","$x_2$","$x_3$","$x_4$"]); if(Dvi>0) dviout(M|keep=Keep); } @@ -9596,7 +9901,7 @@ def mcmgrs(G,P) L=cons(TL,L); } if(Dvi){ - if(Dvi!=-1) dviout(S|eq=0); + if(Dvi!=-1) dviout(S|eq=0,keep=Keep); return S; } return reverse(L); @@ -17211,8 +17516,294 @@ def conf1sp(M) return P; } +/* ((1)(1)) ((1)) 111|11|21 [[ [2,[ [1,[1]],[1,[1]] ]], [1,[[1,[1]]]] ]] */ +/* (11)(1),111 111|21,111 [[[2,[1,1]],[1,[1]]],[1,1,1]] */ +def s2csp(S) +{ + if(type(S)!=7){ + U=""; + if(type(N=getopt(n))>0){ + for(D=0,S=reverse(S);S!=[];S=cdr(S),D++){ + if(D) U=","+U; + T=str_subst(rtostr(car(S)),","," "); + U=str_cut(T,1,str_len(T)-2)+U; + } + V=strtoascii(U); + for(R=[];V!=[];V=cdr(V)){ + if((CC=car(V))==91){ /* [ */ + if(length(V)>1 && V[1]==91) V=cdr(V); + for(I=1;(CC=V[I])!=91&&CC!=93;I++); + if(CC==91){ + R=cons(40,R); /* ( */ + while(I--) V=cdr(V); + }else{ + V=cdr(V); + while(--I) R=cons(car(V),R); + } + }else if(CC==93){ /* ] */ + R=cons(41,R); + if(length(V)>1 && V[1]==93) V=cdr(V); + }else R=cons(CC,R); + } + return asciitostr(reverse(R)); + } + for(;S!=[];S=cdr(S)){ + if(U!="") U=U+","; + for(D=0,TU="",T=car(S);T!=[];D++){ + if(type(car(T))==4){ + R=lpair(T,0); + T=R[0];R1=m2l(R[1]|flat=1); + }else R1=[]; + if(D) TU="|"+TU; + TU=s2sp([T])+TU; + T=R1; + } + U=U+TU; + } + 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){ /* , */ + + } + } + } + } + for(P=TS=[],I=D=0; S!=[]; S=cdr(S)){ + if((C=car(S))==44){ /* , */ + P=cons(D,P);D=0; + }else if(C==124){ /* | */ + D++;C=44; + } + TS=cons(C,TS); + } + S=reverse(TS); + P=reverse(cons(D,P)); + U=s2sp(asciitostr(S)); + + for(R=[];P!=[];P=cdr(P),U=cdr(U)){ + D=car(P);R0=car(U); + while(D--){ + U=cdr(U); + for(U0=car(U),R2=[];U0!=[];U0=cdr(U0)){ + for(R1=[],N=car(U0);N>0;R0=cdr(R0)){ + R1=cons(car(R0),R1); + if(type(car(R0))==4) N-=car(R0)[0]; + else N-=car(R0); + } + R2=cons([car(U0),reverse(R1)],R2); + } + R0=reverse(R2); + } + R=cons(R0,R); + } + return reverse(R); +} + + +def partspt(S,T) +{ + if(length(S)>length(T)) return []; + if(type(Op=getopt(opt))!=1) Op=0; + else{ + 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; + else return []; + }else if(getopt(sort)==1){ + S0=S1=[]; + for(;S!=[]&&car(S)==car(T);S=cdr(S),T=cdr(T)) + S0=cons(car(S),S0); + if(S!=[]&&car(S)0){ + if(X=Y) R=cons(vtol(VT),R); + if(!vnext(VT)) break; + } + } + if(Op){ + for(W=[];R!=[];R=cdr(R)){ + for(I=0,S=VS[0],K=U=[],TR=car(R);TR!=[];TR=cdr(TR)){ + K=cons(car(TR),K); + if(!(S-=car(K))){ + U=cons([VS[I],reverse(K)],U); + K=[]; + S=VS[++I]; + if(I==L){ + U=cons([S,cdr(TR)],U); + break; + } + } + } + W=cons(reverse(U),W); + } + R=W; + if(iand(Op,1)){ + for(R=[];W!=[];W=cdr(W)) + R=cons(reverse(qsort(car(W))),R); + R=lsort(R,[],1); + } + if(Op==3){ + for(W=[];R!=[];R=cdr(R)){ + for(S=[],TR=car(R);TR!=[];TR=cdr(TR)) + S=append(S,car(TR)[1]); + W=cons(S,W); + } + R=reverse(W); + } + } + return R; +} + +#if 0 +def confspt(S,T) +{ + R=[]; + LS=length(S);LT=length(T); + if(LScar(TT)) return R; + if(car(ST)==car(TT){ + S0=cons(car(ST));T0=cons(car(TT)); + LS--;LT--;continue; + } + V=car(TT);D=LS-LT; + for(P=[ST],DD=D;DD>0;){ + VD=V-car(car(ST)); + } + } +} +#endif + +def mcvm(N) +{ + X=getopt(var); + if((Z=getopt(z))!=1) Z=0; + if(type(N)==4){ + if((K=length(N))==1&&isvar(X)) X=[X]; + if(type(X)!=4){ + for(X=[],I=0;I=N[1]+N[3]) return 0; + X=X[0]; + for(R=[],I=1;I0;I--) V=cons(makev(["a0",I]),V); + MI=myinv(M); + V=ltov(V)*MI; + for(R=[],I=0;I