===================================================================
RCS file: /home/cvs/OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v
retrieving revision 1.76
retrieving revision 1.90
diff -u -p -r1.76 -r1.90
--- OpenXM/src/asir-contrib/packages/src/os_muldif.rr	2020/10/13 02:58:16	1.76
+++ OpenXM/src/asir-contrib/packages/src/os_muldif.rr	2022/02/14 08:24:13	1.90
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.75 2020/10/05 03:42:01 takayama Exp $ */
+/* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.89 2022/02/13 07:39:47 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
 */
@@ -6,7 +6,7 @@
 /* #undef USEMODULE */
 
 /*             os_muldif.rr (Library for Risa/Asir) 
- *          Toshio Oshima (Nov. 2007 - Oct. 2020)
+ *          Toshio Oshima (Nov. 2007 - Feb. 2022)
  *
  *   For polynomials and differential operators with coefficients
  *   in rational funtions (See os_muldif.pdf)
@@ -61,6 +61,7 @@ localf mycoef$
 localf mydiff$
 localf myediff$
 localf mypdiff$
+localf difflog$
 localf pTaylor$
 localf pwTaylor$
 localf m2l$
@@ -81,10 +82,15 @@ localf ndict$
 localf nextsub$
 localf nextpart$
 localf transpart$
+localf getCatalan$
+localf pg2tg$
+localf pgpart$
+localf xypg2tg;
 localf trpos$
 localf sprod$
 localf sinv$
 localf slen$
+localf sexps$
 localf sord$
 localf vprod$
 localf dvangle$
@@ -104,7 +110,10 @@ localf mtranspose$
 localf mtoupper$
 localf mydet2$
 localf myrank$
+localf lext2$
 localf meigen$
+localf pf2kz$
+localf mext2$
 localf transm$
 localf vgen$
 localf mmc$
@@ -124,6 +133,9 @@ localf lchange$
 localf llsize$
 localf llbase$
 localf llget$
+localf lcut$
+localf rev$
+localf qsortn$
 localf lsort$
 localf rsort$
 localf lpair$
@@ -163,6 +175,8 @@ localf execdraw$
 localf execproc$
 localf myswap$
 localf mysubst$
+localf sort2$
+localf n2a$
 localf evals$
 localf myval$
 localf myeval$
@@ -185,6 +199,7 @@ localf mylog$
 localf nlog$
 localf mypow$
 localf scale$
+localf catalan$
 localf iceil$
 localf arg$
 localf sqrt$
@@ -213,6 +228,7 @@ localf mmulbys$
 localf appldo$
 localf appledo$
 localf muldo$
+localf caldo$
 localf jacobian$
 localf hessian$
 localf wronskian$
@@ -369,6 +385,7 @@ localf getbyshell$
 localf show$
 localf dviout$
 localf rtotex$
+localf togreek$
 localf mtotex$
 localf ltotex$
 localf texbegin$
@@ -383,6 +400,7 @@ localf confexp$
 localf confspt$
 localf vConv$
 localf mcvm$
+localf s2cspb$
 localf s2csp$
 localf partspt$
 localf pgen$
@@ -427,6 +445,7 @@ localf openGlib$
 localf xyproc$
 localf xypos$
 localf xyput$
+localf xylabel$
 localf xybox$
 localf xyline$
 localf xylines$
@@ -459,12 +478,14 @@ 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$
@@ -517,7 +538,7 @@ extern AMSTeX$
 extern Glib_math_coordinate$
 extern Glib_canvas_x$
 extern Glib_canvas_y$
-Muldif.rr="00201012"$
+Muldif.rr="00220213"$
 AMSTeX=1$
 TeXEq=5$
 TeXLim=80$
@@ -835,6 +856,22 @@ def mypdiff(P,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;
@@ -1055,12 +1092,13 @@ def cmpsimple(P,Q)
 
 def simplify(P,L,T)
 {
-	if(type(P) > 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
@@ -1140,6 +1178,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);
@@ -1154,11 +1193,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;
 		}
 	}
@@ -1169,13 +1208,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;
@@ -1184,6 +1224,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);
@@ -1251,6 +1292,10 @@ def transpart(L)
 
 def trpos(A,B,N)
 {
+	if(!N){
+		N=(A<B)?B:A;
+		N++;
+	}
 	S = newvect(N);
 	for(I = 0; I < N; I++)
 		S[I]=(I==A)?B:((I==B)?A:I);
@@ -1259,20 +1304,43 @@ def trpos(A,B,N)
 
 def sprod(S,T)
 {
-	L = length(S);
+	if(F=isint(S)){
+		S=vtol(ldict(S,0));T=vtol(ldict(T,0));
+	}
+	if((L=length(S))==2){
+		S=trpos(S[0],S[1],0);
+		L=length(S);
+	}
+	if((R=length(T))==2){
+		T=trpos(T[0],T[1],0);
+		R=length(T);
+	}
+	if(L!=R){
+		if(L>R){
+			W=newvect(L);
+			for(I=0;I<L;I++) W[I]=(I<R)?T[I]:I;
+			T=W;
+		}
+		else{
+			W=newvect(R);
+			for(I=0;I<R;I++) W[I]=(I<L)?S[I]:I;
+			S=W;
+			L=R;
+		}
+	}
 	V = newvect(L);
-	while(--L >= 0)
-		V[L] = S[T[L]];
-	return V;
+	while(--L >= 0) V[L] = S[T[L]];
+	return (F)?ndict(V):V;
 }
 
 def sinv(S)
 {
+	if(F=isint(S)) S=ltov(ldict(S,0));
 	L = length(S);
 	V = newvect(L);
 	while(--L >= 0)
 		V[S[L]] = L;
-	return V;
+	return (F)?ndict(V):V;
 }
 
 def slen(S)
@@ -1285,6 +1353,18 @@ def slen(S)
 	return V;
 }
 
+def sexps(S)
+{
+	K=length(S);S=ltov(S);
+	for(R=[],I=0;I<K-1;I++){
+		for(J=I;J>=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);
@@ -1404,6 +1484,664 @@ 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(M<N||M<0||N<0) return 0;
+			if(N==0) return 1;
+			T=fac(M+N);
+			return T/fac(M)/fac(N)-T/fac(M+1)/fac(N-1);
+		}
+		if(length(K)==3){
+			T=K[0];N=K[1];K=K[2];
+			if(T=="P"||T=="01"||T=="H"){
+				if(T=="P") return fac(N)/fac(N-K);
+				if(T=="H") N+=K-1;
+				if(K<0||K>N) 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&&!K[J];J++);
+				if(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;J<M;J++){
+				V=cons(0,V);
+				for( ; Y[J]>0; 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;I<N;I++) V[I]+=R[I];
+				return V;
+			}else{
+				for(K="(",I=0;;I){
+					for(J=R[I];J>0;J--) K=K+")";
+					for(J=V[I];J>0;J--) K=K+"(";
+					if(++I<N) K=K+"*";
+					else break;
+				}
+				K=K+")";
+				if(To!="T") K=getCatalan(K,0|to=To);
+				return K;
+			}
+		}
+			/* 01 list */
+		if(To=="s") return asciitostr(calc(X,["+",48]));
+		if(To=="T"||To=="#"||To=="P"){
+			for(R=["*"],TX=X;TX!=[];TX=cdr(TX)){
+				if(car(TX)==0) R=cons("*",R);
+				else{
+					T=car(R);R=cdr(R);T="("+car(R)+T+")";
+					R=cons(T,cdr(R));
+				}
+			}
+			R=car(R);
+			if(To=="P"||To=="#") R=getCatalan(R,0|to=To);
+			return R;
+		}
+		K=length(X)/2;
+		if(K<2) return 0;
+		if(F){			/* Not checked */
+			for(I=F=0,TX=X;TX!=[];I++,TX=cdr(TX)){
+				if(!car(TX)){
+					F++;continue;
+				}
+				F--;
+				if(F<=0){
+					for(J=1;J<I;J+=2) F+=catalan((J-1)/2)*catalan(K-(J+1)/2);
+					J=lcut(X,1,I-1);J=getCatalan(J,(I-1)/2|opt=1);
+					V=lcut(X,I+1,2*K-1);V=getCatalan(V,K-(I+1)/2|opt=1);
+					F+=J*catalan(K-(I+1)/2)+V;
+					break;
+				}
+			}
+			return F;
+		}
+		for(R=M=N=0,TX=X;TX!=[];TX=cdr(TX)){
+			if(car(TX)==0){
+				if(M>N) 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<N;I++){
+			V=catalan(I)*catalan(N-I-1);
+			if(X>=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;I<S;I++){
+				J=I%M;T=ladd(Or0,Or1,J);
+				J=(I-J)/M;T=ladd(T,Or2,-J);
+				R=cons(T,R);
+			}
+			Or=reverse(R);
+		}
+		if(!Pr&&TikZ){
+			Tb=str_tb("%TikZ0%\n",0);
+			if(Line!=[]) F=ior(F,512);
+			for(I=0;K!=[];K=cdr(K),Or=cdr(Or),I++){
+				T=append([["org",car(Or)],["skip",F],["num",I+1]],Opt);
+				Tb=str_tb("%"+rtostr(I+1)+"%"+"\n",Tb); 
+				Tb=str_tb(xypg2tg(car(K)|option_list=T),Tb);
+				F=ior(F,1);
+			}
+			if(length(Line)>0){
+				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;I<S;I++) L[I]=[Or[0]+M*dcos(D+3.1416*2*I/S),Or[1]+M*dsin(D+3.1416*2*I/S)];
+		if(!iand(F,2)) R=[xylines(L|proc=1,close=1)];
+		if(!iand(F,4)){
+			for(TK=K;TK!=[];TK=cdr(TK)) 
+				R=append([xylines([L[car(TK)[0]],L[car(TK)[1]]]|proc=1,opt=Op)],R);
+		}
+		return R;
+	}
+	if(!TikZ) return "";
+	if(Op!="" && strtoascii(Op)[0]!=91) Op="["+Op+"]";
+
+	Tb=str_tb(0,0);
+	J=iand(F,64)?-1:1;
+	if(iand(F,256)){
+		M1=10;M2=6;
+	}else{
+		M1=6;M2=3;
+	}
+	if(!L[0]){
+			for(I=0;I<S;I++) L[I]=["("+rtostr(I+St)+")", 
+				[Or[0]+M*dcos(D+3.1416*2*I*J/S),Or[1]+M*dsin(D+3.1416*2*I*J/S)]];
+	}else{
+		for(I=0;I<S;I++)
+			if(length(L[I])==2) L[I]=["("+rtostr(I+St)+")", L[I][0],L[I],1];
+	}
+	if(!iand(F,1)){
+		for(I=0;I<S;I++) Tb=str_tb("\\coordinate"+(iand(F,256)?"(P"+rtostr(I+St)+")":car(L[I]))
+			+" at "+xypos(L[I][1])+";\n",Tb);
+		if(iand(F,32)){
+			for(I=0;I<S;I++){
+				if(!iand(F,16))  VV=ladd(0,L[I][1],1+V/M);
+				else{
+					VV=ladd(L[I][1],L[(I+1)%S][1],1);
+					J=dnorm(VV);
+					VV=ladd(0,VV,1/2+V/J);
+				}
+				Tb=str_tb("\\coordinate("+(iand(F,16)?"E":"V")+rtostr(I+St)+")"
+					+" at "+xypos(VV)+";\n",Tb);
+			}
+		}
+		Tb=str_tb("%\n",Tb);
+	}
+	if(!iand(F,4)){
+		Tb=str_tb("\\coordinate(S) at "+xypos(Or)+";\n",Tb);
+		if(iand(F,512)) str_tb("\\node(S"+rtostr(Num)+
+			")[circle,minimum size=" +rtostr(2*(M+W))+ "cm] at (S){};\n",Tb);
+		if(iand(F,256)) for(I=0;I<S;I++)  Tb=str_tb("\\coordinate"+L[I][0]+" at " + 
+			"($(S)+(P"+ rtostr(I+St) +")$);\n",Tb);
+		if(Every!=""){
+			Tb=str_tb(Every,Tb);Tb=str_tb("\n",Tb);
+		}
+	}
+	if(!iand(F,2)){
+		Tb=str_tb("\\draw ",Tb);
+		for(I=0;I<S;I++){
+			if(!(I%M1)) Tb=str_tb("\n",Tb);
+			if(I) Tb=str_tb("--",Tb);
+			Tb=iand(F,256)? str_tb(car(L[I]),Tb):str_tb("($(S)+"+car(L[I])+"$)",Tb);
+		}
+		Tb=str_tb("--cycle;\n",Tb);
+	}
+	if(!iand(F,8)){
+		for(I=0,T=K;T!=[];T=cdr(T),I++){
+			if(length(car(T))>2){
+				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<S;I++)
+		Tb=str_tb("\\node at ($(S)+("+(iand(F,16)?"E":"V") + rtostr(I+St) + 
+			")$) \{" +rtostr(I+St) +"\};\n",Tb);
+	S=str_tb(0,Tb);
+	if(Dvi==1) xyproc(S|dviout=1);
+	else if(Dvi==-1) S=xyproc(S);
+	return S;
+}
+
+/* 
+   F=0 : sort
+     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<0 : circulate -F: [I,J] -> {I-F,J-F] 
+   F=[I,J] => another diagonal (flip option)
+   F=[I] : the other ends of diagonal starting from I
+   ["ext",I] 
+   ["res",I] 
+   ["pair",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<S;I++) R=cons(pgpart(K,["ext",I]),R);
+					return R;
+				}
+				if(F[1]=="sym"){
+					R=pgpart(K,1);
+					for(K=[];R!=[];R=cdr(R)) K=cons(pgpart(car(R),["ext",0]),K);
+					for(R=[];K!=[];K=cdr(K)) R=cons(pgpart(car(K),3),R);
+					return R;
+				}
+				if(F[1]) K=pgpart(K,-F[1]);
+				K=cons([0,S-1],K);
+				return pgpart(K,F[1]);
+			}
+			if(F[0]=="res"){
+				F1=F[1];
+				I=(F1-1+S)%S;J=(F1+1)%S;
+				T=sort2([I,J]);
+				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-2) return [];
+				return pgpart(R,0);
+			}
+			if(F[0]=="pair"){
+				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);
+				}
+				return reverse(R);
+			}
+		}
+	}
+	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;I<S;I++){
+			TR=pgpart(K,-I);
+			if(!I) R0=TR;
+			else if(R0==TR) break;
+			R=cons(TR,R);
+		}
+		if(F>1) 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);
+	}
+    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;I<S;I++){
+				if(L[I]||K[I]) continue;
+				for(J=1;J<S;J++) if(K[T0=(I+J)%S]) break;
+				for(J=S-1;J>0;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==8||F=="-#"){
+		for(R=[];K!=[];K=cdr(K)) R=cons(sort2(car(K)),R);
+		return reverse(R);
+	}
+	if(F==12||F=="res"){
+		K=pgpart(K,7);
+		for(I=0,R=[];K!=[];K=cdr(K),I++) if(!K[I]) 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<R?T:R;
+	}
+	if(F=="T"||F=="TT"){
+		if(F=="T") K=asciitostr(K);
+		L=length(K);
+		for(R=[[0]],I=0,N=1;I<L;I++){
+			if(K[I]==ST)
+				R=cons(n2a(N++|opt="[]",s=-1),R);
+			else{
+				TR=append(R[0],[41]);
+				TR=append(R[1],TR);
+				TR=cons(40,TR);
+				R=cons(TR,cdr(cdr(R)));
+			}
+		}
+		return asciitostr(car(R));
+	}
+}
+
+
+def pg2tg(K)
+{
+	if((F=getopt(verb))!=1) F=0;
+	if(!isint(Al=getopt(all))) Al=0;
+	M=getopt(red);
+	if(!isint(M) || M<32) M=1024;
+	if(isint(K)){
+		R=[];
+		if(K>3){
+			while(K-- > 3) R=pg2tg(R|verb=F,red=M,all=Al);
+			return R;
+		}else if(K<-3){
+			for(RR=[],K=-K-3;K>0;K--) RR=cons(R=pg2tg(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;
@@ -2261,6 +2999,75 @@ def meigen(M)
 	return (F==1)?getroot(P,zz|mult=1):getroot(P,zz);
 }
 
+def lext2(L)
+{
+	if(length(L)==2){
+		for(S=0,I=1;I<L[1];I++){
+			S+=L[1]-I;
+			if(S>L[0]) break;
+		}
+		return [I-1,L[0]+L[1]-S];
+	}
+	if(L[0]==L[1]) return [0,0];
+	if(L[0]<L[1]){
+		L0=L[0];L1=L[1];S=1;
+	}else{
+		L0=L[1];L1=L[0];S=-1;
+	}
+	return [(2*L[2]-L0-3)*L0/2+L1-1,S];
+}
+
+def pf2kz(M)
+{
+	L=length(M);
+	for(S=0,N=1;S<=L;N++){
+		S+=N;
+	}
+	if(S!=L+1) return 0;
+	for(L=[],I=0;I<N-2;I++){
+		for(S=J=0;J<N;J++)
+			if(I!=J) S+=M[car(lext2([I,J,N]))];
+		L=cons(-S,L);
+	}
+	for(S=0,K=M;K!=[];K=cdr(K)) S-=car(K);
+	S0=S;
+	for(I=0;I<N-2;I++) S+=M[car(lext2([I,N-2,N]))];
+	L=cons(-S,L);
+	for(S=0,K=L;K!=[];K=cdr(K)) S+=car(K);
+	L=cons(-S,L);L=reverse(L);
+	if(getopt(all)==1){
+		for(LL=[],I=0;I<N-2;I++){
+			for(J=I+1;J<N;J++) LL=cons(M[car(lext2([I,J,N]))],LL);
+			LL=cons(car(L),LL); L=cdr(L);
+		}
+		LL=cons(S0,LL); LL=cons(car(L),LL);LL=cons(L[1],LL);
+		return reverse(LL);
+	}
+	return cons(S0,L);
+}
+
+def mext2(M)
+{
+	S=size(M);
+	if(S[0]!=S[1]) return 0;
+	S=S[0];
+	SS=S*(S-1)/2;MM=matrix(SS,SS);
+	for(I=0;I<S;I++){
+		for(J=I+1;J<S;J++){
+			II=lext2([I,J,S])[0];
+			for(I0=0;I0<S;I0++){
+				L=lext2([I0,J,S]);
+				MM[II][L[0]]+=L[1]*M[I][I0];
+			}
+			for(J0=0;J0<S;J0++){
+				L=lext2([I,J0,S]);
+				MM[II][L[0]]+=L[1]*M[J][J0];
+			}
+		}
+	}
+	return MM;
+}
+
 def transm(M)
 {
 	if(type(M)!=6) M=s2m(M);
@@ -2401,15 +3208,46 @@ def vgen(V,W,S)
 
 def mmc(M,X)
 {
+	if(getopt(full)==1){
+		M=mmc(M,X|option_list=delopt(getopt(),"full"));
+		if(type(M)<4) return -1;
+		L=length(M);
+		Mt=getopt(mult);
+		if((L>=6 && Mt!=0)||(L==3&&Mt==1)){
+			for(SS=2,I=3; I<L; I+=(++SS));
+			if(I==L) Mt=1;
+			else Mt=0;
+		}
+		if(Mt!=1){
+			for(R=[],I=S=0;I<L;I++){
+				S=radd(S,M[I]);
+				R=cons([[0,I+1],M[I]],R);
+			}
+			R=cons([[0,I+1],-S],R);
+			return reverse(R);
+		}
+		for(R=[],I=S=0;I<SS;I++)
+			for(J=I+1;J<=SS;J++,S++) R=cons([[I,J],M[S]],R);
+		for(I=0;I<=SS;I++){
+			for(J=S=0;J<=SS;J++){
+				if(I==J) continue;
+				S=radd(S,delopt(R,(I<J)?[I,J]:[J,I]|get=1));
+			}
+			R=cons([[I,SS+1],-S],R);
+		}
+		return qsort(R);
+	}
+
 	Mt=getopt(mult);
 	if(type(M)==7) M=s2sp(M);
-	if(type(M)!=4) return 0;
+	if(type(M)!=4&&type(M)!=5) return 0;
 	if(type(M[0])<=3){
+		if(type(M)==5) M=vtol(M);
 		for(RR=[];M!=[];M=cdr(M)) RR=cons(mat([car(M)]),RR);
 		M=reverse(RR);
 	}
 	if(type(M[0])!=6){			/* spectre type -> 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);
@@ -2422,17 +3260,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;I<SS;I++) M=cons(MZ,M);
 		for(RR=R; RR!=[]; RR=cdr(RR)){
 			RT=car(RR)[0];
 			if(type(RT)==4){
 				if(RT[0]!=0) M=mmc(M,[RT[0]]|simplify=Simp);
-				M=mmc(M,[cdr(RT)]);
+				for(TT=cdr(RT);TT!=[];TT=cdr(TT)){
+					if(car(TT)!=0){
+						M=mmc(cdr(M),cdr(RT));
+						break;
+					}
+				}
 			}
 		}
-/*		for(R=cdr(R);R!=[];R=cdr(R)) M=mmc(M,[car(R)[0]]|mult=Mt); */
 	}
 	if(X==0) return M;
 	L=length(M);
@@ -2443,9 +3285,13 @@ def mmc(M,X)
 	}else{
 		SS=L;Mt=0;
 	}
+	if(type(X[0])==4){
+		for(;X!=[];X=cdr(X)) M=mmc(M,car(X));
+		return M;
+	}
 	if(length(X)==SS+1){
-		if(car(X)!=0&&(M=mmc(M,[car(X)]|mult=Mt))==0) return M;
-		return mmc(M,cdr(X)|mult=Mt);
+		if(car(X)!=0) M=mmc(M,[car(X)]|simplify=Simp);
+		return mmc(M,cdr(X));
 	}
 	for(I=X;I!=[];I=cdr(I)) if(I[0]!=0) break;
 	if(I==[]) return M;
@@ -2454,10 +3300,10 @@ def mmc(M,X)
 	N=newvect(L);
 	for(I=0;I<L;I++) N[I]=dupmat(M[I]);
 	S=size(N[0])[0];
-	if(type(X)==4&&length(X)>SS){	/* addition */
-		for(I=0;I<SS;I++,X=cdr(X)) if(X[I] != 0) N[I] = radd(N[I],car(X));
+	if(type(X)==4&&length(X)>=SS){	/* addition */
+		for(I=0;I<SS;I++,X=cdr(X)) if(car(X) != 0) N[I] = radd(N[I],diagm(S,[car(X)]));
 	}
-	if(length(X)!=1) return 0;
+	if(length(X)!=1||!X[0]) return N;
 	X=X[0];
 	MZ = newmat(S,S);
 	MM = newvect(L);
@@ -2502,7 +3348,11 @@ def mmc(M,X)
 	if(length(KE) == 0) return MM;
 	KK = mtoupper(lv2m(KE),0);
 	for(I=0;I<L;I++) MM[I] = mmod(MM[I],KK);
-	if(Simp!=0) MM = mdsimplify(MM|type=Simp);
+	if(Simp!=0){
+		MM = mdsimplify(MM|type=Simp,show=1);
+		if(getopt(verb)) show([size(MM[0][0]),MM[1]]);
+		MM=MM[0];
+	}
 	return MM;
 }
 
@@ -2774,7 +3624,7 @@ def mdivisor(M,X)
 			P=M[0][0]; M[0][0]=1;
 			for(J=0;J<S1;J++){	/* (1,1) -> 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;
@@ -3069,6 +3919,7 @@ def mdsimplify(L)
 	return L;
 }
 
+#if 1
 def m2mc(M,X)
 {
 	if(type(M)<2){
@@ -3159,8 +4010,11 @@ def m2mc(M,X)
 			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")
+				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"){
@@ -3256,9 +4110,209 @@ def m2mc(M,X)
 	KK = mtoupper(lv2m(KE),0);
 	for(I=0;I<5;I++)
 		MM[I] = mmod(MM[I],KK);
+	if(Simp!=0){
+		MM = mdsimplify(MM|type=Simp);
+		if(getopt(verb)) show([size(MM[0][0]),MM[1]]);
+		MM=MM[0];
+	}
+	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)
 {
@@ -3498,6 +4552,12 @@ def llbase(VV,L)
 		X = var(L[J]); N = deg(L[J],X);
 		for(I = LV; I < S; I++){
 			if((C2=coef(V[I],N,X)) != 0){
+				if(type(C2)==1){
+					for(K=I+1;K<S;K++){
+						if(!(C1=coef(V[K],N,X))||type(C1)!=1) continue;
+						if(abs(C2)<abs(C1)) I=K;
+					}
+				}
 				if(I > LV){
 					Temp = V[I];
 					V[I] = V[LV];
@@ -3516,6 +4576,10 @@ def llbase(VV,L)
 	return V;
 }
 
+def rev(A,B){return A>B?-1:(A<B?1:0);}
+
+def qsortn(X) {return qsort(X,os_md.rev);}
+
 def rsort(L,T,K)
 {
 	for(R=[];L!=[];L=cdr(L))
@@ -3525,6 +4589,15 @@ def rsort(L,T,K)
 	return (iand(K,1))? reverse(R):R;
 }
 
+def lcut(L,M,N)
+{
+	for(I=0,R=[];L!=[]&&I<=N;I++,L=cdr(L)){
+		if(I<M) continue;
+		R=cons(car(L),R);
+	}
+	return reverse(R);
+}
+
 def llget(L,LL,LC)
 {
 	if(type(LL)==4){
@@ -4089,7 +5162,14 @@ def ladd(X,Y,M)
 	}
 	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)
@@ -5422,8 +6502,8 @@ def appldo(P,F,L)
 		L = vweyl(L);
 		X = L[0]; DX = L[1];
 		for(I=mydeg(P,DX);I>0;I--){
-			if(!(TP=mycoef(P,D,DX))) continue;
-			P=red(P+TP*(muldo(D^(I-1),F,L)-D^I));
+			if(!(TP=mycoef(P,I,DX))) continue;
+			P=red(P-TP*DX^I+TP*muldo(DX^(I-1),F,L));
 		}
 		return P;
 	}
@@ -5467,6 +6547,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;
@@ -5591,6 +6691,13 @@ def mce(P,L,V,R)
 {
 	L = vweyl(L);
 	X = L[0]; DX = L[1];
+	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);
 }
@@ -7195,6 +8302,7 @@ def fromeul(P,L,V)
 		R += mycoef(P,J,DX)*S;
 	}
 	if(getopt(raw)!=1){
+		R=nm(R);
 		while(mycoef(R,0,X) == 0)
 			R = tdiv(R,X);
 	}
@@ -7206,11 +8314,10 @@ def fromeul(P,L,V)
 def sftexp(P,L,V,N)
 {
 	L = vweyl(L); DX = L[1];
-	P = mysubst(toeul(P,L,V|opt_list=getpt()),[DX,DX+N]);
+	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);
@@ -8217,9 +9324,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;
 	}
@@ -8229,7 +9337,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];
@@ -8242,9 +9350,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);
@@ -8262,7 +9370,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;
 }
@@ -8466,7 +9574,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];
@@ -8487,7 +9594,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++){
@@ -8668,34 +9774,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<V);
-			else if(L0=="=")	X=(X==V);
-			else if(L0==">=")   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<V);
+		else if(L0=="=")	X=(X==V);
+		else if(L0==">=")   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;
 }
 
@@ -10067,11 +11175,12 @@ def mc2grs(G,P)
 			}
 		}
 		if(F=="rest"||F=="eigen"||F=="rest0"||F=="rest1"){
-			if(F!="eigen") G=mc2grs(G,"homog");
+			if((Hg=getopt(homog))!=0) Hg=1;
+			if(F!="eigen"&&Hg) G=mc2grs(G,"homog");
 			if(length(P)==1){
 				for(R=[],I=0;I<4;I++){
 					for(J=I+1;J<5;J++){
-						S=mc2grs(G,[F,[I,J]]);
+						S=mc2grs(G,[F,[I,J]]|homog=Hg);
 						if(S!=[]) R=cons(cons([I,J],S),R);
 					}
 				}
@@ -10856,6 +11965,10 @@ def mcmgrs(G,P)
 
 def delopt(L,S)
 {
+	if(getopt(get)==1){
+		for(;L!=[];L=cdr(L)) if(car(L)[0]==S) return car(L)[1];
+		return [];
+	}
 	if((Inv=getopt(inv))!=1&&Inv!=2) Inv=0;
 	if(Inv&&type(S)==4&&type(car(S))==4){
 		for(R=[];L!=[];L=cdr(L)){
@@ -12573,6 +13686,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;	*/
@@ -12960,6 +14092,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];
@@ -12972,6 +14107,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";
@@ -14310,7 +15470,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){
@@ -15751,6 +16911,39 @@ def dext(P,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;
@@ -16866,57 +18059,74 @@ def xybox(L)
 def xyang(S,P,Q,R)
 {
 	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(type(Ar=getopt(ar))!=1) Ar=0;
-		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);
+		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()); /* 垂線 */
+			}
+			O=pt5center(P,Q,R);			
+			if(S==2) H=P;	/* 外心 */
 			else{
-				ANG=[0.7854,0.5236,1.0472];
-				X=(AR==0)?1.5708:ANG[AR-2];
+				if(S>2) S++; /* 内心，傍心 */
+				H=ptcommon([P,Q],[O[S],0]);
 			}
-			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;
+			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(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){
@@ -17393,8 +18603,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);
@@ -17467,6 +18685,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;
@@ -19350,13 +20637,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;I<K;){
+				J=str_char(L,I,",");
+				if(J<0) J=str_len(L)+1;
+				R=cons(s2cspb(str_cut(L,I,J-1)|sub=1),R);
+				I=J+1;
+			}
+			L=reverse(R);
+		}
+		if(Sub==-1) return L;
+	}
+	if(!Sub){
+		for(R=[];L!=[];L=cdr(L)) R=cons(s2cspb(car(L)|sub=1),R);
+		return reverse(R);
+	}
+	if(type(car(L))<2){
+		if(Sub==1) return L;
+		else for(N=0,TL=L;TL!=[];TL=cdr(TL)) N+=car(TL);
+		return [N,L];
+	}
+	for(R=[];L!=[];L=cdr(L)) R=cons(s2cspb(car(L)|sub=2),R);
+	if(Sub==2){
+		for(N=0,T=R;T!=[];T=cdr(T)) N+=car(T[0]);
+		R=[R,N];
+	}
+	return reverse(R);
+}
+
+/* ((1)(1)) ((1))   111|21|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){
+			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)),","," ");
@@ -19396,27 +20729,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=[];J<L;J++){
-			if(S[J]==72) I=J;			/* ( */
-			else if(S[J]>47&&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;
@@ -19448,18 +20769,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<N;I++){
+		if(S[I][I]) continue;
+		for(J=0;J<N;J++){
+			if(I==J) continue;
+			if(S[I]==S[J]){
+				P[I][I]++; P[J][J]--;
+			}
+			R=partsp(S[I],S[J]|opt=2);
+			if(R!=[]) P[I][J]=R;
+		}
+	}
+	for(TT=[];I=N-1;I--){
+		for(R=[],J=0;J<N;J++){
+			if(I==J) continue;
+			if(S[I][J]!=0) R=cons(J,R);
+		}
+		TT=cons(R,TT);
+	}
+	for(TP=I=0;I<N;I++) if(S[I][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))<S) break;
+			for(TT0=TT[car(T0)];TT0!=[];TT0=cdr(TT0)){
+				TP=cons(cons(car(TT0),T0),TP);
+				F=1;
+			}
+	}
+	
+}
+*/
 
+
 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));
-	}
+    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=[];
@@ -20219,17 +21578,16 @@ def pfrac(F,X)
 
 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]);
@@ -20374,11 +21732,72 @@ def sqrt2rat(X)
 
 def cfrac2n(X)
 {
+	if(((Q=getopt(q))==1||Q==-1)&&isall(os_md.isint,X)){
+		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(isvar(car(X))&&length(X)==4&&isint(X[1])){
+		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){
+			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;
 	}
+	Sg=getopt(neg)==1?-1:1;
 	if(L>1){
 		for(Y=[];L>1;L--){
 			Y=cons(car(X),Y);
@@ -20387,21 +21806,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;
 }
@@ -21811,129 +23230,92 @@ 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(type(T)==4&&(length(T)==K+3||length(T)==2)){
+		for(U=[],I=K+2;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<K;I++) R[I]=X[I];
-		for(I=11,J=5; I<K; I+=J++){
-			R[I]=X[I+1]; R[I+1]=X[I];
-		}
 		return R;
 	}
-	if(T==1){
-		if(D==1){
-			for(R=[];X!=[];X=cdr(X)) R=cons(1-car(X),R);
-			return reverse(R);
+	if(!S) S=getopt(tr);
+	if(type(S)==4&&length(S)==K+3){
+		D=2;
+	}else if(S==1) for(S=[],I=K+2;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;I<K;I++) R[I]=X[I];
-		for(I=10, J=5; I<K; I+=J++){
-			R[I]=X[I+J-2]; R[I+J-2]=X[I];
-		}
 		return R;
-	}
-	if(T==2){
-		if(D==1){
-			for(R=[]; X!=[]; X=cdr(X)) R=cons(red(1/car(X)),R);
-			return reverse(R);
-		}
-		R=newvect(K,[X[5],X[1],X[2],X[6],X[4],X[0],X[3],X[9],X[8],X[7]]);
-		for(I=11;I<K;I++) R[I]=X[I];
-		for(I=10,J=5; I<K; I+=J++){
-			R[I]=X[I+J-1]; R[I+J-1]=X[I];
-		}
-		return R;
-	}
-	if(T==3){
-		if(D==1){
-			T=car(X);
-			for(R=[T],X=cdr(X); X!=[]; X=cdr(X))
-				R=cons(red(T/car(X)),R);
-			return reverse(R);
-		}
-		R=newvect(K,[X[7],X[4],X[2],X[6],X[1],X[9],X[3],X[0],X[8],X[5]]);
-		for(I=10,J=5; I<K; I+=J++){
-			R[I]=X[I+J-1]; R[I+1]=X[I+J-2]; R[I+J-2]=X[I+1]; R[I+J-1]=X[I];
-		}
-		return R;
-	}
-	if(T==-1){
-		if(D==1){
-			return append([X[1],X[2],X[0]],cdr(cdr(cdr(X))));
-		}
-		R=newvect(K,[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]]);
-		for(I=11;I<K;I++) R[I]=X[I];
-		for(I=17,J=5; I<K; I+=J++){
-			R[I]=X[I+1]; R[I+1]=X[I];
-		}
-		return R;
-	}
-	if(T<0){
-		if(D==1){
-			for(R=[],I=0; X!=[]; X=cdr(X),I--){
-				if(I==T){
-					R=cons(X[1],R);
-					R=cons(X[0],R);
-					X=cdr(X);
-				}
-				else R=cons(car(X),R);
-			}
-			return reverse(R);
-		}
-		T=3-T;
-		R=newvect(K);
-		for(I=0;I<K;I++) R[I]=X[I];
-		for(I=10,J=5;J<T;I+=J++);
-		for(II=0; II<J-2; II++){
-			R[I]=X[I+J]; R[I+J]=R[I];
-		}
-		for( ; II<J; II++){
-			R[I]=X[I+J+1]; R[I+J+1]=X[I];
-		}
-		return R;
-	}
-	return 0;
+	}else if(T==1){ /* (x_1=0, x_2=1) : 1 */
+		for(R=[];X!=[];X=cdr(X)) R=cons(1-car(X),R);
+		if(S) S=mperm(S,[[1,2]],0);
+	}else if(T==2){ /* (x_1=0, x_{K+2}=infty) */ 
+		for(R=[]; X!=[]; X=cdr(X)) R=cons(red(1/car(X)),R);
+		if(S) S=mperm(S,[[1,K+2]],0);
+	}else if(T==3){ /* (x_0=x, x_2=1) */ 
+		T=car(X);
+		for(R=[red(1/T)],X=cdr(X); X!=[]; X=cdr(X)) R=cons(red(car(X)/T),R);
+		if(S) S=mperm(S,[[0,2]],0);
+	}else if(T==4){ /* (x_0=x,x_1=0) : 0 */
+		T=car(X);
+		for(R=[red(T/(T-1))],X=cdr(X); X!=[]; X=cdr(X)) R=cons(red((T-car(X))/(T-1)),R);
+		if(S) S=mperm(S,[[0,1]],0);
+	}else if(T==5){ /* (x_2=1,x_3=y) : 2 */
+		T=X[1];
+		for(R=[1/T,red(X[0]/T)],X=cdr(cdr(X));X!=[]; X=cdr(X)) R=cons(red(car(X)/T),R);
+		if(S) S=mperm(S,[[2,3]],0);
+	}else if(T==6){ /* (x_{K+1}=y_{K-1}, x_{K+2}=infty) : K+1 */
+		T=X[K-1];
+		for(R=[];length(X)>1;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 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<L; K++,I+=J++);
-		if(I!=L) return 0;
+	if(type(X)==4){
+		K=length(X)-2;
+		if(type(car(X))==4){
+			for(U=[],I=K+4;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;
@@ -21948,7 +23330,7 @@ def linfrac01(X)
 			}
 		}
 	}
-	return L;
+	return reverse(L);
 }