===================================================================
RCS file: /home/cvs/OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v
retrieving revision 1.16
retrieving revision 1.35
diff -u -p -r1.16 -r1.35
--- OpenXM/src/asir-contrib/packages/src/os_muldif.rr	2017/06/08 06:41:51	1.16
+++ OpenXM/src/asir-contrib/packages/src/os_muldif.rr	2018/09/26 04:49:44	1.35
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.15 2017/06/02 05:36:32 takayama Exp $ */
+/* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.34 2018/09/25 00:13:52 takayama Exp $ */
 /* The latest version will be at ftp://akagi.ms.u-tokyo.ac.jp/pub/math/muldif
  scp os_muldif.[dp]* ${USER}@lemon.math.kobe-u.ac.jp:/home/web/OpenXM/Current/doc/other-docs
 */
@@ -6,7 +6,7 @@
 /* #undef USEMODULE */
 
 /*             os_muldif.rr (Library for Risa/Asir) 
- *          Toshio Oshima (Nov. 2007 - June 2017)
+ *          Toshio Oshima (Nov. 2007 - Sep. 2018)
  *
  *   For polynomials and differential operators with coefficients
  *   in rational funtions (See os_muldif.pdf)
@@ -116,6 +116,7 @@ localf lchange$
 localf llsize$
 localf llbase$
 localf lsort$
+localf lpair$
 localf lmax$
 localf lmin$
 localf lgcd$
@@ -181,6 +182,8 @@ localf eta$
 localf jell$
 localf frac$
 localf erfc$
+localf orthpoly$
+localf schurpoly$
 localf fouriers$
 localf todf$
 localf f2df$
@@ -238,6 +241,8 @@ localf expower$
 localf seriesHG$
 localf seriesMc$
 localf seriesTaylor$
+localf mulpolyMod$
+localf taylorODE$
 localf evalred$
 localf toeul$
 localf fromeul$
@@ -272,6 +277,7 @@ localf okuboetos$
 localf heun$
 localf fspt$
 localf abs$
+localf sgn$
 localf calc$
 localf isint$
 localf israt$
@@ -342,6 +348,7 @@ localf getbygrs$
 localf mcop$
 localf shiftop$
 localf conf1sp$
+localf confexp$
 localf pgen$
 localf diagm$
 localf mgen$
@@ -358,6 +365,7 @@ localf fimag$
 localf trig2exp$
 localf intpoly$
 localf integrate$
+localf rungeKutta$
 localf simplog$
 localf fshorter$
 localf isshortneg$
@@ -375,6 +383,7 @@ localf varargs$
 localf ptype$
 localf pfargs$
 localf average$
+localf tobig$
 localf sint$
 localf frac2n$
 localf xyproc$
@@ -399,10 +408,14 @@ localf areabezier$
 localf saveproc$
 localf xygraph$
 localf xy2graph$
+localf addIL$
+localf xy2curve$
+localf xygrid$
 localf xyarrow$
 localf xyarrows$
 localf xyang$
 localf xyoval$
+localf xypoch$
 localf ptcommon$
 localf ptcopy$
 localf ptaffine$
@@ -452,7 +465,7 @@ extern LQS$
 static S_Fc,S_Dc,S_Ic,S_Ec,S_EC,S_Lc$
 static S_FDot$
 extern AMSTeX$
-Muldif.rr="00170607"$
+Muldif.rr="00180925"$
 AMSTeX=1$
 TeXEq=5$
 TeXLim=80$
@@ -551,9 +564,7 @@ def makenewv(L)
 	if((V=getopt(var))<2) V="z_";
 	else if(isvar(V)) V=rtostr(V);
 	if(type(N=getopt(num))!=1) N=0;
-	Var=vars(L);
-	for(Va=Var;Va!=[];Va=cdr(Va))
-		if(vtype(car(Va))==2) Var=append(vars(args(car(Va))),Var);
+	Var=varargs(L|all=2);
 	for(XX=[],I=J=0;;I++){
 		X=strtov(V+rtostr(I));
 		if(findin(X,Var)<0){
@@ -617,6 +628,12 @@ def mycat(L)
 def fcat(S,X)
 {
 	if(type(S)!=7){
+		if(type(DIROUTD)!=7){
+			DIROUTD=str_subst(DIROUT,["%HOME%","%ASIRROOT%","\\"],
+				[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";
@@ -624,16 +641,18 @@ def fcat(S,X)
 		if(S!=0&&access(T)) remove_file(T);
 		S=T;
 	}
-	output(S);
+	R=output(S);
 	print(X);
 	output();
 	if(getopt(exe)==1) shell("\""+S+"\"");
+	return R;
 }
 
 def mycat0(L,T)
 {
 	Opt = getopt(delim);
 	Del = (type(Opt) >= 0)?Opt:"";
+	if(type(L)!=4) L=[L];
 	while(L != []){
 		if(Do==1)
 			print(Del,0);
@@ -729,7 +748,7 @@ def mydiff(P,X)
 		for(;X!=[];X=cdr(X)) P=mydiff(P,car(X));
 		return P;
 	}
-	if(deg(dn(P),X) == 0) 
+	if(ptype(dn(P),X)<2) 
 		return red(diff(nm(P),X)/dn(P));
 	return red(diff(P,X));
 }
@@ -1245,22 +1264,134 @@ def mulseries(V1,V2)
 
 def scale(L)
 {
-	T=0;LS=1;
-	if(type(L)!=4){
+	T=F=0;LS=1;
+	Pr=getopt(prec);
+	Inv=getopt(inv);
+	Log10=dlog(10);
+	if(type(L)==7){
+		V=findin(L,["CI","DI","CIF","CIF'","DIF","DIF'","SI","TI1","TI2","STI"]);
+		if(V>=0){
+			L=["C","D","CF","CF'","DF","DF'","S","T1","T2","ST"];
+			Inv=1;L=L[V];
+		}
+		V=findin(L,["C","A","K","CF","CF'","S","T1","T2","ST","LL0","LL1","LL2","LL3","LL00",
+			"LL01","LL02","LL03"])+1;
+		if(V==0) V=findin(L,["D","B","K","DF","DF'"])+1;
+		if(V>0) L=V;
+	}
+	if(type(OL=L)!=4){
 		if(L==2){
-			L=[[[1,2,1/20],[2,5,1/10],[5,10,1/5], [10,20,1/2],[20,50,1],[50,100,2]],
+			L=(Pr==0)?
+			  [[[1,2,1/20],[2,5,1/10],[5,10,1/5], [10,20,1/2],[20,50,1],[50,100,2]],
     		   [[1,2,1/10],[2,5,1/2], [10,20,1],[20,50,5]],
-    		   [[1,2,1/2],[2,10,1], [10,20,5],[20,100,10]]];
-			LS=2;
+    		   [[1,2,1/2],[2,10,1], [10,20,5],[20,100,10]]]:
+			   [[[1,2,1/50],[2,5,1/20],[5,10,1/10], [10,20,1/5],[20,50,1/2],[50,100,1]],
+				[[1,5,1/10],[5,10,1/2], [10,20,1],[50,100,5]],
+				[[1,5,1/2],[5,10,1], [10,50,5],[50,100,10]]];
+			LS=2;M2=[[1,10,1],[10,100,10]];
 		}else if(L==3){
-			L=[[[1,2,1/20],[2,5,1/10],[5,10,1/5], [10,20,1/2],[20,50,1],[50,100,2],
-                [100,200,5],[200,500,10],[500,1000,20]],
+			L=(Pr==0)?
+			  [[[1,2,1/20],[2,5,1/10],[5,10,1/5], [10,20,1/2],[20,50,1],[50,100,2],
+                  [100,200,5],[200,500,10],[500,1000,20]],
     			[[1,2,1/10],[2,5,1/2], [10,20,1],[20,50,5], [100,200,10],[200,500,50]],
-    			[[1,2,1/2],[2,10,1], [10,20,5],[20,100,10], [100,200,50],[200,1000,100]]];
-			LS=3;
+    			[[1,2,1/2],[2,10,1], [10,20,5],[20,100,10], [100,200,50],[200,1000,100]]]:
+			  [[[1,2,1/50],[2,5,1/20],[5,10,1/10],[10,20,1/5],[20,50,1/2],[50,100,1], 
+			      [100,200,2],[200,500,5],[500,1000,10]],
+				[[1,5,1/10],[5,10,1/2], [10,50,1],[50,100,5], [100,500,10],[500,1000,50]],
+				[[1,5,1/2],[5,10,1],[10,50,5],[50,100,10], [100,500,50],[500,1000,100]]];
+			LS=3;M2=[[1,5,1],[10,50,10],[100,500,100],[500,1000,500]];
+		}else if(L>9&&L<18){
+			if(L<18){	/* LL0 - LL3, LL00 - LL03 */
+				if(L==10){
+					L=[ [[1.001,1.002,0.00001],[1.002,1.005,0.00002],[1.005,1.0105,0.00005]],
+						[[1.001,1.002,0.00005],[1.002,1.005,0.0001], [1.005,1.0105,0.0001]],
+						[[1.001,1.002,0.0001],[1.002,1.005,0.0005], [1.005,1.0105,0.0005]]];
+					M2=[1.001,1.0015,1.002,1.003,1.004,1.005,1.006,1.007,1.008,1.009,1.01];
+				}
+				if(L==11){
+					L=[ [[1.01,1.02,0.0001],[1.02,1.05,0.0002],[1.05,1.105,0.0005]],
+						[[1.01,1.02,0.0005],[1.02,1.05,0.001], [1.05,1.105,0.001]],
+						[[1.01,1.02,0.001],[1.02,1.05,0.005], [1.05,1.105,0.005]]];
+					M2=[1.01,1.015,1.02,1.03,1.04,1.05,1.06,1.07,1.08,1.09,1.10];
+				}else if(L==12){
+					L=[ [[1.105,1.2,0.001],[1.2,1.4,0.002],[1.4,1.8,0.005],[1.8,2.5,0.01],
+					      [2.5,2.72,0.02]],
+						[[1.105,1.2,0.005],[1.2,1.4,0.01],[1.4,1.8,0.01],[1.8,2.5,0.05],
+					      [2.5,2.72,0.1]],
+						[[1.105,1.2,0.01],[1.2,1.4,0.05],[1.4,1.8,0.05],[1.8,2.5,0.1],
+					      [2.5,2.72,0.1]]];
+					M2=[1.11,1.15,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.5];
+				}else if(L==13){
+					L=[ [[2.72,4,0.02],[4,6,0.05],[6,10,0.1],[10,15,0.2],[15,30,0.5],[30,50,1],
+						 [50,100,2],[100,200,5],[200,400,10],[400,500,20],[500,1000,50],
+						 [1000,2000,100],[2000,5000,200],[5000,10000,500],[10000,22000,1000]],
+						[[2.7,4,0.1],[4,6,0.1],[6,10,0.5],[10,15,1],[15,30,1],[30,50,5],
+						 [50,100,10],[100,200,10],[200,400,50],[400,500,100],[500,1000,100],
+						 [1000,2000,500],[2000,5000,1000],[5000,10000,1000],[10000,22000,5000]],
+						[[3,4,0.5],[4,6,0.5],[6,10,1],[10,15,5],[15,30,5],[30,50,10],
+						 [50,100,50],[100,200,50],[200,400,100],[400,500,100],[500,1000,500],
+						 [1000,2000,1000],[2000,5000,3000],[5000,10000,5000],[10000,22000,10000]]];
+					M2=[3,4,5,6,7,8,9,10,15,20,30,40,50,100,200,500,1000,2000,5000,10000,20000];
+				}else if(L==14){
+					L=[ [[0.998,0.999,0.00001],[0.995,0.998,0.00002],[0.99,0.995,0.00005]],
+						[[0.998,0.999,0.00005],[0.995,0.998,0.0001],[0.99,0.995,0.0001]],
+						[[0.998,0.999,0.0001],[0.995,0.998,0.0005],[0.99,0.995,0.0005]]];
+					M2=[0.999,0.9985,0.998,0.997,0.996,0.995,0.994,0.993,0.992,0.991,0.99];
+				}else if(L==15){
+					L=[ [[0.98,0.9901,0.0001],[0.95,0.98,0.0002],[0.905,0.95,0.0005]],
+						[[0.98,0.99,0.0005],[0.95,0.98,0.001], [0.905,0.95,0.001]],
+						[[0.98,0.99,0.001],[0.95,0.98,0.005], [0.91,0.95,0.005]]];
+					M2=[0.99,0.985,0.98,0.97,0.96,0.95,0.94,0.93,0.92,0.91];
+				}else if(L==16){
+					L=[ [[0.8,0.906,0.001],[0.6,0.8,0.002],[0.37,0.6,0.005]],
+						[[0.8,0.906,0.005],[0.6,0.8,0.01],[0.37,0.6,0.01]],
+						[[0.8,0.9,0.01],[0.6,0.8,0.05],[0.4,0.6,0.05]]];
+					M2=[0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.55,0.5,0.45,0.4];
+				}else{
+					L=[ [[0.05,0.37,0.002],[0.02,0.05,0.001],[0.01,0.02,0.0005],
+						 [0.005,0.01,0.0002],[0.001,0.005,0.0001],
+						 [0.0005,0.001,0.00002],[0.0001,0.0005,0.00001],[0.00005,0.0001,0.000002]],
+						[[0.05,0.37,0.01],[0.02,0.05,0.002],[0.01,0.02,0.001],
+						 [0.005,0.01,0.001],[0.001,0.005,0.0002],
+						 [0.0005,0.001,0.0001],[0.0001,0.0005,0.00002],[0.00005,0.0001,0.00001]],
+						[[0.05,0.37,0.05],[0.02,0.05,0.01],[0.01,0.02,0.005],
+						 [0.005,0.01,0.005],[0.002,0.005,0.001],
+						 [0.0005,0.001,0.0005],[0.0001,0.0005,0.0001],[0.00005,0.0001,0.00005]]];
+					M2=[0.3,0.2,0.1,0.05,0.03,0.02,0.01,0.005,0.002,0.001,0.0005,0.0002,0.0001];
+				}
+			}
 		}else{
-			L=[[[1,2,1/50],[5,10,1/2],[5,10,1/10]],[[1,5,1/10],[5,10,1/2]],
-				[[1,5,1/2],[5,10,1]]];
+			if(L==6){	/* S */
+				L=[ [[6-3/12,15,1/12],[15,30,1/6],[30,50,1/3],[50,70,1/2],[70,80,1],[80,90,5]],
+					[[6-1/6,15,1/6],[15,30,1/2],[30,70,1],[70,80,5],[80,90,10]],
+					[[6,15,1/2],[15,30,1],[30,70,5],[70,90,10]] ];
+				M2=[6,7,8,9,10,15,20,30,40,50,60,70,90];
+			}else if(L==7){ /* T1 */
+				F=log(tan(x*3.1416/180))/Log10+1;
+				L=[ [[6-1/3,15,1/12],[15,45,1/6]],
+					[[6-1/3,15,1/6],[15,45,1/2]],
+					[[6,45,1]] ];
+				M2=[6,7,8,9,10,15,20,30,40,45];
+			}else if(L==8){	/* T2 */
+				L=[ [[45,75,1/6],[75,84+1/6,1/12]],
+					[[45,75,1],[75,84+1/6,1/6]],
+					[[45,84,1]] ];
+				M2=[45,50,60,70,75,80,81,82,83,84];
+			}else if(L==9){	/* ST */
+				L=[ [[35/60,1,1/120],[1,2,1/60],[2,5+9/12,1/30]],
+					[[35/60,1,1/60],[1,2,1/6],[2,5+9/12,1/6]],
+					[[40/60,1,1/6],[1,2,1/2],[2,5+9/12,1]] ];
+				M2=[1,2,3,4,5];
+			}else{
+				M2=(L==4||L==5)?[[1,2,1/2],[2,9,1]]:[[1,2,1/2],[2,10,1]];
+				L=(Pr==0)?
+			  	[ [[1,2,1/50],[2,5,1/20],[5,10,1/10]],
+				  [[1,5,1/10],[5,10,1/2]],
+				  [[1,5,1/2],[5,10,1]] ]:
+				[[[1,2,1/100],[2,5,1/50],[5,10,1/20]],
+			      [[1,2,1/20],[2,10,1/10]],
+			      [[1,2,1/10],[2,10,1/2]] ];
+			}
 		}
 	}else if(type(L[0])!=4){
 		L=[L];
@@ -1285,41 +1416,99 @@ def scale(L)
 	SC=getopt(scale);
 	if(type(SC)==4){
 		S0=SC[0];S1=SC[1];
-	}else if(type(S)==1){
-		S0=S;S1=0;
+	}else if(type(SC)==1){
+		S0=SC;S1=0;
 	}else return T;
 	if(type(D=getopt(shift))==4){
 		D0=D[0];D1=D[1];
+	}else if(type(D)<2&&type(D)>=0){
+		D0=0;D1=D;
+	};
+	if(Inv==1){
+		D0+=S0;S0=-S0;
 	}
-	if(type(F=getopt(f))>1) F=f2df(F);
-	else F=0;
-	for(M=[],I=length(T);T!=[];T=cdr(T),I--){
+	if(type(TF=getopt(f))>1) F=TF;
+	if(F) F=f2df(F);
+	if(type(I=getopt(ol))==1&&OL>3) OL=I;
+	for(M=M0=[],I=length(T);T!=[];T=cdr(T),I--){
 		for(S=car(T);S!=[];S=cdr(S)){
-			V=((!F)?dlog(car(S))/dlog(10)/LS:myfdeval(F,car(S)))*S0;
+			VS=car(S);
+			if(F) V=myfdeval(F,car(S));
+			else if(OL==4) V=frac(dlog(VS)/Log10+0.5);
+			else if(OL==5) V=frac(dlog(VS*3.1416)/Log10);
+			else if(OL>5&&OL<10){
+				VS=VS*3.1416/180;
+				if(OL==6) V=dlog(dsin(VS))/Log10+1;
+				else if(OL==9) V=dlog(VS)/Log10+2;
+				else V=dlog(dtan(VS))/Log10+8-OL; 
+			}
+			else if(OL>9&&OL<14) V=dlog(dlog(VS))/Log10+13-OL;
+			else if(OL>13&&OL<18) V=dlog(-dlog(VS))/Log10+17-OL;
+			else V=dlog(VS)/Log10/LS;
+			V*=S0;
 			if(S1!=0){
 				M=cons([V+D0,D1],M);
-				M=cons([V+D0,I*((length(SC)>2)?SC[I]:S1)+D1],M);
+				M=cons([V+D0,((length(SC)>2)?SC[I]:(I*S1))+D1],M);
 				M=cons(0,M);
-			}else M=cons(V+D0,M);
+			}else M0=cons(V+D0,M0);
 		}
+		if(S1==0) M=cons(reverse(M0),M);
 	}
 	if(S1!=0) M=cdr(M);
+	if(S1==0||getopt(TeX)!=1) return M;
 	M=reverse(M);
-	if(getopt(TeX)==1){
-		if(type(U=getopt(line))==4)
-			M=cons([U[0]+D0,D1],cons([U[1]+D0,D1],cons(0,M)));
-		S=xylines(M);
-		if(type(Mes=getopt(mes))==4){
-			S3=car(Mes);
-			V=car(scale(cdr(Mes)));
-			if(!F) Mes=scale(cdr(Mes)|scale=[S0/LS,0],shift=[D0,D1]);
-			else Mes=scale(cdr(Mes)|f=F,scale=[S0,0],shift=[D0,D1]);
-			for(M=Mes;M!=[];M=cdr(M),V=cdr(V))
-				S+=xyput([car(M),S3,deval(car(V))]);
+	if(type(U=getopt(line))==4){
+		if(Inv==1) U=[U[0]+S0,U[1]+S0];
+		M=cons([U[0]+D0,D1],cons([U[1]+D0,D1],cons(0,M)));
+	}
+	if((VT=getopt(vert))==1){
+		for(N=[];M!=[];M=cdr(M)){
+			if(type(TM=car(M))==4) N=cons([TM[1],TM[0]],N);
+			else N=cons(TM,N);
 		}
-		return S;
+		M=reverse(N);
 	}
-	return M;
+	if(type(Col=getopt(col))<1) S=xylines(M);
+	else S=xylines(M|opt=Col);
+	if(type(Mes=getopt(mes))==4){
+		if(length(Mes)==1&&type(M2)==4) Mes=cons(car(Mes),M2);
+		S3=car(Mes);
+		if(type(S3)==4){
+			Col=S3[1];
+			S3=car(S3);
+		}else Col=0;
+		V=car(scale(cdr(Mes)));
+		if(!F) Mes=scale(cdr(Mes)|scale=[S0/LS,0],shift=[D0,D1],ol=OL);
+		else Mes=scale(cdr(Mes)|f=F,scale=[S0,0],shift=[D0,D1]);
+		for(M=car(Mes);M!=[];M=cdr(M),V=cdr(V)){
+			TV=deval(car(V));
+			if(Col!=0) TV=[Col,TV];
+			S+=(VT==1)?xyput([S3+D1,car(M),TV]):xyput([car(M),S3+D1,TV]);
+		}
+	}
+	if(type(Mes=getopt(mes2))==4){
+		if(type(car(Mes))!=4) Mes=[Mes];
+		for(;Mes!=[];Mes=cdr(Mes)){
+			TM=car(Mes);
+			if(!F) V=scale([car(TM)]|scale=[S0/LS,0],shift=[D0,D1],ol=OL);
+			else V=scale([car(TM)]|f=F,scale=[S0,0],shift=[D0,D1]);
+			V=car(car(V));
+			TM=cdr(TM);
+			if(type(Col=car(TM))==4){
+				C0=Col[0];C1=Col[1];
+				if(length(Col)==3){
+					S+=(VT==1)?xyline([D1+C0,V],[D1+C1,V]|opt=Col[2])
+						:xyline([V,D1+C0],[V,D1+C1]|opt=Col[2]);
+				}else S+=(VT==1)?xyline([D1+C0,V],[D1+C1,V]):xyline([V,D1+C0],[V,D1+C1]);
+			}
+			if(type(TM[1]<2)){
+				TM=cdr(TM);
+				S3=car(TM);
+			}
+			S+=(VT==1)?xyput([S3+D1,V,TM[1]]):xyput([V,S3+D1,TM[1]]);
+		}
+	}
+	return S;
 }
 
 def pluspower(P,V,N,M)
@@ -3247,7 +3436,7 @@ def lsort(L1,L2,T)
 			}else if(L2=="transpose") return mtranspose(L1);
 			else if(L2=="subst"||L2=="adjust"){
 				Null=(!K)?"":car(K);
-				if(L2=="adjust) C1=[];
+				if(L2=="adjust") C1=[];
 				R=lv2m(L1|null="");
 				for(;C1!=[];C1=cdr(C1)) R[car(C1)[0]][car(C1)[1]]=car(C1)[2]; 
 				return m2ll(R);
@@ -3492,6 +3681,20 @@ def msort(L,S)
 	return qsort(L,os_md.mqsub);
 }
 
+def lpair(A,B)
+{
+	if(B==0){
+		for(S=T=[];A!=[];A=cdr(A)){
+			S=cons(car(A)[0],S);T=cons(car(A)[1],T);
+		}
+		return [reverse(S),reverse(T)];
+	}else{
+		for(R=[];A!=[];A=cdr(A),B=cdr(B))
+			R=cons([car(A),car(B)],R);
+		return reverse(R);
+	}
+}
+
 def lmax(L)
 {
 	if(type(L)==4){
@@ -3625,16 +3828,34 @@ def lnsol(VV,L)
 
 def ladd(X,Y,M)
 {
-	if(type(X)==4) X=ltov(X);
 	if(type(Y)==4) Y=ltov(Y);
+	if(type(X)==4) X=ltov(X);
 	return vtol(X+M*Y);
 }
 
 def mrot(X)
 {
+	if(type(X)==4){
+		if(getopt(deg)==1)
+			X=[deval(@pi*X[0]/180),deval(@pi*X[1]/180),deval(@pi*X[2]/180)];
+		if(getopt(conj)==1)
+			return mrot([-X[2],-X[1],0])*mrot([X[0],X[1],X[2]]);
+		if(X[1]==0){
+			X=[X[0]+X[2],0,0];
+			if(X[0]==0) return diagm(3,[1]);
+		}
+		if(X[0]!=0){
+			M=mat([dcos(X[0]),-dsin(X[0]),0],[dsin(X[0]),dcos(X[0]),0],[0,0,1]);
+			if(X[1]==0) return M;
+		}
+		N=mat([dcos(X[1]),0,-dsin(X[1])],[0,1,0],[dsin(X[1]),0,dcos(X[1])]);
+		if(X[0]!=0) N=M*N;
+		if(X[2]==0) return N;
+		return N*mrot([X[2],0,0]);
+	}
 	if(getopt(deg)==1) X=@pi*X/180;
 	X=deval(X);
-	return mat([dcos(X),dsin(X)],[-dsin(X),dcos(X)]);
+	return mat([dcos(X),-dsin(X)],[dsin(X),dcos(X)]);
 }
 
 def m2v(M)
@@ -4728,6 +4949,7 @@ def myswap(P,L)
 def mysubst(P,L)
 {
 	if(P==0) return 0;
+	if(getopt(lpair)==1||(type(L[0])==4&&length(L[0])>2)) L=lpair(L[0],L[1]);
 	Inv=getopt(inv);
 	if(type(L[0]) == 4){
 		while((L0 = car(L))!=[]){
@@ -4923,11 +5145,11 @@ def muldo(P,Q,L)
 def jacobian(F,X)
 {
 	F=ltov(F);X=ltov(X);
-	N=length(F);
-	M=newmat(N,N);
+	N=length(F);L=length(X);
+	M=newmat(N,L);
 	for(I=0;I<N;I++)
-		for(J=0;J<N;J++) M[I][J]=red(diff(F[I],X[J]));
-	if(getopt(mat)==1) return M;
+		for(J=0;J<L;J++) M[I][J]=red(diff(F[I],X[J]));
+	if(N!=L||getopt(mat)==1) return M;
 	return mydet(M);
 }
 
@@ -5780,7 +6002,7 @@ def pol2sft(F,A)
 
 def binom(P,N)
 {
-	if(type(N)!=1 || N<0) return 1;
+	if(type(N)!=1 || N<=0) return 1;
 	for(S=1;N>0;N--,P-=1)	S*=P/N;
 	return red(S);
 }
@@ -5796,6 +6018,7 @@ def expower(P,R,N)
 
 def seriesHG(A,B,X,N)
 {
+	if(N==0) return 1;
 	if(type(N)!=1 || N<0) return 0;
 	if(type(X)<4){
 		for(K=0,S=S0=1;K<N;K++){
@@ -5830,14 +6053,31 @@ def evalred(F)
 		 Opt=[];
 	}else if(length(Opt)==2 && type(Opt[0])!=4)	Opt=[Opt];
 	for(;;){
-		G=mysubst(F,[[sin(0),0],[tan(0),0],[asin(0),0],[atan(0),0],[sinh(0),0],[tanh(0),0],
-			[log(1),0],[cos(0),1],[cosh(0),1],[exp(0),1]]);
+		G=mysubst(F,[[tan(0),0],[asin(0),0],[atan(0),0],[sinh(0),0],[tanh(0),0],
+			[log(1),0],[cosh(0),1],[exp(0),1]]);
 		for(Rep=Opt; Rep!=[]; Rep=cdr(Rep))
 			G=subst(G,car(Rep)[0],car(Rep)[1]);
 		Var=vars(G);
 		for(V=Var; V!=[]; V=cdr(V)){
-			if(functor(car(V))!=pow || (P=args(car(V))[0])!=1) continue;
-			G=subst(G,car(V),1);
+			if(!(VV=args(CV=car(V)))) continue;
+			if((functor(CV)==sin||functor(CV)==cos)){
+				P=2*red(VV[0]/@pi);
+				if(functor(CV)==sin) P=1-P;
+				if(isint(P)){
+					if(iand(P,1)) G=subst(G,CV,0);
+					else if(!iand(P,3)) G=subst(G,CV,1);
+					else G=subst(G,CV,-1);
+					continue;
+				}
+				if(isint(P*=3/2)){
+					if(iand(P,3)==1) G=subst(G,CV,1/2);
+					else G=subst(G,CV,-1/2);
+				}
+			}
+			for(;VV!=[];VV=cdr(VV))
+				if(car(VV)!=(TV=evalred(car(VV)))) G=subst(G,car(VV),TV);
+			if(functor(CV)!=pow || (args(CV)[0])!=1) continue;
+			G=subst(G,CV,1);
 		}
 		if(G==F) return F;
 		F=G;
@@ -5912,6 +6152,110 @@ def seriesTaylor(F,N,V)
 	return F;
 }
 
+def mulpolyMod(P,Q,X,N)
+{
+	Red=(type(P)>2||type(Q)>2)?1:0;
+	for(I=R=0;I<=N;I++){
+		P0=mycoef(P,I,X);
+		for(J=0;J<=N-I;J++){
+			R+=P0*mycoef(Q,J,X)*X^(I+J);
+			if(Red) R=red(R);
+		}
+	}
+	return R;
+}
+
+def taylorODE(D){
+	Dif=(getopt(dif)==1)?1:0;
+	if(D==0) return Dif?f:f_00;
+	if(type(T=getopt(runge))!=1||ntype(T)!=0) T=0;
+	if(type(F=getopt(f))!=7&&type(F)<2) F="f_";
+	if(type(D)!=1||ntype(D)!=0||D<0||D>30) return 0;
+	if(type(H=getopt(taylor))==4&&length(H)==2){
+		if(type(Lim=getopt(lim))==2) DD=D;
+		else if(type(Lim)==4){
+			DD=Lim[1];Lim=Lim[0];
+		}else Lim=0;
+		for(R=I=0;I<=D;I++){
+			if(I){
+				if(Lim) H0=mulpolyMod(H0,H[0],Lim,DD);
+				else H0*=H[0];
+			}else  H0=1;
+			if(type(F)!=7) G=I?mydiff(G,x):F;
+			for(J=0;J<=D-I;J++){
+				if(J){
+					if(Lim) H1=mulpolyMod(H1,H[1],Lim,DD);
+					else H1*=H[1];
+				}else H1=H0;
+				if(type(F)==7) G=makev([F,I,J]);
+				else if(J) G=mydiff(G,y);
+				R+=G*H1/fac(I)/fac(J);
+			}
+		}
+		if(Lim) R=os_md.polcut(R,DD,Lim);
+		return R;
+	}else{
+		if(type(H=getopt(series))>=0||getopt(list)==1){
+			if(type(F)!=7){
+				for(PP=[F],I=1;I<D;I++)
+					PP=cons(mydiff(car(PP),x)+mydiff(car(PP),y)*F,PP);
+				if(type(H)<0) return PP;
+				for(R=0,DD=D;DD>=1;DD--,PP=cdr(PP)) R+=car(PP)*H^DD/fac(DD);
+				return red(R);
+			}
+			if(type(H)>=0) D--;
+			PP=taylorODE(D-1|list=1);
+			if(type(PP)!=4) PP=[PP];
+			P=car(PP);
+		}else P=taylorODE(D-1);
+		for(R=I=0;I<D;I++){
+			for(J=0;J<D-I;J++){
+				Q=diff(P,makev([F,I,J]));
+				if(Q!=0) R+=Q*(f_00*makev([F,I,J+1])+makev([F,I+1,J]));
+			}
+		}
+		if(getopt(list)==1){
+			R=cons(R,PP);
+			if(Dif!=1) return R;
+		}else if(type(H)>=0){
+			R=y+R*H^(D+1)/fac(D+1);
+			for(DD=D;DD>0;PP=cdr(PP),DD--) R+=car(PP)*H^(DD)/fac(DD);
+			if(T){
+				if(T<0){
+					Dif=0;TT=-T;
+				}else TT=T;
+				K=newvect(TT);K[0]=Dif?f:f_00;
+				if(getopt(c1)==1) K[0]=taylorODE(D|taylor=[c_1*H,0]);
+				for(I=1;I<TT;I++){
+					for(S=J=0;J<I;J++) S+=makev(["a_",I+1,J+1])*K[J];
+					K[I]=taylorODE(D|taylor=[makev(["c_",I+1])*H,S*H],lim=[H,D]);
+				}
+				for(S=I=0;I<TT;I++) S+=makev(["b_",I+1])*K[I];
+				S=S*H+y;
+				R=S-R;
+				if(T<0){
+					for(V=[H],I=0;I<=D;I++)
+						for(J=0;J<=D-I;J++) V=cons(makev([F,I,J]),V);
+					return os_md.ptol(R,reverse(V)|opt=0);
+				}
+			}else T=0;
+		}
+	}
+	if(Dif){
+		for(I=0;I<=D;I++){
+			for(J=0;J<=D;J++){
+				if(I==0&&J==0){
+					R=subst(R,f_00,f);
+					continue;
+				}
+				V=makev([F,str_times("x",I),str_times("y",J)]);
+				R=subst(R,makev([F,I,J]),V);
+			}
+		}
+	}
+	return R;
+}
+
 def toeul(F,L,V)
 {
 	L = vweyl(L); 
@@ -7399,6 +7743,22 @@ def abs(X)
 	return X;
 }
 
+def sgn(X)
+{
+	if(X==0) return 0;
+	if(type(X)==1){
+		return (X>0)?1:-1;
+	}
+	if(type(X)==5) X=vtol(X);
+	if(type(X)==4){
+		for(W=0,Y=X;Y!=[];Y=cdr(Y))
+			for(Z=cdr(Y);Z!=[];Z=cdr(Z))
+				if(car(Y)>car(Z)) W++;
+		if(getopt(val)==1) return W;
+		return (iand(W,1))?-1:1;
+	}
+}
+
 def calc(X,L)
 {
 	if(type(X)<4||type(X)==7){
@@ -7432,6 +7792,12 @@ def calc(X,L)
 	return X;
 }
 
+def tobig(X)
+{
+	if((type(X)==1 && ntype(X)==3)||type(X)>3) return X;
+	return eval(X*exp(0));
+}
+
 def isint(X)
 {
 	if(X==0||(type(X)==1 && ntype(X)==0 && dn(X)==1)) return 1;
@@ -7677,6 +8043,8 @@ def spgen(MO)
 	if(F!=1&&F!=-1) F=0;
 	if(type(LP)==4){
 		L0=LP[0]; L1=LP[1];
+	}else if(type(LP)==1){
+		L0=L1=LP;
 	}else{
 		L0=0; L1=MO+1;
 	}
@@ -7701,6 +8069,7 @@ def spgen(MO)
 					B=cons(sub_str(S,P0,P1-1),B);
 				}
 			}
+			close_file(Id);
 		}else{
 			MO/=2;
 			if(L1<=1) L1=MO+4;
@@ -7966,7 +8335,7 @@ def chkspt(M)
 	Opt= getopt(opt);
 	Mat= getopt(mat);
 	if(type(M)==7) M=s2sp(M);
-	if(type(Opt) >= 0){
+	if(type(Opt) >= 0&&Opt!="idx"){
 		if(type(Opt) == 7)
 			Opt = findin(Opt, ["sp","basic","construct","strip","short","long","sort","root"]);
 		if(Opt < 0){
@@ -7975,7 +8344,6 @@ def chkspt(M)
 		}
 		return fspt(M,Opt);
 	}
-	MR = fspt(M,1);
 	P  = length(M);
 	OD = -1;
 	XM = newvect(P);
@@ -7998,8 +8366,8 @@ def chkspt(M)
 		if(OD < 0)
 			OD = SM;
 		else if(OD != SM){
-			print("irregal partitions");
-			return 0;
+			if(getopt(dumb)!=1) print("irregal partitions");
+			return -1;
 		}
 		XM[I] = JM;
 	}
@@ -8016,13 +8384,14 @@ def chkspt(M)
 		SM += MV;
 	}
 	SM -= (P-2)*OD;
+	if(Opt=="idx") return SSM;
 	if(SM > SMM && SM != 2*OD){
-		print("not realizable");
-		return -1;
+		if(getopt(dumb)!=1) print("not realizable");
+		return 0;
 	}
 	if(JM==1 && Mat!=1)
 		Fu -= OD - SSM/2;
-	return [P, OD, SSM, Fu, SM, XM, MR];
+	return [P, OD, SSM, Fu, SM, XM, fspt(M,1)];
 }
 
 def cterm(P)
@@ -8137,7 +8506,7 @@ def redgrs(M)
 				}
 				L = cons([VM,EV], L);
 /*
-				if(R[2] >= 2){ */ /* digid */
+				if(R[2] >= 2){ */ /* rigid */
 /*          P = dx^(R[1]);
 				} */
 			}
@@ -8164,30 +8533,51 @@ def mcgrs(G, R)
 {
 	NP = length(G);
 	Mat = (getopt(mat)==1)?0:1;
+	if(Mat==1 && type(SM=getopt(slm))==4){
+		SM0=SM[0];SM1=anal2sp(SM[1],["*",-1]);
+		if(findin(0,SM0)>=0){
+			for(SM=[],I=length(G)-1;I>0;I--)
+				if(findin(I,SM0)<0) SM=cons(I,SM);
+			SM=[SM,SM1];
+			G=mcgrs(G,R|slm=SM);
+			return [G[0],anal2sp(G[1],["*",-1])];
+		}
+	}else SM0=0;
 	for(R = reverse(R) ; R != []; R = cdr(R)){
 		GN = [];
 		L = length(G)-1;
 		RT = car(R);
 		if(type(RT) == 4){
-			RT = reverse(RT); S = 0;
+			if(length(RT)==L&&RT[0]!=0){
+				R=cons(cdr(RT),R);
+				R=cons(RT[0],R);
+				continue;
+			}		/* addition */
+			RT = reverse(RT); S = ADS = 0;
 			for(G = reverse(G); G != []; G = cdr(G), L--){
 				AD = car(RT); RT = cdr(RT);
-				if(L > 0)
+				if(L > 0){
 					S += AD;
-				else
+					if(SM && findin(L,SM0)>=0) ADS+=AD;
+				}else
 					AD = -S;
 				for(GTN = [], GT = reverse(car(G)); GT != []; GT = cdr(GT))
 					GTN = cons([car(GT)[0],car(GT)[1]+AD], GTN);
 				GN = cons(GTN, GN);
 			}
 			G = GN;
+			if(SM0){
+				for(ST=reverse(SM1),SM1=[]; ST!=[]; ST=cdr(ST))
+					SM1 = cons([car(ST)[0],car(ST)[1]+ADS], SM1);
+			}
 			continue;
 		}
-		VP = newvect(L+1); GV = ltov(G);
+		if(RT==0) continue;
+		VP = newvect(L+1); GV = ltov(G);	/* middle convolution */
 		for(I = S = OD = 0; I <= L; I++){
 			RTT = (I==0)?(Mat-RT):0;
 			VP[I] = -1;
-			for(J = M = 0, GT = GV[I]; GT != []; GT = cdr(GT), J++){
+			for(J = M = K = 0, GT = GV[I]; GT != []; GT = cdr(GT), J++){
 				if(I == 0)
 					OD += car(GT)[0];
 				if(car(GT)[1] == RTT && car(GT)[0] > M){
@@ -8195,29 +8585,45 @@ def mcgrs(G, R)
 					VP[I] = J;
 				}
 			}
-			S -= (L-1)*OD;
-			for(GN = [] ; L >= 0; L--){
-				GT = GV[L];
-				RTT = (L==0)?(-RT):RT;
-				FTN = (VP[L] >= 0 || S == 0)?[]:[-S,(L==0)?(Mat-RT):0];
-				for(J = 0; GT != []; GT = cdr(GT), J++){
-					if(J != VP[L]){
-						GTN = cons([car(GT)[0],car(GT)[1]+RTT], GTN);
-						continue;
-					}
-					K = car(GT)[0] - S;
-					if(K < 0){
-						print("Not realizable");
-						return;
-					}
-					GTN = cons([K,(L==0)?(Mat-RT):0], GTN);
+		}
+		S -= (L-1)*OD;
+		for(GN = []; L >= 0; L--){
+			GT = GV[L];
+			RTT = (L==0)?(-RT):RT;
+			GTN = [(VP[L] >= 0 || S == 0)?[]:[-S,(L==0)?(Mat-RT):0]];
+			for(J = 0; GT != []; GT = cdr(GT), J++){
+				if(J != VP[L]){
+					GTN = cons([car(GT)[0],car(GT)[1]+RTT], GTN);
+					continue;
 				}
-				GN = cons(reverse(GTN), GN);
+				K = car(GT)[0] - S;
+				if(K < 0){
+					print("Not realizable");
+					return;
+				}
+				GTN = cons([K,(L==0)?(Mat-RT):0], GTN);
 			}
+			if(VP[L]<0) GTN=cons([-S,(L==0)?(Mat-RT):0],GTN);
+			GN = cons(reverse(GTN), GN);
 		}
+		if(SM0){
+			for(M=0,L=length(G)-1;L>0;L--)
+				if(findin(L,SM0)>=0&&VP[L]>=0) M+=GV[L][VP[L]][0];
+			Mx=sp2anal(SM1,["max",1,0]);
+			for(SM2=[],J=0;SM1!=[];J++,SM1=cdr(SM1)){
+				if(J!=Mx[0]) SM2=cons([car(SM1)[0],car(SM1)[1]+RT],SM2);
+				else if((V=car(SM1)[0]-M)!=0) SM2=cons([V,RT],SM2);
+			}
+			Mx=sp2anal(SM1,["max",1,0]);
+			for(J=0;SM2!=[];J++,SM2=cdr(SM2)){
+				if(J!=Mx[0]) SM1=cons(car(SM2),SM1);
+				else if((V=car(SM1)[0]-S+M)!=0) SM1=cons([V,0],SM2);
+			}
+			if(Mx[0]<0) SM1=cons([-S,0],SM1);
+		}
 		G = cutgrs(GN);
 	}
-	return G;
+	return SM0?[G[0],SM1]:G;
 }
 
 /*
@@ -8231,6 +8637,8 @@ def mcgrs(G, R)
   F=["put",F,V] : 
   F=["get1",F,V] : 
   F=["put1",F,V] : 
+  F=["max"] : 
+  F=["max",F.V] : 
   F=["put1"] : 
   F=["val",F];
   F=["swap"];
@@ -8275,6 +8683,18 @@ def anal2sp(R,F)
 		return G;
 	}
 	if(F[0]=="add") return append(R,F[1]);
+	if(F[0]=="max"){
+		if(length(F)==3) C=1;
+		else C=0;
+		M=-10^10;K=[-1];
+		for(I=0;R!=[];R=cdr(R),I++){
+			if(C>0&&car(R)[F[1]]!=F[2]) continue;
+			if(M<car(R)[0]){
+				M=car(R)[0];K=[I,car(R)];
+			}
+		}
+		return K;
+	}
 	R=reverse(R);
 	if(F[0]=="sub"){
 		for(S=F[1];S!=[];S=cdr(S))
@@ -8287,13 +8707,19 @@ def anal2sp(R,F)
 		return G;
 	}
 	if(F[0]=="+"){
-		for(G=[];R!=[];R=cdr(R))
-			G=cons([car(R)[0],car(R)[1]+F[1],car(R)[2]+F[2]],G);
+		L=length(F);
+		for(G=[];R!=[];R=cdr(R)){
+			for(S=[],I=L-1;I>0;I--) S=cons(car(R)[I]+F[I],S);
+			G=cons(cons(car(R)[0],S),G);
+		}
 		return G;
 	}
 	if(F[0]=="*"){
-		for(G=[];R!=[];R=cdr(R))
-			G=cons([car(R)[0],car(R)[1]*F[1]+car(R)[2]*F[2]],G);
+		L=length(F);
+		for(G=[];R!=[];R=cdr(R)){
+			for(S=0,I=1;I<L;I++) S+=car(R)[I]*F[I];
+			G=cons([car(R)[0],S],G);
+		}
 		return G;
 	}
 	if(F[0]=="mult"){
@@ -8351,8 +8777,11 @@ def anal2sp(R,F)
  P=["get",L]
     L=n        for variable x_n
     L=[m,n]    for residue [m,n]
+    L=[m,n,l]  for residue [m,n,l]
     L=[[m,n],[m',n']] for common spct
+ P=["eigen",I]   decomposition of A_I
  P=["get0",[m,n],[m',n']] for the sum of residues
+ P=["rest",[m,n]] restriction 
  P=["swap",[m,n]] for symmetry
  P=["perm",[...]] for symmetry
  P=["deg"]
@@ -8422,7 +8851,9 @@ def mc2grs(G,P)
 	}
 	if(type(P)<2) return G;
 	F=0;
-	if(type(P)==7||(type(P)==4&&type(P[0])<4)) P=[P];
+	if(type(P)==7||(type(P)==4&&
+		(type(P[0])<4||(type(P[0])==4&&length(P[0])==2&&type(P[0][0])<4&&type(P[1])<4))
+	  )) P=[P];
 	if((Dvi=getopt(dviout))!=1&&Dvi!=2&&Dvi!=-1) Dvi=0;
 	Keep=(Dvi==2)?1:0;
 	if(type(P)==4&&type(F=car(P))==7){
@@ -8452,6 +8883,56 @@ def mc2grs(G,P)
 			return R;
 		}
 		if(F=="show0"){
+			if(type(Fig=getopt(fig))>0){
+				PP=[[-1.24747,-5.86889],[1.24747,-5.86889],[3.52671,-4.8541],[5.19615,-3],
+				  [5.96713,-0.627171],[5.70634,1.8541],[4.45887,4.01478],[2.44042,5.48127],
+				  [0,6],[-2.44042,5.48127],[-4.45887,4.01478],[-5.70634,1.8541],
+				  [-5.96713,-0.627171],[-5.19615,-3],[-3.52671,-4.8541]];
+				PL=[[1.8,-5.2],[5.7,-1.7],[3.2,5],[-3.6,4.7],[2.2,3],[-2.8,2.8],
+				  [-1.5,-1.4],[-3.2,-2.5],[0.76,-1.4],[-2,0.2]];
+				PC=["black,dashed","green,dashed","red,dashed","blue,dashed",
+					"black","cyan","green","blue","red","magenta"];
+				N=["1","2","3","4","5","6","7","8","9","a","b","c","d","e","f"];
+				LL=[[1,2,3],[4,5,6],[7,8,9],[10,11,12],[7,10,13],[4,11,14],[5,8,15],[1,12,15],
+				  [2,9,14],[3,6,13]];
+				TB=str_tb("\\draw\n",TB);
+				if(type(Fig)==4){
+					if(type(car(Fig))==1){
+						PP=ptaffine(car(Fig)/12,PP);PL=ptaffine(car(Fig)/12,PL);
+						Fig=cdr(Fig);
+					}
+					if(Fig!=[]&&length(Fig)==10) PC=Fig;
+				}
+				for(R=mc2grs(G,"show0"|dviout=-1),I=0;R!="";I++){	/* 頂点 */
+					J=str_chr(R,0,",");
+					if(J>0){
+						S=str_cut(R,0,J-1);
+						R=str_cut(R,J+1,1000);
+					}else{
+						S=R;R="";
+					}
+					T=(str_chr(S,0,"1")==0)?"":"[red]";
+					str_tb(["node",T,"(",N[I],") at ",xypos(PP[I]),"{$",S,"$}\n"],TB);
+				}
+				for(S=PC,P=PL,I=0;I<4;I++){
+					for(J=I+1;J<5;J++,S=cdr(S),P=cdr(P)){			/* 線の番号 */
+						SS=car(S);
+						if((K=str_chr(SS,0,","))>0) SS=sub_str(SS,0,K-1);
+						str_tb(["node[",SS,"] at ",xypos(car(P)),
+							"{$[",rtostr(I),rtostr(J),"]$}\n"],TB);
+					}
+				}
+				str_tb(";\n",TB);
+				for(I=0;I<10;I++){		/* 線 */
+					S=car(PC);P0=car(PC);L0=car(LL);PC=cdr(PC);LL=cdr(LL);
+					C=[N[L0[0]-1],N[L0[1]-1],N[L0[2]-1]];
+					str_tb(["\\draw[",S,"] (", C[0],")--(",C[1],") (",
+						C[0],")--(",C[2],") (",C[1],")--(",C[2],");\n"],TB);
+				}
+				R=str_tb(0,TB);
+				if(TikZ==1&&Dvi!=-1) dviout(xyproc(R)|dviout=1,keep=Keep);
+				return R;
+			}
 			for(S="",L=[];G!=[];G=cdr(G)){
 				for(TL=[],TG=cdr(car(G));TG!=[];TG=cdr(TG)) TL=cons(car(TG)[0],TL);
 				TL=msort(TL,[-1,0]);
@@ -8529,6 +9010,7 @@ def mc2grs(G,P)
 						else S="A_{"+rtostr(T[0][0])+rtostr(T[0][1])+"}&"+S;
 					}
 					L=ltotex(R|opt="GRS",pre=S);
+					if(type(D=getopt(div))==1 || type(D)==4) L=divmattex(L,D);
 					if(Dvi>0) dviout(L|eq=0,keep=Keep);
 				}
 				return L; /* get all spct */
@@ -8540,14 +9022,14 @@ def mc2grs(G,P)
 					if(I[0]>I[0]){S=I;I=J;J=S;};
 					K=lsort(I,J,0);
 					if(length(K)==4){
-						S=sp2grs(G,["get0",[I,J]]);
+						S=mc2grs(G,["get0",[I,J]]);
 						return anal2sp(S,[["*",1,1],0]);
 					}
 					I=lsort(K,lsort(I,J,2),1);
 					S=lsort([0,1,2,3,4],K,1);
-					D=sp2grs(G,"deg");
+					D=mc2grs(G,"deg");
 					if(findin(4,S)<0) D=-D;
-					J=sp2grs(G,["get0",[I,S]]);
+					J=mc2grs(G,["get0",[I,S]]);
 					if(I[0]>S[0]) J=sp2grs(J,"swap");
 					return anal2sp(J,[["+",0,D],["*",-1,1]]);
 				}
@@ -8559,6 +9041,10 @@ def mc2grs(G,P)
 						if(car(PG)[0]==T) return (F=="get")?car(PG):cdr(car(PG));
 					return [];	/* get common spct */
 				}
+				if(length(T)==3){
+					T0=T;T=lsort([0,1,2,3,4],T,1);
+					if(length(T)!=2) return [];
+				}else T0=0;
 				if(T[0]>T[1]) T=[T[1],T[0]];
 				for(FT=0,PG=G;PG!=[];PG=cdr(PG)){
 					if(car(PG)[0][0]==T){
@@ -8570,9 +9056,118 @@ def mc2grs(G,P)
 				}
 				if(!FT) return [];
 				L=anal2sp(cdr(car(PG)),[["get1",FT],0]);
+				if(T0!=0){
+					if((K=mc2grs(G,"deg"))!=0){
+						if(T[1]!=4) K=-K;
+						R=reverse(L);
+						for(L=[];R!=[];R=cdr(R)) L=cons([car(R)[0],car(R)[1]+K],L);
+					}
+					T=T0;
+				}
 				return (F=="get")?cons(T,L):L;
 			}
 		}
+		if(F=="rest"||F=="eigen"||F=="rest0"||F=="rest1"){
+			if(F!="eigen") 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]]);
+						if(S!=[]) R=cons(cons([I,J],S),R);
+					}
+				}
+				R=reverse(R);
+				if(Dvi){
+					TB=str_tb(0,0);
+					if(F=="rest0"||F=="rest1"){
+						for(T=R;;){
+							TT=car(T);
+							S=rtostr(car(TT)[0])+rtostr(car(TT)[1]);
+							str_tb(["[",S,"]","&: "],TB);
+							for(TR=[],TT=cdr(TT);TT!=[];TT=cdr(TT))
+								TR=cons(car(TT)[1],TR);
+							for(TR=qsort(TR);TR!=[];TR=cdr(TR))
+								str_tb([s2sp(car(TR)|short=1,std=-1),"\\ \\ "],TB);
+							if((T=cdr(T))==[]) break;
+							str_tb("\\\\\n",TB);
+						}
+					}else{
+						TB=str_tb(0,0);
+						for(T=R;;){
+							TT=car(T);
+							S=rtostr(car(TT)[0])+rtostr(car(TT)[1]);
+							str_tb(["[",S,"]",":\\ "],TB);
+							for(TR=[],TT=cdr(TT);;){
+								T0=car(TT);
+								str_tb(["&",my_tex_form(car(T0)),"&&\\to\\ \n",
+									ltotex(cdr(T0)|opt="GRS")],TB);
+								if((TT=cdr(TT))==[]) break;
+								str_tb("\\\\\n",TB);
+							}
+							if((T=cdr(T))==[]) break;
+							str_tb("\\allowdisplaybreaks\\\\\n",TB);
+						}
+					}
+					R=texbegin("align*",str_tb(0,TB));
+					if(Dvi!=-1) dviout(R|keep=Keep);
+				}
+				return R;
+			}
+			I=P[1];
+			if(I[0]>I[1]) I=[I[1],I[0]];
+			L=lsort([0,1,2,3,4],I,1);
+			if(F=="rest"&&length(P)==3){
+				J=P[2];if(J[0]>J[1]) J=[J[1],J[0]];
+				L=lsort(L,J,1);
+				if(length(L)!=1) return 0;
+				return [mc2grs(G,["get0",I]),mc2grs(G,["get0",[I[0],J[0]],[I[1],J[1]]]),
+					mc2grs(G,["get0",[I[0],J[1]],[I[1],J[0]]]),mc2grs(G,["get0",[I[0],I[1],L[0]]])];
+			}
+			L=[[L[0],L[1]],[L[0],L[2]],[L[1],L[2]]];
+			if(F!="eigen"){
+				if(I==[0,4]) L=reverse(L);
+				else{
+					for(V=[],J=2;J>=0;J--){
+						if(L[J][0]==0) V=cons([L[J][1],J],V);
+						else{
+							for(K=4;K>=0;K--){
+								if(findin(K,L[J])<0){
+									V=cons([K,J],V);break;
+								} 
+							}
+						}
+					}
+					V=qsort(V);
+					L=[L[V[0][1]],L[V[1][1]],L[V[2][1]]];
+				}
+			}
+			for(LL=[],T=L;T!=[];T=cdr(T))
+				LL=cons(mc2grs(G,["get0",[I,car(T)]]),LL);
+			LL=reverse(LL);
+			for(R=[],Q=mc2grs(G,["get0",I]);Q!=[];Q=cdr(Q)){
+				for(T=[],J=2;J>=0;J--){
+					V=anal2sp(LL[J],["get1",(I[0]<L[J][0])?1:2,car(Q)[1]]);
+					if(F=="rest"){
+						if(I[0]==0){
+							if(I[1]!=4){
+								if(L[J][1]!=4) V=anal2sp(V,["+",-car(Q)[1]]);
+							}else if (L[J][0]!=2) V=anal2sp(V,["+",-car(Q)[1]]);
+						}else if(L[J][0]!=0) V=anal2sp(V,["+",-car(Q)[1]]);
+					}
+					T=cons(V,T);
+				}
+				R=cons(cons(car(Q)[1],T),R);
+			}
+			if(F=="rest0"||F=="rest1"){
+				for(L=[];R!=[];R=cdr(R)){
+					TR=cdr(car(R));
+					if(F=="rest1"&&chkspt(TR|opt="idx")==2) continue;
+					L=cons([car(R)[0],s2sp(chkspt(TR|opt=6))],L);
+				}
+				R=reverse(L);
+			}
+			return R;
+		}
 		if(F=="deg"){
 			for(S=I=0;I<3;I++){
 				for(J=I+1;J<4;J++){
@@ -8583,9 +9178,10 @@ def mc2grs(G,P)
 			}
 			return S/L[0];
 		}
-		if(F=="spct"){
+		if(F=="spct"||F=="spct1"){
+			K=(F=="spct")?5:6;
 			G=mc2grs(G,"get");
-			M=newmat(5,5);
+			M=newmat(5,K);
 			for(;G!=[];G=cdr(G)){
 				GT=car(G);I=GT[0][0];J=GT[0][1];
 				for(S=0,L=[],GT=cdr(GT);GT!=[];GT=cdr(GT)){
@@ -8602,19 +9198,29 @@ def mc2grs(G,P)
 					for(L=M[I][J];L!=[];L=cdr(L)) S+=car(L)^2;
 				}
 				M[I][I]=S;
+				if(K==6){
+					for(S=[],J=4;J>=0;J--)
+						if(I!=J) S=cons(M[I][J],S);
+					R=chkspt(S|opt=2);
+					M[I][5]=((L=length(R))>1)?s2sp(R[L-2]|short=1):"";
+				}
 			}
 			if(Dvi){
 				S=[];
 				for(I=4;I>=0;I--){
-					L=[M[I][I]];
+					L=(K==6)?[M[I][5]]:[];
+					L=cons(M[I][I],L);
 					for(J=4;J>=0;J--){
 						if(I==J) L=cons("",L);
 						else L=cons(s2sp([M[I][J]]),L);
 					}
 					S=cons(L,S);
 				}
-				S=cons([x0,x1,x2,x3,x4,"idx"],S);
-				M=ltotex(S|opt="tab",hline=[0,1,z],vline=[0,1,z-1,z],left=["","$x_0$","$x_1$","$x_2$","$x_3$","$x_4$"]);
+				T=(K==6)?["reduction"]:[];
+				S=cons(append([x0,x1,x2,x3,x4,"idx"],T),S);
+				M=ltotex(S|opt="tab",hline=[0,1,z],
+					vline=(K==6)?[0,1,z-2,z-1,z]:[0,1,z-2,z-1,z],
+					left=["","$x_0$","$x_1$","$x_2$","$x_3$","$x_4$"]);
 				if(Dvi>0) dviout(M|keep=Keep);
 			}
 			return M;
@@ -8877,7 +9483,7 @@ def mcmgrs(G,P)
 	Keep=(Dvi==2)?1:0;
 	if(type(P)==4 && type(F=car(P))==7){
 		if(F=="mult"){
-			for(P=cdr(P);P!=[];P=cdr(P)) G=os_md.mc2grs(G,car(P)|option_list=getopt());
+			for(P=cdr(P);P!=[];P=cdr(P)) G=mc2grs(G,car(P)|option_list=getopt());
 			return G;
 		}
 		if(F=="get"||F=="get0"){
@@ -10058,6 +10664,7 @@ def divmattex(S,T)
 	if(length(L0)>0) L=cons(reverse(L0),L);
 	L=lv2m(reverse(L));	/* get matrix */
 	if(T==0) return L;
+	if(type(T)==1) T=[T];
 	Size=size(L);S0=Size[0];
 	if(type(T[0])!=4){
 		S1=Size[1];
@@ -10332,6 +10939,7 @@ def tocsv(L)
 	S=str_tb(0,Tb);
 	if(type(EXE=getopt(exe))!=1&&EXE!=0&&type(EXE)!=7) return S;
 	if(type(F)!=7){
+		fcat(-1,0);
 		F="risaout";
 		if(EXE>=2&&EXE<=9) F+=rtostr(EXE);
 		F=DIROUTD+F+".csv";
@@ -11083,7 +11691,6 @@ def xyline(P,Q)
 
 def xylines(P)
 {
-/*	mycat([P,getopt()]);	*/
 	Lf=getopt(curve);
 	if(type(Lf)!=1) Lf=0;
 	SS=getopt(opt);
@@ -11242,9 +11849,357 @@ def saveproc(S,Out)
 	}
 }
 
+def xygrid(X,Y)
+{
+	for(RR=[],I=0,Z=X;I<2;I++){
+		U=Z[2];L=LL=[];M=Z[3];
+		if(Z[1]==1||Z[1]==-1){
+			if(type(M)==4) L=M;
+			else{
+				if(U*(-dlog(1-1/20)/dlog(10))>=M){
+					L=cons([1,2,1/10],L);
+					LL=cons([1,2,1/2],LL);
+				}else if(U*(-dlog(1-1/10)/dlog(10))>=M)
+					L=cons([1,2,1/5],L);
+				else if(U*(-dlog(1-1/4)/dlog(10))>=M)
+					L=cons([1,2,1/2],L);
+				if(U*(-dlog(1-1/50)/dlog(10))>=M){
+					L=cons([2,5,1/10],L);
+					LL=cons([2,5,1/2],LL);
+				}else if(U*(-dlog(1-1/25)/dlog(10))>=M)
+					L=cons([2,5,1/5],L);
+				else if(U*(-dlog(1-1/10)/dlog(10))>=M)
+					L=cons([2,5,1/2],L);
+				if(U*(-dlog(1-1/100)/dlog(10))>=M){
+					L=cons([5,10,1/10],L);
+					LL=cons([5,10,1/2],LL);
+				}
+				else if(U*(-dlog(1-1/50)/dlog(10))>=M)
+					L=cons([5,10,1/5],L);
+				else if(U*(-dlog(1-1/20)/dlog(10))>=M)
+					L=cons([5,10,1/2],L);
+				L=cons(L,cons(LL,[[[1,10,1]]]));
+			}
+			R=scale(L|scale=U);
+			if(Z[1]==-1){
+				for(LL=[];R!=[];R=cdr(R)){
+					for(L=[],T=car(R);T!=[];T=cdr(T)) L=cons(U-car(T),L);
+					LL=cons(reverse(L),LL);
+				}
+				R=reverse(LL);
+			}
+		}else if(Z[1]==0){
+			if(type(M)==4){
+				R=scale(M|f=x,scale=U);
+			}else{
+				V=0;
+				if(U/10>=M) V=1/10;
+				else if(U/5>=M) V=1/5;
+				else if(U/2>=M) V=1/2;
+				R=[];
+				if(V>0){
+					UU=U*V;
+					for(R=[],J=UU;J<U;J+=UU) R=cons(J,R);
+				}
+				if(V==1/10) L=[U/2];
+				else L=[];
+				R=cons(R,cons(L,[[0,U]]));
+			}
+		}else if(type(Z[1])==4){
+			R=Z[1];
+			if(length(R)==0||type(R[0])!=4) R=[[],[],R];
+		}else return 0;
+		K=length(R);
+		S=newvect(K);
+		for(J=0;J<K;J++){
+			for(S[J]=[],JJ=0;JJ<=Z[0];JJ+=U){
+				for(P=R[J];P!=[];P=cdr(P))
+					if(car(P)+JJ<=Z[0]) S[J]=cons(car(P)+JJ,S[J]);
+			}
+		}
+		for(J=0;J<K;J++) S[J]=lsort(S[J],[],1);
+		for(U=[],J=K-1;J>0;J--){
+			U=lsort(S[J],U,0);S[J-1]=lsort(S[J-1],U,1);
+		}
+		RR=cons(vtol(S),RR);
+		Z=Y;
+	}
+	if((Raw=getopt(raw))==1) return RR;
+	SS=[];
+	if(type(Sf=getopt(shift))==7){
+		Sx=Sf[0];Sy=Sf[1];
+	}else Sx=Sy=0;
+	for(I=0;I<2;I++){
+		for(S0=[],L=RR[I];L!=[];L=cdr(L)){
+			for(S=[],T=car(L);T!=[];T=cdr(T)){
+				if(S!=[]) S=cons(0,S);
+				if(I==0){
+					S=cons([X[0]+Sx,car(T)+Sy],S);
+					S=cons([Sx,car(T)+Sy],S);
+				}else{
+					S=cons([car(T)+Sx,Y[0]+Sy],S);
+					S=cons([car(T)+Sx,Sy],S);
+				}
+			}
+			S0=cons(S,S0);
+		}
+		SS=cons(reverse(S0),SS);
+	}
+	SS=reverse(SS);
+	if(Raw==2) return SS;
+	if(length(Y)<5) T=[["",""]];
+	else if(type(Y[4])==4) T=[Y[4]];
+	else T=[Y[4],Y[4]];
+	if(length(X[4])==4) T=cons([""],T);
+	else if(type(X[4])==4) T=cons(X[4],T);
+	else T=cons([X[4]],T);
+	for(Sx=Sy=[],I=0;I<2;I++){
+		TT=T[I];
+		for(V=SS[I];V!=[];V=cdr(V)){
+			Op=car(TT);
+			if(length(TT)>1) TT=cdr(TT);
+			if(car(V)==[]) continue;
+			if(Op=="") S=xylines(car(V));
+			else S=xylines(car(V)|opt=Op);
+			if(I==0) Sx=cons(S,Sx);
+			else Sy=cons(S,Sy);
+		}
+	}
+	for(S="",Sx=reverse(Sx), Sy=reverse(Sy);Sx!=[]&&Sy!=[];){
+		if(Sx!=[]){
+			S+=car(Sx);Sx=cdr(Sx);
+		}
+		if(Sy!=[]){
+			S+=car(Sy);Sy=cdr(Sy);
+		}
+	}
+	return S;
+}
+
+
+def addIL(I,L)
+{
+	if(I==0){
+		for(R=[];L!=[];L=cdr(L)) R=addIL(car(L),R);
+		return reverse(R);
+	}
+	if(type(In=getopt(in))==1){
+		if(In==-1){
+			J=JJ=I[1];I=I[0];
+			for(R=[];L!=[];L=cdr(L)){
+				J=lmin([car(L)[0],JJ]);
+				if(J>I) R=cons([I,J],R);
+				I=lmax([car(L)[1],I]);
+			}
+			if(I<JJ) R=cons([I,JJ],R);
+			return reverse(R);
+		}else{
+			for(;L!=[];L=cdr(L)){
+				if(car(L)[0]>I) return 0;
+				if(car(L)[1]>=I){
+					if(In==3) return car(L);
+					if(In==1||(I!=car(L)[0]&&I!=car(L)[1])) return 1;
+					return 2;
+				}
+			}
+			return 0;
+		}
+	}
+	I0=car(I);I1=I[1];
+	for(F=0,R=[];L!=[];L=cdr(L)){
+		if(I0>car(L)[1]){
+			R=cons(car(L),R);
+			continue;
+		}
+		if(I0<=car(L)[1]){
+			I0=lmin([I0,car(L)[0]]);
+			if(I1<car(L)[0]){
+				R=cons([I0,I1],R);
+				for( ;L!=[];L=cdr(L)) R=cons(car(L),R);
+				F=1;
+				break;
+			}
+			I1=lmax([I1,car(L)[1]]);
+		}
+	}
+	if(!F) R=cons([I0,I1],R);
+	return reverse(R);
+}
+
+def xy2curve(F,N,Lx,Ly,Lz,A,B)
+{
+	Raw=getopt(raw);
+	if(type(Gap=getopt(gap))==4){
+		MG=Gap[1];Gap=car(Gap);
+	}else MG=3;
+	if(type(Gap)!=1 && Gap!=0) Gap=0.7;
+	if(type(Dvi=getopt(dviout))<1) Dvi=0;
+	OL=[["dviout",Dvi]];
+	if(type(Opt=getopt(opt))<1) Opt=0;
+	else OL=cons(["opt",Opt],OL);
+	if(type(Sc=getopt(scale))!=1 && type(Sc)!=4) Sc=[1,1,1];
+	else if(type(Sc)!=4) Sc=[Sc,Sc,Sc];
+	else if(length(Sc)!=3) Sc=[Sc[0],Sc[1],Sc[1]];
+	M=diagm(3,Sc);
+	if(A!=0||B!=0){
+		if(type(A)==6) M=A;
+		else M=mrot([0,-B,-A]|deg=1)*M;
+		V=M*newvect(3,[x,y,z]);
+		Fx=compdf(V[0],[x,y,z],F);Fy=compdf(V[1],[x,y,z],F);Fz=compdf(V[2],[x,y,z],F);
+	}else{
+		for(I=0;I<3;I++){
+			if(type(T=F[I])!=4) T=f2df(T);
+			if(type(T)==4) T=cons(car(T)*Sc[I],cdr(T));
+			else T*=Sc[I];
+			if(I==0) Fx=T;
+			else if(I==1) Fy=T;
+			else Fz=T;
+		}
+	}
+	if(Raw==5||!Gap)
+		return (Dvi||!Gap)? xygraph([Fy,Fz],N,Lx,Ly,Lz|option_list=OL):[Fx,Fy,Fz];
+	R=xygraph([Fy,Fz],N,Lx,Ly,Lz|raw=2);
+	R0=cdr(car(R));R1=R[1];
+	for(LT=[];R0!=[];R0=cdr(R0),R1=cdr(R1))
+		if(car(R0)!=0) LT=cons([R1[0],R1[1]],LT);
+	LT=reverse(LT);
+	if(N<0){
+		Be=xylines(car(R)|curve=1,proc=3,close=-1);
+		LT=reverse(cdr(LT));
+		LT=reverse(cdr(LT));
+	}
+	else Be=xylines(car(R)|curve=1,proc=3);
+	Be=cdr(cdr(Be));
+	Be=lbezier(car(Be));
+	if(Raw==4) return [Be,LT,Lx];
+	X=ptcombz(Be,0,0);
+	Var=(length(Lx)==3)?car(Lx):x;
+	if(type(Eq=getopt(eq))!=1) Eq=0.01;
+	if(TikZ==1){
+		Gap/=10;Eq/=10;
+	}
+	for(R=[],XT=X;XT!=[];XT=cdr(XT)){
+		V=car(XT);
+		U=LT[V[0][0]];
+		T=U[0]*V[1][0]+U[1]*(1-V[1][0]);
+		VV=myfdeval(Fx,[Var,T]);
+		U=LT[V[0][1]];
+		T=U[0]*V[1][1]+U[1]*(1-V[1][1]);
+		VV-=myfdeval(Fx,[Var,T]);
+		if(abs(VV)<Eq) continue;
+		I=(VV<0)?0:1;
+		R=cons([V[0][I],V[1][I],V[0][1-I],V[1][1-I]],R);
+	}
+	R=qsort(R);
+	if(Raw==3) return [Be,R];
+    Db=newvect(L=length(Be));
+	for(I=0;I<L;I++) Db[I]=[];
+	for(TR=R;TR!=[];TR=cdr(TR)){
+    	V1=ptbezier(Be,[I=car(TR)[0],P=car(TR)[1]])[1];
+		V2=ptbezier(Be,[car(TR)[2],car(TR)[3]])[1];
+		T=dsqrt(1-dvangle(V1,V2)^2);
+		if(T<1/MG) T=MG;
+		GP=Gap/T;
+		W=GP/dnorm(V1);
+		Db[I]=addIL([P-W,P+W],Db[I]);
+		if(P-W<0 && I>0) Db[I-1]=addIL([P-W+1,1],Db[I-1]);
+		if(P+W>1 && I+1<L) Db[I+1]=addIL([0,P+W-1],Db[I+1]);
+	}
+	Db=vtol(Db);
+	for(Bf=[];Be!=[];Be=cdr(Be),Db=cdr(Db)){
+		if(car(Db)==[]) Bf=cons(car(Be),Bf);
+		else{
+			D=addIL([0,1],car(Db)|in=-1);
+			for(;D!=[];D=cdr(D))
+				Bf=cons(tobezier(car(Be)|inv=car(D)),Bf);
+		}
+	}
+	Bf=reverse(Bf);
+	if(Raw==2) return Bf;
+	OL=[];
+	if(Opt){
+		if(type(Opt)==4&&length(Opt)>1) OL=[["opt",Opt[0]],["cmd",Opt[1]]];
+		else OL=[["opt",Opt]];
+	}else OL=[];
+	S=xybezier(lbezier(Bf|inv=1)|option_list=OL);
+	if(Raw==1||!Dvi) return S;
+	return xyproc(S|dviout=Dvi);
+}
+
+def rungeKutta(F,N,Lx,Y,IY)
+{
+	if((Pr=getopt(prec))==1){
+		One=eval(exp(0));
+	}else{
+		One=1;Pr=0;
+	}
+	if((FL=getopt(last))!=1) FL=0;
+	if(length(Lx)>2){
+		V=car(Lx);Lx=cdr(Lx);
+	}else V=x;
+	if(Pr==0) Lx=[deval(Lx[0]),deval(Lx[1])];
+	else Lx=[eval(Lx[0]),eval(Lx[1])];
+	if(type(Y)==4){
+		if((Sing=getopt(single))==1||type(F)!=4)
+			F=append(cdr(Y),[F]);
+		L=length(Y);
+		for(TF=[];F!=[];F=cdr(F))
+			TF=cons(f2df(car(F)),TF);
+		F=reverse(TF);
+	}else{
+		L=1;
+		F=f2df(F);
+	}
+	if(getopt(val)==1) V1=1;
+	else V1=0;
+	H=(Lx[1]-Lx[0])/N;H2=H/2;
+	FV=findin(V,vars(F));
+	K=newvect(4);
+	if(L==1){
+		R=[[T=Lx[0],S=IY]];
+		if(!H) return R;
+		for(;;){
+			for(I=0;I<4;I++){
+				if(I==0)      W=[[V,T],[Y,S]];
+				else if(I==3) W=[[V,T+H],[Y,S+H*K[2]]];
+				else          W=[[V,T+H2],[Y,S+H2*K[I-1]]];
+				if(FV<0) W=cdr(W);
+				K[I]=Pr?myfeval(F,W)*One:myfdeval(F,W);
+			}
+			S+=(K[0]+2*K[1]+2*K[2]+K[3])*H/6;T+=H;
+			if(!FL) R=cons([deval(T),S],R);
+			if((T+H-Lx[1])*H>0) break;
+		}
+	}else{
+		T=Lx[0];
+		R=[cons(T,V1?[car(IY)]:IY)];
+		S=ltov(IY);
+		if(!H) return R;
+		for(;;){
+			for(I=0;I<4;I++){
+				if(I==0)      W=cons([V,T   ],lpair(Y,vtol(S)));
+				else if(I==3) W=cons([V,T+H ],lpair(Y,vtol(S+H*K[2])));
+				else          W=cons([V,T+H2],lpair(Y,vtol(S+H2*K[I-1])));
+				if(FV<0) W=cdr(W);
+				for(TK=[],TF=F;TF!=[];TF=cdr(TF)){
+					TK=cons(Pr?myfeval(car(TF),W)*One:myfdeval(car(TF),W),TK);
+				}
+				K[I]=ltov(reverse(TK));
+			}
+			S+=(K[0]+2*K[1]+2*K[2]+K[3])*H/6;T+=H;
+			TS=vtol(S);
+			if(V1) TS=[car(TS)];
+			if(!FL) R=cons(cons(deval(T),TS),R);
+			if((T+H-Lx[1])*H>0) break;
+		}
+	}
+	return FL?(V1?S[0]:S):reverse(R);
+}
+
 def xy2graph(F0,N,Lx,Ly,Lz,A,B)
 {
-	/* (x,y,z) -> ( -x sin A + y cos A, z cos B - x cos A sin B - y sin A sin B) */
+	/* (x,y,z) -> (z sin B + x cos A cos B + y sin A cos B, 
+	    -x sin A + y cos A, z cos B - x cos A sin B - y sin A sin B) */
 	if((Proc=getopt(proc))==1||Proc==2){
 		OPT0=[["proc",3]];
 	}else{
@@ -11678,22 +12633,100 @@ def xy2graph(F0,N,Lx,Ly,Lz,A,B)
 	if(Dvi<0) return Lout;
 }
 
+def orthpoly(N)
+{
+	F=0;
+	if(type(P=getopt(pol))==7){
+		for(L=["Le","Ge","Tc","2T","Ja","He","La","Se"];L!=[];L=cdr(L),F++)
+			if(str_str(P,car(L)|end=2)==0) break;
+	}else P=0;
+	if(type(D=N)==4) D=N[0];
+	if(!isint(D)||D<0) return 0;
+	if(F==0) return seriesHG([-D,D+1],[1],(1-x)/2,D);
+	if(F==1) return red(seriesHG([-D,D+2*N[1]],[N[1]+1/2],(1-x)/2,D)*binom(D+2*N[1]-1,D));
+	if(F==2) return seriesHG([-D,D],[1/2],(1-x)/2,D);
+	if(F==3){
+		if(D==0) return 0;
+		return orthpoly([D-1,1]|pol="Ge");
+	}
+	if(F==4) return red(seriesHG([-D,D+N[1]],[N[2]],x,D));
+	if(F==5){
+		for(S=I=1;I<=D;I+=2) S*=I;
+		if(iand(D,1)) return seriesHG([-(D-1)/2],[3/2],x^2/2,D-1)*x*S*(-1)^((D-1)/2);
+		else return seriesHG([-D/2],[1/2],x^2/2,D)*S*(-1)^(D/2);
+	}
+	if(F==6){
+		NN=(type(N)==4)?N[1]:0;
+		return red(seriesHG([-D],[NN+1],x,D)*binom(D+NN,D)); 
+	}
+	if(F==7){
+		NN=N[1];
+		for(S=1,I=1;I<=D;I++) S+=(-1)^I*binom(D,I)*binom(D+I,I)*sftpow(x,I)/sftpow(NN,I);
+		return S;
+	}
+	return 0;
+}
+
+def schurpoly(L)
+{
+	N=length(L);
+	for(R=[],I=1;L!=[];L=cdr(L),I++) R=cons(car(L)+N-I,R);
+	L=reverse(R);
+	if(type(X=getopt(var))!=4){
+		V=(type(X)>1)?X:"x";
+		for(X=[],I=0;I<N;I++) X=cons(makev([V,N-I]),X);
+	}
+	M=newmat(N,N);
+	for(I=0;I<N;I++)
+		for(J=0;J<N;J++) M[I][J]=X[I]^L[J];
+	P=det(M);
+	for(I=0;I<N;I++)
+		for(J=I+1;J<N;J++) P=sdiv(P,X[I]-X[J]);
+	return P;
+}
+
 def fouriers(A,B,X)
 {
+	if((Y=getopt(y))==0||type(Y)>0) Y=deval(Y);
+	else Y=0;
+	if((V=getopt(const))==0||type(V)>0){
+		V=myfeval(V,Y);
+		K=1;
+	}else K=0;
 	if(A!=[]&&type(car(A))>1){
-		for(C=[],I=A[1];I>=0;I--) C=cons(myfeval(car(A),I),C);
+		for(C=[],I=A[1];I>=K;I--) C=cons(myf2eval(car(A),I,Y),C);
+		if(K) C=cons(0,C);
 		A=C;
 	}
+	if(K){
+		if(A!=[]) A=cdr(A);
+		A=cons(V,A);
+	}
 	if(B!=[]&&type(car(B))>1){
-		for(C=[],I=B[1];I>0;I--) C=cons(myfeval(car(B),I),C);
+		for(C=[],I=B[1];I>0;I--) C=cons(myf2eval(car(B),I,Y),C);
 		B=C;
 	}
-	R=0;
+	L=length(B)+1;
+	if(length(A)>=L) L=length(A)+1;
+	if(type(Sum=getopt(sum))>0){
+		if(Sum==1) Sum=1-x;
+		else if(Sum==2) Sum=[(z__)/(3.1416*x),[z__,os_md.mysin,3.1416*x]];
+ 		else Sum=f2df(Sum);
+		C=[];
+		if(A!=[]){
+			C=cons(car(A),C);
+			A=cdr(A);
+		}
+		for(I=1;A!=[];A=cdr(A),I++) C=cons(car(A)*myf2eval(Sum,I/L,L),C);
+		A=reverse(C);
+		for(C=[],I=1;B!=[];B=cdr(B),I++) C=cons(car(B)*myf2eval(Sum,I/L,L),C);
+		B=reverse(C);
+	}
 	if(getopt(cpx)==1){
-	if(type(X=eval(X))>1) return todf([os_md.fouriers,[["cpx",1]]],[[A],[B],[X]]);
+		if(type(X=eval(X))>1) return todf([os_md.fouriers,[["cpx",1]]],[[A],[B],[X]]);
 		V=dexp(@i*X);
 		for(C=A,P=1,I=0;C!=[];C=cdr(C),I++){
-			R+=car(C)*P;
+			R+=S*car(C)*P;
 			P*=V;
 		}
 		V=dexp(-@i*X);
@@ -11738,7 +12771,7 @@ def mysin(Z)
 def mytan(Z)
 {
 	if(type(Z=eval(Z))>1) return todf(os_md.mytan,[Z]);
-	if((Im=imag(Z))==0) return dsin(Z);
+	if((Im=imag(Z))==0) return dtan(Z);
 	V=myexp(2*Z*@i);
 	return @i*(1-V)/(1+V);
 }
@@ -12209,10 +13242,6 @@ def compdf(F,V,G)
 {
 	FL=["abs","floor","rint","zeta","gamma","arg","real","imag","conj"];
 	FS=[os_md.abs,floor,rint,os_md.zeta,os_md.gamma,os_md.myarg,real,imag,conj];
-	if(type(V)==4){
-		for(;V!=[];V=cdr(V),G=cdr(G)) F=compdf(F,car(V),car(G));
-		return F;
-	}
 	if(type(F)==7){
 		if(str_str(F,"|")==0){
 			F="abs("+str_cut(F,1,str_len(F)-2)+")";
@@ -12233,7 +13262,20 @@ def compdf(F,V,G)
 	}
 	if(type(F)!=4) F=f2df(F);
 	if(type(G)!=4) G=f2df(G);
+	if(V==G) return F;	/* subst(F(V),V,G) */ 
 	VF=vars(F);VG=vars(G);
+	if(type(V)==4){
+		for(VT=[],VV=V;VV!=[];VV=cdr(VV)){
+			if(findin(car(VV),VF)>=0){
+				X=makenewv(append(VF,VG));
+				VF=cons(X,VF);
+				F=mysubst(F,[car(VV),X]);
+				VT=cons(X,VT);
+			}else VT=cons(car(VV),VT);
+		}
+		for(V=reverse(VT);V!=[];V=cdr(V),G=cdr(G)) F=compdf(F,car(V),car(G));
+		return F;
+	}
 	for(E=I=0;I<30;I++){
 		for(J=0;J<30;J++){
 			X=makev(["z__",I,J]);
@@ -12244,7 +13286,6 @@ def compdf(F,V,G)
 		if(E) break;
 	}
 	if(!E) return 0;
-	if(V==G) return F;	/* subst(F(V),V,G) */ 
 	if(type(G)<4) return mysubst(F,[V,G]);
 	if(type(F)<4) F=[F]; /* return compdf([X,[X,0,F]],V,G); */
 	F=mysubst(F,[V,X]);
@@ -12692,7 +13733,8 @@ def xygraph(F,N,LT,LX,LY)
 		}
 		V=reverse(NV);
 	}
-	if(getopt(raw)==1) return V;
+	if((Raw=getopt(raw))==1) return V;
+	if(Raw==2) return [V,LT];
 	OL=[["curve",1]];OLP=[];
 	if(type(C=getopt(ratio))==1){
 		OL=cons(["ratio",C],OL);OLP=cons(["ratio",C],OLP);
@@ -12950,7 +13992,11 @@ def polroots(L,V)
 	Lim=Lim2=[];
 	if(type(L)<4){
 		if(type(Lim=getopt(lim))==4){
-			if(type(Lim[0])!=4) Lim=[Lim];
+			if(type(Lim[0])!=4){
+				if(!isvar(Lim[0])) Lim=cons(V,[Lim]);
+				Lim=[Lim];
+			}
+			if(!isvar(Lim[0][0])) Lim=[cons(V,Lim)];
 			Lim=delopt(Lim,V|inv=1);
 			if(Lim!=[]){
 				Lim=Lim[0];
@@ -13012,7 +14058,7 @@ def polroots(L,V)
 	if(SS==0&&INIT==1){
 		SS=polroots(L,V|option_list=OL);
 		if(SS!=0) return SS;
-		for(C=0;SS==0&&C<4;C++){
+		for(C=0;SS==0&&C<5;C++){
 			I=(C==0)?1:(iand(random(),0xff)-0x80);
 			for(LL=[],K=length(L)-1;K>=0;K--){
 				for(Q=0,J=length(L)-1;J>=0;J--)
@@ -13038,6 +14084,7 @@ def polroots(L,V)
 	for(SS=[];R!=[];R=cdr(R)){
 		RS=(N==2)?[car(R)]:car(R);
 		for(I=0,L0=L[0];I<N-1;I++) L0=mysubst(L0,[V1[I],RS[I]]);
+		if(L0==0) return 0;
 		S0=polroots(L0,V[0]|option_list=OL);
 		if(type(S0)<2) return S0;
 		for(S=S0;S!=[];S=cdr(S)){
@@ -13230,7 +14277,7 @@ def cutf(F,X,VV)
 			if(car(V)!=[] && car(V)[0]<X) return myfeval(car(V)[1],Y);
 			return myfeval(F,Y);
 		}
-		if(X>car(V)[0]) continue;
+		if(car(V)==[]||X>car(V)[0]) continue;
 		if(X==car(V)[0]) return car(V)[1];
 		return myfeval(F,Y);
 	}
@@ -13257,8 +14304,8 @@ def periodicf(F,L,X)
 	if(type(L)==4) L=[eval(L[0]),eval(L[1])];
 	else L=eval(L);
 	if(isvar(X)){
-		Y=makenewv([X,V]);
-		Z=makenewv([X,Y,V]);
+		Y=makenewv([X,F]);
+		Z=makenewv([X,Y,F]);
 		return [Z,[Z,os_md.periodicf,[mysubst(F,[x,Y])],(type(L)==4)?[L]:L,[[Y,X]]]];
 	}
 	if(type(X)==4){
@@ -13666,6 +14713,38 @@ def draw_bezier(ID,IDX,B)
 	return 0;
 }
 
+
+/*
+def redbezier(L)
+{
+	V=newvect(4);ST=0;
+	for(R=[],I=0,T=L;T=[];T=cdr(T){
+		if(type(car(T))<4){
+			F=0;
+			if(I==3)
+			if(car(T)==0){
+			}else if(car(T)==1){
+			}else if(car(T)==-1){
+				if(I<3) V[I++]=ST;
+			}
+		}else if(I==3){
+			if(R==[] || car(R)!=1){
+				R=cons(V[0],R);
+				if(ST==0) ST=V[0];
+			}
+			for(J=1;J<3;J++) R=cons(V[J],R);
+			while((T=cdr(T))!=[]){
+				R=cons(car(T),R);
+				if(type(car(R))<4)
+			}
+		}else{
+			if(ST==0) ST=car(T);
+			V[I++]= car(T);
+		}
+	}
+}
+*/
+
 def lbezier(L)
 {
 	if((In=getopt(inv))==1||In==2||In==3){
@@ -14054,6 +15133,39 @@ def xycirc(P,R)
     return S+"}};\n";
 }
 
+def xypoch(W,H,R1,R2)
+{
+ 	if(H>R1||2*H>R2){
+		errno(0);
+		return;
+ 	}
+	if(type(Ar=getopt(ar))!=1) Ar=TikZ?0.25:2.5;
+	T1=dasin(H/R1);S1=R1*dcos(T1); 
+	T2=dasin(H/R2);S2=R2*dcos(T2);
+	T3=dasin(2*H/R2);S3=R2*dcos(T3);
+	S=xyline([R1,0],[W-R1,0]);
+	S+=xyang(R1,[W,0],-@pi,@pi-T1);
+	S+=xyline([S2,H],[W-S1,H]);
+	S+=xyang(R2,[0,0],T2,2*@pi-T3);
+	S+=xylines([[S3,-2*H],[W-H-R2,-2*H],[W-H-R2,2*H],[W-S3,2*H]]); 
+	S+=xyang(R2,[W,0],-@pi+T2,@pi-T3);
+	S+=xyline([W-T2,-H],[W-T2,-H]);
+	S+=xyang(R1,[0,0],0,2*@pi-T1);
+	S+=xyline([W-S2,-H],[S1,-H]);
+	if(Ar>0){
+		S+=xyang(Ar,[W/2,0],[0,0],8);
+		S+=xyang(Ar,[W/2,-2*H],[0,-2*H],8);
+		S+=xyang(Ar,[W/2-Ar,-H],[W,-H],8);
+		S+=xyang(Ar,[W/2-Ar,H],[W,H],8);
+		S+=xyang(Ar,[W-S3,2*H],[W-H-R2,2*H],8);
+	}
+	S+=xyput([R1,0,"$\\bullet$"]);
+	S+=xyput([0,0,"$\\times$"]);
+	S+=xyput([W,0,"$\\times$"]);
+	if(TikZ) S=str_subst(S,";\n\\draw","\n");
+	return S;
+}
+
 def ptaffine(M,L)
 {
 	if(type(L)!=4&&type(L)!=5){
@@ -14342,16 +15454,27 @@ def ptcopy(L,V)
 
 def	average(L)
 {
-	L=os_md.m2l(L|flat=1);
-	M0=M1=car(L);
-	for(I=SS=0, LT=L; LT!=[]; LT=cdr(LT), I++){
-		S+=(V=car(LT));
-		SS+=V^2;
-		if(V<M0)		M0=V;
-		else if(V>M1)	M1=V; 
+	if(getopt(opt)=="co"){
+		S0=average(L[0]);V0=car(S0);
+		S1=average(L[1]);V1=car(S1);
+		L0=os_md.m2l(L[0]|flat=1);
+		L1=os_md.m2l(L[1]|flat=1);
+		for(S=0;L0!=[];L0=cdr(L0),L1=cdr(L1))
+			S+=(car(L0)-V0)*(car(L1)-V1);
+		S/=S0[1]*S1[1]*S0[2];
+		S=[S,S0,S1];
+	}else{
+		L=os_md.m2l(L|flat=1);
+		M0=M1=car(L);
+		for(I=SS=0, LT=L; LT!=[]; LT=cdr(LT), I++){
+			S+=(V=car(LT));
+			SS+=V^2;
+			if(V<M0)		M0=V;
+			else if(V>M1)	M1=V; 
+		}
+		SS=dsqrt(SS/I-S^2/I^2);
+		S=[deval(S/I),SS,I,M0,M1];
 	}
-	SS=dsqrt(SS/I-S^2/I^2);
-	S=[deval(S/I),SS,I,M0,M1];
 	if(isint(N=getopt(sint))) S=sint(S,N);
 	return S;
 }
@@ -16088,6 +17211,30 @@ def conf1sp(M)
 	return P;
 }
 
+def confexp(S)
+{
+	if(type(S[0])==4){
+		for(E=[];S!=[];S=cdr(S))
+			E=cons(confexp(car(S),E));
+		return reverse(E);
+	}
+	V=x;E=[];
+	for(P=0,Q=[],ST=S;ST!=[];ST=cdr(ST)){
+		Q=cons(car(ST)[0],Q);
+		P+=car(ST)[1]/(V-car(ST)[0]);
+		P=red(P);
+	}
+	P=red(P*polbyroot(Q,V));
+	Q=cdr(reverse(Q));
+	for(I=(length(W=Q));I>=0;I--){
+		C=mycoef(P,I,V);
+		P-=C*polbyroot(W,V);
+		W=cdr(W);
+		E=cons(red(C),E);
+	}
+	return reverse(E);
+}
+
 def pgen(L,VV)
 {
 	if(type(L[0])<4) L=[L];
@@ -17213,7 +18360,7 @@ def integrate(P,X)
 			if(S!=RR) R=cons([[1,RR=S]],R);
 			for(V=FR=[];R!=[];R=cdr(R))
 				if(car(R)!=FR) V=cons(FR=car(R),V);
-			Var=varargs(V|all=1)[1];
+			Var=varargs(V|all=2);
 			for(S0=[x0,x1,x2,x3],S=[t,s,u,v,w];S0!=[]&&S!=[];){
 				if(findin(car(S0),Var)<0){
 					S0=cdr(S0); continue;
@@ -18344,10 +19491,10 @@ def linfrac01(X)
 
 def varargs(P)
 {
-	if((All=getopt(all))!=1) All=0;
+	if((All=getopt(all))!=1&&All!=2) All=0;
 	V=vars(P);
 	for(Arg=FC=[];V!=[];V=cdr(V)){
-		if(vtype(CV=car(V))==0&&All==1){
+		if(vtype(CV=car(V))==0&&All!=0){
 			Arg=lsort([CV],Arg,0);
 		}
 		if(vtype(CV)!=2) continue;
@@ -18364,7 +19511,8 @@ def varargs(P)
 			}
 		}
 	}
-	return [FC,Arg];
+	Arg=reverse(Arg);
+	return (All==2)?Arg:[reverse(FC),Arg];
 }
 
 def pfargs(P,X)
@@ -18790,7 +19938,7 @@ def init() {
 	}
 	if(Id>=0){
 		while((S=get_line(Id))!=0){
-			if(type(P=str_str(S,LS))==4 && (P0=str_char(S,P[1]+5,"="))>0){
+			if(type(P=str_str(S,LS))==4 && (P0=str_char(S,P[1]+4,"="))>0){
 				if(P[0]<5){
 					P0=str_chr(S,P0+1,"\"");
 					if(P0>0){
@@ -18814,15 +19962,12 @@ def init() {
 					else if(P[0]==7)	TikZ=SV;
 					else if(P[0]==8)	XYPrec=SV;
 					else if(P[0]==9)	XYcm=SV;
-					else if(P[0]==10)	XYcm=Canvas;
+					else if(P[0]==10)	Canvas=SV;
 				}
 			}
 		}
 		close_file(Id);
 	}
-	DIROUTD=str_subst(DIROUT,["%HOME%","%ASIRROOT%","%TMP%","%TEMP%"],
-		[getenv("HOME"),get_rootdir(),getenv("TMP"),getenv("TEMP")])
-		+((isMs())?"\\":"/");
 	chkfun(1,0);
 }