===================================================================
RCS file: /home/cvs/OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v
retrieving revision 1.81
retrieving revision 1.90
diff -u -p -r1.81 -r1.90
--- OpenXM/src/asir-contrib/packages/src/os_muldif.rr	2021/07/04 22:31:52	1.81
+++ 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.80 2020/11/23 10:23:51 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 - July 2021)
+ *          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,6 +82,10 @@ localf ndict$
 localf nextsub$
 localf nextpart$
 localf transpart$
+localf getCatalan$
+localf pg2tg$
+localf pgpart$
+localf xypg2tg;
 localf trpos$
 localf sprod$
 localf sinv$
@@ -105,7 +110,10 @@ localf mtranspose$
 localf mtoupper$
 localf mydet2$
 localf myrank$
+localf lext2$
 localf meigen$
+localf pf2kz$
+localf mext2$
 localf transm$
 localf vgen$
 localf mmc$
@@ -125,6 +133,9 @@ localf lchange$
 localf llsize$
 localf llbase$
 localf llget$
+localf lcut$
+localf rev$
+localf qsortn$
 localf lsort$
 localf rsort$
 localf lpair$
@@ -164,6 +175,8 @@ localf execdraw$
 localf execproc$
 localf myswap$
 localf mysubst$
+localf sort2$
+localf n2a$
 localf evals$
 localf myval$
 localf myeval$
@@ -186,6 +199,7 @@ localf mylog$
 localf nlog$
 localf mypow$
 localf scale$
+localf catalan$
 localf iceil$
 localf arg$
 localf sqrt$
@@ -214,6 +228,7 @@ localf mmulbys$
 localf appldo$
 localf appledo$
 localf muldo$
+localf caldo$
 localf jacobian$
 localf hessian$
 localf wronskian$
@@ -385,6 +400,7 @@ localf confexp$
 localf confspt$
 localf vConv$
 localf mcvm$
+localf s2cspb$
 localf s2csp$
 localf partspt$
 localf pgen$
@@ -429,6 +445,7 @@ localf openGlib$
 localf xyproc$
 localf xypos$
 localf xyput$
+localf xylabel$
 localf xybox$
 localf xyline$
 localf xylines$
@@ -461,6 +478,7 @@ localf xypoch$
 localf xycircuit$
 localf ptline$
 localf ptcommon$
+localf ptinversion$
 localf ptcontain$
 localf ptcopy$
 localf ptaffine$
@@ -520,7 +538,7 @@ extern AMSTeX$
 extern Glib_math_coordinate$
 extern Glib_canvas_x$
 extern Glib_canvas_y$
-Muldif.rr="00210702"$
+Muldif.rr="00220213"$
 AMSTeX=1$
 TeXEq=5$
 TeXLim=80$
@@ -838,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;
@@ -1144,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);
@@ -1158,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;
 		}
 	}
@@ -1173,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;
@@ -1188,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);
@@ -1255,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);
@@ -1263,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)
@@ -1420,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;
@@ -2277,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);
@@ -2833,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;
@@ -3761,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];
@@ -3779,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))
@@ -3788,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){
@@ -4352,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)
@@ -5730,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;
@@ -5854,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);
 }
@@ -7458,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);
 	}
@@ -7469,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);
@@ -8730,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];
@@ -8751,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++){
@@ -8932,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;
 }
 
@@ -13263,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";
@@ -14601,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){
@@ -16042,6 +16911,39 @@ def dext(P,Q)
 	return P[0]*Q[1]-P[1]*Q[0];
 }
 
+def ptinversion(P)
+{
+	if(type(P)==4&&type(P[1])==4){
+		for(R=[];P!=[];P=cdr(P))
+			R=cons(ptinversion(car(P)|option_list=getopt()),R);
+		return reverse(R);
+	}
+	if(type(V=getopt(org))!=0) V=[0,0];
+	if(P==[0,0]) return 0;
+	if(type(P[0])==4){
+		R=P[1];P=P[0];
+	}
+	X=P[0]-V[0];Y=P[1]-V[1];N=X^2+Y^2;
+	if(type(T=getopt(sc))==1||T==0){
+		S=(T<0)?(-T^2):T^2;
+	}else S=-1;
+	if(!R){
+		if(!N) return 0;
+		return [X/N+V[0],S*Y/N+V[1]];
+	}
+	N-=R^2;
+	if(!N){
+		if(X+R!=0) X0=X+R;
+		else X0=X-R;
+		S=[];
+		S=cons(ptinversion([X0,Y]|option_list=getopt()),S);
+		if(Y+R!=0) Y0=Y+R;
+		else Y0=Y-R;
+		return cons(ptinversion([X,Y0]|option_list=getopt()),S);
+	}
+	return [[X/N+V[0],S*Y/N+V[1]],T^2*R/N];
+}
+
 def ptcommon(X,Y)
 {
 	if(length(X)!=2 || length(Y)!=2) return 0;
@@ -17160,55 +18062,71 @@ def xyang(S,P,Q,R)
 	if(type(S)>2) S=dnorm([S,P]);
 	if(type(Prec=getopt(prec))!=1) Prec=0;
 	if(type(Q)>2){
-		if(type(Ar=getopt(ar))!=1) Ar=0;
-		if(R==1||R==-1){				/* 直角 */
-			P1=ptcommon([Q,P],[-S,0]);
-			S*=R;
-			P2=ptcommon([P,P1],[S,@pi/2]);
-			P3=ptcommon([P1,P2],[S,@pi/2]);
-			return xylines([P1,P2,P3]|option_list=Opt);
-		}else if((AR=abs(R))==0||AR==2||AR==3||AR==4||AR>=10){	/* 矢印 */
-			Ang=myarg([Q[0]-P[0],Q[1]-P[1]]);
-			if(R<0) Ang+=3.14159;
-			if(AR>10) X=deval(@pi/180*AR);
+		if(isint(S)&&S<0&&S>-8){
+			if((S=-S)==6||S==7){
+				H=ptcommon([Q,R],[P,0]);
+				if(S==6) return xyang(H,P,0,0|option_list=getopt()); /* 円 */
+				return xylines([P,H]|option_list=getopt()); /* 垂線 */
+			}
+			O=pt5center(P,Q,R);			
+			if(S==2) H=P;	/* 外心 */
 			else{
-				ANG=[0.7854,0.5236,1.0472];
-				X=(AR==0)?1.5708:ANG[AR-2];
+				if(S>2) S++; /* 内心，傍心 */
+				H=ptcommon([P,Q],[O[S],0]);
 			}
-			U=[P[0]+S*dcos(Ang+X),P[1]+S*dsin(Ang+X)];
-			V=[P[0]+S*dcos(Ang-X),P[1]+S*dsin(Ang-X)];	/* 矢先 */
-			L=(X==0)?[U,V]:[U,P,V];
-			if(X&&iand(Ar,2)){
-				L=append([V],L);
-				if((X=ptcommon([P,Q],[U,V]|in=1))!=0) P=X;
+			return xyang(H,O[S],0,0|option_list=getopt());
+		}
+		if(type(Ar=getopt(ar))!=1) Ar=0;
+		if(isint(R)){
+			if(R==1||R==-1){				/* 直角 */
+				P1=ptcommon([Q,P],[-S,0]);
+				S*=R;
+				P2=ptcommon([P,P1],[S,@pi/2]);
+				P3=ptcommon([P1,P2],[S,@pi/2]);
+				return xylines([P1,P2,P3]|option_list=Opt);
+			}else if((AR=abs(R))==0||AR==2||AR==3||AR==4||AR>=10){	/* 矢印 */
+				Ang=myarg([Q[0]-P[0],Q[1]-P[1]]);
+				if(R<0) Ang+=3.14159;
+				if(AR>10) X=deval(@pi/180*AR);
+				else{
+					ANG=[0.7854,0.5236,1.0472];
+					X=(AR==0)?1.5708:ANG[AR-2];
+				}
+				U=[P[0]+S*dcos(Ang+X),P[1]+S*dsin(Ang+X)];
+				V=[P[0]+S*dcos(Ang-X),P[1]+S*dsin(Ang-X)];	/* 矢先 */
+				L=(X==0)?[U,V]:[U,P,V];
+				if(X&&iand(Ar,2)){
+					L=append([V],L);
+					if((X=ptcommon([P,Q],[U,V]|in=1))!=0) P=X;
+				}
+				if(iand(Ar,1))
+					L=append([Q,P,0],L);		/* 心棒 */
+					return xylines(L|option_list=Opt);
+				}else if(AR>4&&AR<9){
+					Ang=myarg([Q[0]-P[0],Q[1]-P[1]]);
+				ANG=[0.7854,0.5236,0.3927,0.2618];
+				X=ANG[AR-5];
+				U=[P[0]+S*dcos(Ang+X),P[1]+S*dsin(Ang+X)];
+				V=[P[0]+S*dcos(Ang-X),P[1]+S*dsin(Ang-X)];
+				W=ptcommon([P,U],[P,Q]|in=-2);
+				W1=[(U[0]+P[0]+W[0])/3,(U[1]+P[1]+W[1])/3];
+				W2=[(V[0]+P[0]+W[0])/3,(V[1]+P[1]+W[1])/3];
+				L=iand(Ar,2)?[V,U,1,W1,P,1,W2,V]:[U,W1,P,1,W2,V];
+				if(iand(Ar,1)){
+					if(iand(Ar,2)) P=ptcommon([P,Q],[U,V]);
+					L=append([Q,P,0],L);
+				};
+				if(type(Sc=getopt(scale))>0){
+					if(type(Sc)==1) Sc=[Sc,Sc];
+					L=ptaffine(diagm(2,Sc),L);
+				}
+				Opt=delopt(Opt,"proc");
+				if(getopt(proc)==1)	return append([2,Opt],L);
+				S=xybezier(L|option_list=Opt);
+				if(getopt(dviout)!=1) return S;
+				dviout(xyproc(S));
+				return 1;
 			}
-			if(iand(Ar,1))
-				L=append([Q,P,0],L);		/* 心棒 */
-			return xylines(L|option_list=Opt);
-		}else if(AR>4&&AR<9){
-			Ang=myarg([Q[0]-P[0],Q[1]-P[1]]);
-			ANG=[0.7854,0.5236,0.3927,0.2618];
-			X=ANG[AR-5];
-			U=[P[0]+S*dcos(Ang+X),P[1]+S*dsin(Ang+X)];
-			V=[P[0]+S*dcos(Ang-X),P[1]+S*dsin(Ang-X)];
-			W=ptcommon([P,U],[P,Q]|in=-2);
-			W1=[(U[0]+P[0]+W[0])/3,(U[1]+P[1]+W[1])/3];
-			W2=[(V[0]+P[0]+W[0])/3,(V[1]+P[1]+W[1])/3];
-			L=iand(Ar,2)?[V,U,1,W1,P,1,W2,V]:[U,W1,P,1,W2,V];
-			if(iand(Ar,1)){
-				if(iand(Ar,2)) P=ptcommon([P,Q],[U,V]);
-				L=append([Q,P,0],L);
-			};
-			if(type(Sc=getopt(scale))>0){
-				if(type(Sc)==1) Sc=[Sc,Sc];
-				L=ptaffine(diagm(2,Sc),L);
-			}
-			Opt=delopt(Opt,"proc");
-			if(getopt(proc)==1)	return append([2,Opt],L);
-			S=xybezier(L|option_list=Opt);
-			if(getopt(dviout)!=1) return S;
-			dviout(xyproc(S));
-			return 1;
 		}
 	}
 	if(type(Q)<3){
@@ -17685,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);
@@ -17761,6 +18687,41 @@ def ptwindow(L,X,Y)
 
 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]);
@@ -17772,6 +18733,24 @@ def pt5center(P,Q,R)
 	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);
 }
 
@@ -19658,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)),","," ");
@@ -19704,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;
@@ -19756,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=[];
@@ -20527,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]);
@@ -20682,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);
@@ -20695,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;
 }