===================================================================
RCS file: /home/cvs/OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v
retrieving revision 1.70
retrieving revision 1.90
diff -u -p -r1.70 -r1.90
--- OpenXM/src/asir-contrib/packages/src/os_muldif.rr	2020/08/23 00:39:59	1.70
+++ 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.69 2020/05/17 23:15:26 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 - Aug. 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,15 +82,21 @@ 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$
 localf dvprod$
 localf dnorm$
+localf dext$
 localf mulseries$
 localf pluspower$
 localf vtozv$
@@ -97,12 +104,16 @@ localf dupmat$
 localf matrtop$
 localf mytrace$
 localf mydet$
+localf permanent$
 localf mperm$
 localf mtranspose$
 localf mtoupper$
 localf mydet2$
 localf myrank$
+localf lext2$
 localf meigen$
+localf pf2kz$
+localf mext2$
 localf transm$
 localf vgen$
 localf mmc$
@@ -117,10 +128,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$
@@ -160,6 +175,8 @@ localf execdraw$
 localf execproc$
 localf myswap$
 localf mysubst$
+localf sort2$
+localf n2a$
 localf evals$
 localf myval$
 localf myeval$
@@ -182,6 +199,8 @@ localf mylog$
 localf nlog$
 localf mypow$
 localf scale$
+localf catalan$
+localf iceil$
 localf arg$
 localf sqrt$
 localf gamma$
@@ -209,6 +228,7 @@ localf mmulbys$
 localf appldo$
 localf appledo$
 localf muldo$
+localf caldo$
 localf jacobian$
 localf hessian$
 localf wronskian$
@@ -345,6 +365,8 @@ localf s2euc$
 localf s2sjis$
 localf r2ma$
 localf evalma$
+localf evalcoord$
+localf readTikZ$
 localf ssubgrs$
 localf verb_tex_form$
 localf tex_cuteq$
@@ -363,6 +385,7 @@ localf getbyshell$
 localf show$
 localf dviout$
 localf rtotex$
+localf togreek$
 localf mtotex$
 localf ltotex$
 localf texbegin$
@@ -375,7 +398,9 @@ localf shiftPfaff;
 localf conf1sp$
 localf confexp$
 localf confspt$
+localf vConv$
 localf mcvm$
+localf s2cspb$
 localf s2csp$
 localf partspt$
 localf pgen$
@@ -420,6 +445,7 @@ localf openGlib$
 localf xyproc$
 localf xypos$
 localf xyput$
+localf xylabel$
 localf xybox$
 localf xyline$
 localf xylines$
@@ -449,15 +475,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$
@@ -506,7 +538,7 @@ extern AMSTeX$
 extern Glib_math_coordinate$
 extern Glib_canvas_x$
 extern Glib_canvas_y$
-Muldif.rr="00200812"$
+Muldif.rr="00220213"$
 AMSTeX=1$
 TeXEq=5$
 TeXLim=80$
@@ -677,7 +709,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";
@@ -825,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;
@@ -1045,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
@@ -1130,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);
@@ -1144,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;
 		}
 	}
@@ -1159,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;
@@ -1174,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);
@@ -1241,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);
@@ -1249,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)
@@ -1275,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);
@@ -1308,6 +1398,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;
@@ -1343,6 +1434,8 @@ 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];
@@ -1391,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;
@@ -1730,6 +2481,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<S;I++){
+		if(!(K=M[I][0])) continue;
+		for(TL=[],SL=LL,J=0;J<S;J++,SL=cdr(SL))
+			if(I!=J) TL=cons(cdr(car(SL)),TL);
+		if(K) V+=K*permanent(lv2m(TL));
+	}
+	return V;
+}
+
 def mperm(M,P,Q)
 {
 	if(type(M) == 6){
@@ -2224,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);
@@ -2364,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);
@@ -2385,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);
@@ -2406,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;
@@ -2417,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);
@@ -2465,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;
 }
 
@@ -2737,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;
@@ -3032,6 +3919,7 @@ def mdsimplify(L)
 	return L;
 }
 
+#if 1
 def m2mc(M,X)
 {
 	if(type(M)<2){
@@ -3122,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"){
@@ -3219,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)
 {
@@ -3461,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];
@@ -3479,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))
@@ -3488,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){
@@ -4052,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)
@@ -4493,7 +5610,8 @@ def fctrtos(P)
 					if(L>0)	str_tb(CR,Out);
 					if(LL>Lim){
 						if(TOC==7)	OC=texlim(OC,Lim|cut=[CR,CR2]);
-						PW+=CR; L=0;
+						if(length(Tm)!=1) PW+=CR; 
+						L=0;
 					}else L=LL;
 				}else L+=LL;
 			}else if(length(Tm)!=1){
@@ -5172,7 +6290,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++;
@@ -5381,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;
 	}
@@ -5426,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;
@@ -5550,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);
 }
@@ -6568,13 +7716,22 @@ def baseODE(L)
 	}
 	if(type(To=getopt(to))<2||type(To)>4) To=0;
 	if(Ord<0){ 	/* cancell y1, z1,... by baseODE0() */
-		if(!++Ord) Ord=2;
+		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];
-			S0=baseODE(L0|TeX=1,f=-1);
-			V=baseODE0(L|step=-1,to=To);
+			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);
@@ -6624,7 +7781,7 @@ def baseODE(L)
 			}
 		}
 	}
-	if(F==-3) return [Var,L];
+	if(F==-3&&!TeX) return [Var,L];
 	for(I=0;I<M;I++) L=subst(L,Var[I],makev([SV[I],0]));
 	if(TeX){
 		for(TL=L,I=0;I<M;I++) 
@@ -6752,8 +7909,8 @@ def eqs2tex(P,L)
 	if(type(L)!=4) L=[];
 	Sgn=0;
 	if(L!=[]){
-		if(car(L)==0) L=[[L]];
-		else if(length(L)>1 && isvar(L[1])) L=[[L]];
+		if(car(L)==0) L=[L];
+		else if(length(L)>1 && isvar(L[1])) L=[L];
 		R=car(L);L=cdr(L);Sgn=1;
 	}else R=[];
 	if(type(R)==4&&car(R)==0){
@@ -6771,6 +7928,7 @@ def eqs2tex(P,L)
 			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;
 	}
@@ -6791,14 +7949,28 @@ def eqs2tex(P,L)
 	return S;
 }
 
-
+/* Opt: var, opt, dbg */
 def res0(P,Q,X)
 {
-	if(!isvar(X)) return -1;
-	if(type(Var=getopt(var))!=4) Var=0;
-	if(!isint(Opt=getopt(opt))) Opt=0;
+	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))<nmono(coef(Q,DQ,X))){
+		R=P;P=Q;Q=R;
+		R=DP;DP=DQ;DQ=R;
+	}
 	while(DQ>0){
 		if(DP<DQ){
 			R=P;P=Q;Q=R;
@@ -6806,37 +7978,91 @@ def res0(P,Q,X)
 			if(Opt==-1) return [P,Q,DP,DQ];
 			if(DQ<1) break;
 		}
+		if(Dbg){
+			if(Dbg>=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);
-		T=gcd(TP,TQ);
-		P=red(TQ/T)*P-red(TP/T)*Q*X^(DP-DQ);
-		if(Var){
-			for(S=fctr(P),P=1;S!=[];S=cdr(S)){
-				TV=vars(car(S)[0]);
-				if(type(TV)==4&&lsort(TV,Var,2)!=[]) P*=car(S)[0];
+		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) return [P,Q,DP,DQ];
+		if(Opt==-2||(type(Opt)==4&&Opt[0]==DP&&Opt[1]==DQ)) return [P,Q,DP,DQ,W];
 	}
-	if(Opt==1) Q=[P,Q,DP,DQ];
+	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 */
+/* Opt : f, var, ord, ord, step, f, to */
 def baseODE0(L)
 {
-	if(!isint(Ord=getopt(ord))) Ord=0;
-	if(!isint(Step=getopt(step))) Step=0;
+	if(!isint(Ord=getopt(ord))) Ord=-1;
+	if(Ord==-1) Ord=2;
+	if(Ord<O) Ord++;
 	if(!isint(F=getopt(f))) F=0;
+	if(!isint(Dbg=getopt(dbg))) Dbg=0;
+	if(type(Step=getopt(step))==4) Dstep=Step;
+	else Dstep=0;
+	if(!isint(Step)) Step=0;
 	if(F<0) Step=1;
 	if(Step>0&&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=0;V=0;
-	}
+	if(!isvar(To)) To=V=0;
 	if(type(SV=Var=getopt(var))!=4){
 		SV=SVORG;
 		if(N>10){
@@ -6878,9 +8104,8 @@ def baseODE0(L)
 			TR=R=reverse(R);
 			if(length(R)>1){	/* reduce common factor */
 				P=car(TR);TR=cdr(TR);
-				for(;TR!=[]&&P!=1;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);
@@ -6894,27 +8119,35 @@ def baseODE0(L)
 				TR=cons(car(R),TR);
 				R=cdr(R);
 			}
-			R0=cons(car(R),R0);
+			R0=(F==2)?append(R,R0):cons(car(R),R0);
 			if(R!=[]){
-				for(P=car(R),R=cdr(R); R!=[]; R=cdr(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){
-						TQ=res0(P,car(R),TV|var=V0,opt=1);
+						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);
-						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];
+						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;
 		}
-		return (F==1)?car(R):cons(car(R),R0);
+		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];
@@ -7069,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);
 	}
@@ -7080,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);
@@ -7681,12 +8914,14 @@ def pmaj(P)
 		return S;
 	}
 	V=vars(P);
-	if(!(K=length(V))) return abs(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):abs(Q))*X^D;
+		if(Q!=0) R+=((type(Q)>1)?pmaj(Q|var=Abs):(Y==1?1:abs(Q)))*X^D;
 	}
-	if(isvar(Y=getopt(var))) for(;V!=[];V=cdr(V)) R=subst(R,car(V),Y);
+	if(isvar(Y)) for(;V!=[];V=cdr(V)) R=subst(R,car(V),Y);
 	return R;
 }
 
@@ -8089,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;
 	}
@@ -8101,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];
@@ -8114,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);
@@ -8134,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;
 }
@@ -8338,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];
@@ -8359,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++){
@@ -8540,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;
 }
 
@@ -9939,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);
 					}
 				}
@@ -10728,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)){
@@ -11203,6 +12444,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 || J<I) return [0,[]];
+	for(F=1,K=I+1;K<J;K++){
+		C=L[K];
+		if(C>32&&(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();
@@ -12020,7 +13387,7 @@ def show(P)
 		if(N==1){
 			if(type(Var=getopt(var))>1){
 				if(isvar(Var)) Var=[0,Var];
-				else if(type(Var)==4) Var=cons(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=[];
@@ -12103,13 +13470,20 @@ 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)){
-				S=P+"\n\n";
-				if(Raw) return S;
-				dviout(S); 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;
+		}
 	}
 	if(Raw) return "\\begin{align}\\begin{split}\n &"+S+"\\end{split}\\end{align}";
 	else dviout(S|eq=5);
@@ -12312,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;	*/
@@ -12542,54 +13935,86 @@ def frac2n(N)
 #endif
 }
 
+/* Option : opt */
 def ptconvex(L)
 {
-	if(!(Opt=isint(getopt()))) Opt=0;
-/*
-	if(Opt==1){
-		P1=P2=P3=P4=car(L)[0];
-		for(TL=cdr(L);TL!=[];TL=cdr(TL)){
-		}
-	}
-*/
-	L0=car(L);X=car(L0);Y=car(L0)[1];
+	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(X>car(TL)[0]||(X==car(TL)[0]&&Y>car(TL)[1])) continue;
+		if(X<car(TL)[0]||(X==car(TL)[0]&&Y<car(TL)[1])) continue;
 		L0=car(TL);X=car(L0);
 	}
+	if(Opt==3) return L0;
 
 	R=[];	/* find a polygone through all points */
 	X0=L0[0];Y0=L0[1];
 	for(TL=L;TL!=[];TL=cdr(TL)){
-		X=car(TL)[0]-X0;Y=car(TL)[1]-Y0;S=X^2+Y^2;
-		if(!S) R=cons(-2,cons(0,R));
-		else R=cons((Y>0?Y^2:-Y^2)/S,cons(S,R));
+		L0=car(TL);
+		X=L0[0]-X0;Y=L0[1]-Y0;S=X^2+Y^2;
+		L0=(!S)? append([-8,0],L0):append([(Y>0?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))
-		R=cons(cdr(cdr(car(TL))),R);
-	R=reverse(R);	
-	if(Opt==1) return R;
 
-	X0=car(L)[0];Y0=car(L)[1];R=[cons(V0=-4,car(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)){
-		X=(X1=car(TL)[0])-X0;Y=(Y1=car(TL)[1])-Y0;S=X^2+Y^2;
-		if(!S) continue;
-		V=(Y>0?Y^2:-Y^2)/S;
+		V=darg(L0,L1=car(TL));
+		if(V<-4) continue;
 		while(V<V0){
 			R=cdr(R);
-			X=X1-car(R)[1];Y=Y1-car(R)[2];S=X^2+Y^2;V0=cdr(R)[0];
-			V=(Y>0?Y^2:-Y^2)/S;
-			if(X<0) V=2-V;		/* -1 - (3-0)  */
+			V0=car(car(R));
+			V=darg(cdr(car(R)),L1);
 		}
 		if(V==V0) R=cdr(R);
-		R=cons(cons(V0=V,car(L)),R);
+		R=cons(cons(V0=V,L0=L1),R);
 	}
-	for(L=[],TL=R;TL!=[];TL=cdr(TL)) L=[cdr(car(TL)),L];
+	for(L=[],TL=R;TL!=[];TL=cdr(TL)) L=cons(cdr(car(TL)),L);
 	return L;
 }
 
+def darg(P,Q)
+{
+	if(type(car(P))==4){
+		if((V=darg(Q[0],Q[1]))<-1) return -8;
+		if((V-=darg(P[0],P[1]))>2){
+			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="";
@@ -12667,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];
@@ -12679,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";
@@ -14017,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){
@@ -14773,8 +16226,8 @@ def xyplot(L,LX,LY)
 	if(type(Vw)!=1 && type(Vw)!=7 && Vw!=0) Vw=-1;
 	if(!LX){
 		L0=llget(L,1,[0]|flat=1);
-		LX=[lmin(L0),lmax(L0)];
-		S=LX[1]-LX[0];
+		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];
@@ -14782,8 +16235,8 @@ def xyplot(L,LX,LY)
 	LX=map(deval,LX);
 	if(!LY){
 		L0=llget(L,1,[1]|flat=1);
-		LY=[lmin(L0),lmax(L0)];
-		S=LY[1]-LY[0];
+		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];
@@ -14799,6 +16252,27 @@ def xyplot(L,LX,LY)
 		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"]);
@@ -15416,6 +16890,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;
@@ -15458,18 +16986,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;
-			R=T[0]?T[0]:T[1];R0=subst(R,x_,0);R1=subst(R,x_,1);
-			F=subst(R0,y_,0);
-			if(subst(R0,y_,1)*F>0&&subst(R1,y_,0)*F>0&&subst(R1,y_,1)*F>0) return 0;
-			if(X[0]==X[1]) return X[0];
-			if(Y[0]==Y[1]) return Y[0];
-			return 1;
+			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];
+			}
+			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(X0<Y0) X0=Y0;
+			if(Y0>Y1) X1=Y1;
+			if(X0>X1) return 0;
+			if(X0<X1) return [X0,X1];
+			return X0;
 		}else if(Y[1]==0){ /* orth */
 			T=[Y[0][0]+(X[1][1]-X[0][1])*y_-S[0],
 				Y[0][1]-(X[1][0]-X[0][0])*y_-S[1]];
@@ -15526,6 +17065,39 @@ def ptcommon(X,Y)
 	return 0;
 }
 
+
+def ptcontain(P,L)
+{
+	if(type(car(P))==4){
+		if((C=getopt(common))!=1) C=0;
+		if((F0=ptcontain(P[0])&&!C)) return F0;
+		if((F1=ptcontain(P[1])&&!C)) return F1;
+		if(F0&&F1) return P;	/* include */
+		L=cons(L[2],L);		/* outside part exists */
+		for(I=1,R=[];I<4;I++,L=cdr(L)){
+			if(!(F[I]=ptcotain(P,[L[0],L[1]]))){
+				if(C) continue;
+				return -1;
+			}
+			if(type(F[I])==4&&length(F[I])==2)	/* infinite points */
+				return F[I];
+			else R=cons(F[I],R);
+		}
+		if(R==[]) return 0; 		/* no intersection */
+		if(F1==1) return [P[0],car(R)];
+		if(F2==1) return [P[1],car(R)];
+		if(length(R)>1 && 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){
@@ -15935,6 +17507,16 @@ def polstrum(P)
 	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);
@@ -15948,7 +17530,7 @@ def polradiusroot(P)
 			TV=eval((abs(coef(P,TD,X))*ND)^(1/(D-TD)));
 			if(V<TV) V=TV;
 		}
-		return (Int==1)?ceil(V):V;
+		return (Int==1)? iceil(X):X;
 	}
 	for(N0=N1=0,TD=0;TD<D;TD++){
 		if(!(C=coef(P,TD,X))) continue;
@@ -15969,35 +17551,36 @@ def polradiusroot(P)
 			if(V1<TV) V1=TV; 
 		}
 	}
-	return (Int==1)?[-ceil(V1),ceil(V2)]:[-V1,V2];
+	return Int?[-iceil(V1),iceil(V2)]:[-V1,V2];
 }
 
-
+/* step, num, strum */
 def polrealroots(P)
 {
 	if(type(MC=getopt(step))==4){
 		MC1=MC[1];MC=car(MC);
 	}else if(isint(MC)&&MC>1&&MC<10001) MC1=MC;
-	else MC1=MC=20;
-	if(type(I=getopt(num))!=4){
-		I=polradisuroot(P)*eval(eval(exp(0)));
+	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 [];
+	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,I0);
-			V1=subst(P,X,I1);
+			V0=subst(P,X,X0);
+			V1=subst(P,X,X1);
 			while(C++<MC1){
-				V2=subst(P,X,X2=(X0+X1)/2);
+				V2=subst(P,X,X2=(X0+X1)/2*Z);
 				if((V0>0&&V2>0)||(V0<0&&V2<0)) X0=X2;
 				else X1=X2;
 			}
@@ -16005,68 +17588,48 @@ def polrealroots(P)
 			continue;
 		}
 		while(++C<MC){
-			N2=sgnstrum(L,X2=(X0+X1)/2);
+			N2=sgnstrum(L,X2=(X0+X1)/2*Z);
 			if(N0>N2){
-				if(N2>N1) LT=cons([C,X2,X1,N2,N1]);
+				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=-1;
+					C=MC+1;
 				}
 			}else{
 				X0=X2;
 				N0=N2;
 			}
 		}
-		if(C>0) R=cons([X0,X1,N0-N1],R);
+		if(C!=MC+2) R=cons([X0,X1,N0-N1],R);
 	}
-	return reverse(cons(P,R));
-}
-
-#if 0
-def polrealroots(P,X0,W)
-{
-	X=vars(P0=P);
-	if(length(X)!=0) return [];
-	X=car(X);
-	if(isfuctr(P)){
-		Q=gcd(P,diff(P,X));
-		P=sdv(P,Q);
-	}
-	D=deg(P,X);
-	P=P/coef(P,deg(P,X),X);
-	Q=diff(P,X)/D;
-	for(L=[Q,P],I=D;I>=0;I--){
-		R=urem(P,Q);
-		if((TD=deg(R,X))>0){
-			C=coef(R,TD,X);
-			if(C<0) C=-C;
-			R/=C;
+	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)<abs(V1))
+				X0=car(TR)[0];
+			else{
+				X0=car(TR)[1];V0=V1;
+			}
+			for(Tn=Nt;Tn>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);
 		}
-		L=cons(R,L);
+		for(TR=R;TR!=[];TR=cdr(TR))
+			if(car(TR)[2]>1) S=cons(car(TR),S);
+		return reverse(S);
 	}
-mycat(["In",[X0,X0+W]]);
-	Div=10;
-	X=var(P);
-	if(!isint(Step=getopt(step))) Step=6;
-	W=W*eval(exp(0));X1=X0+W;
-	Q=diff(subst(P,X,X+X0+W/2),X);
-	Mx=subst(pmaj(Q),X,W/2);
-	W/=Div;
-	M=Mx*W/2;
-	for(R=[],X0+=W/2;X0<X1;X0+=W){
-		if((V=subst(P,X,X0))<0) V=-V;
-/* mycat([X0,V,M]); */
-		if(V-M<=0) R=cons(X0-W/2,R);
-	}
-	if(R!=[]) mycat(["get",reverse(R)]);
-	if(!Step) return reverse(R);
-	for(L=[];R!=[];R=cdr(R))
-		L=append(L,polrealroots(P,car(R),W|step=Step-1));
-	return L;
+	return reverse(cons(P,R));
 }
-#endif
 
 /*
 def ptcombezier0(P,Q)
@@ -16286,7 +17849,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);
@@ -16305,7 +17868,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!=[]){
@@ -16327,6 +17890,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){
@@ -16494,49 +18058,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);
@@ -16698,6 +18289,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){
@@ -16840,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);
@@ -16914,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;
@@ -18797,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)),","," ");
@@ -18843,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;
@@ -18895,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=[];
@@ -19001,6 +20913,23 @@ 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+J<I+1) return 0;
+	if(K+J==I+1) return 1;
+	else
+#if 1
+	L=I-K<J-2?I-K+1:J;
+	for(S=0,M=0;M<L;M++) S+=(makev([X,K+M])-makev([X,J-M-1]))*vConv(K+M,I,J-M-1|var=X);
+	return S;
+#else
+	return  vConv(K+1,I,J-1|var=X)+(makev([X,K])-makev([X,J-1]))*vConv(K,I,J-1|var=X);
+#endif
+}
+
 def mcvm(N)
 {
   X=getopt(var);
@@ -19008,26 +20937,12 @@ def mcvm(N)
   if(type(N)==4){
     if((K=length(N))==1&&isvar(X)) X=[X];
     if(type(X)!=4){
-      for(X=[],I=0;I<K;I++) X=cons(asciitostr([97+I]),X);
+      for(X=[],I=0;I<K;I++) X=cons(asciitostr([97+I]),X);	/* a,b,... */
       X=reverse(X);
     }
-	if(getopt(e)==1){
-	  if(length(N)==4){
-		N=ltov(N);
-		if(N[1]<N[3]){
-			I=N[1];N[1]=N[3];N[3]=I;
-		}
-		if(N[2]<N[3]||N[2]>=N[1]+N[3]) return 0;
-		X=X[0];
-		for(R=[],I=1;I<N[3];I++) R=cons(makev([X[0],I]),R);
-		for(L=[],I=N[1];I<=N[2];I++) L=cons(makev([X[0],I]),L);
-		for(S=0,I=N[1];I<=N[2];I++){
-		  V=makev([X[0],I]);
-          S+=polbyroot(R,V)/polbyroot(lsort(L,V,1),V);
-		  S=red(S);
-		}
-		return S;
-      }
+	if((E=getopt(e))==1||E==2){
+	  if(length(N)==4) N=cdr(N);
+	  if(length(N)==3) return vConv(N[0],N[1],N[2]|var=X,e=E);
 	}
 	for(M=[],I=S=0;I<K;Z=0,I++){
 		M=cons(mcvm(N[I]|var=X[I],z=Z),M);
@@ -19083,7 +20998,7 @@ def mcvm(N)
     MV=confexp([F0,V]|sym=3);
     RR=newvect(Mx);
     for(K=0;K<Mx;K++) for(RR[K]=0,I=0;I<Mx;I++) RR[K]=map(red,RR[K]+MV[I][K]*A0[I]);
-    RR0=mysubst(RR,[append(cdr(V1),cdr(V2)),vtol(newvect(Mx-2))]|lpair=1);
+	for(RR0=RR,VV=append(cdr(V1),cdr(V2));VV!=[];VV=cdr(VV)) RR0=subst(RR0,car(VV),0);
     RR0=vtol(RR0);
     return (Get==3)?[RR,RR0]:RR0;
   }
@@ -19663,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]);
@@ -19818,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);
@@ -19831,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;
 }
@@ -21255,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;
@@ -21392,7 +23330,7 @@ def linfrac01(X)
 			}
 		}
 	}
-	return L;
+	return reverse(L);
 }