[BACK]Return to os_muldif.rr CVS log [TXT][DIR] Up to [local] / OpenXM / src / asir-contrib / packages / src

Diff for /OpenXM/src/asir-contrib/packages/src/os_muldif.rr between version 1.47 and 1.59

version 1.47, 2019/03/05 01:52:51 version 1.59, 2020/03/03 01:57:03
Line 1 
Line 1 
 /* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.46 2019/03/04 02:02:04 takayama Exp $ */  /* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.58 2020/02/25 02:47:35 takayama Exp $ */
 /* The latest version will be at ftp://akagi.ms.u-tokyo.ac.jp/pub/math/muldif  /* 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   scp os_muldif.[dp]* ${USER}@lemon.math.kobe-u.ac.jp:/home/web/OpenXM/Current/doc/other-docs
 */  */
 #define USEMODULE 1  #define USEMODULE 1
 /* #undef USEMODULE */  /* #undef USEMODULE */
   
 /*             os_muldif.rr (Library for Risa/Asir)  /*             os_muldif.rr (Library for Risa/Asir)
  *          Toshio Oshima (Nov. 2007 - Feb. 2019)   *          Toshio Oshima (Nov. 2007 - Mar. 2020)
  *   *
  *   For polynomials and differential operators with coefficients   *   For polynomials and differential operators with coefficients
  *   in rational funtions (See os_muldif.pdf)   *   in rational funtions (See os_muldif.pdf)
Line 59  localf countin$
Line 59  localf countin$
 localf mycoef$  localf mycoef$
 localf mydiff$  localf mydiff$
 localf myediff$  localf myediff$
   localf mypdiff$
   localf pTaylor$
   localf pwTaylor$
 localf m2l$  localf m2l$
 localf m2ll$  localf m2ll$
 localf mydeg$  localf mydeg$
Line 147  localf getel$
Line 150  localf getel$
 localf ptol$  localf ptol$
 localf rmul$  localf rmul$
 localf mtransbys$  localf mtransbys$
   localf trcolor$
 localf drawopt$  localf drawopt$
 localf execdraw$  localf execdraw$
 localf execproc$  localf execproc$
Line 171  localf myasin$
Line 175  localf myasin$
 localf myacos$  localf myacos$
 localf myatan$  localf myatan$
 localf mylog$  localf mylog$
   localf nlog$
 localf mypow$  localf mypow$
 localf scale$  localf scale$
 localf arg$  localf arg$
Line 260  localf expat$
Line 265  localf expat$
 localf polbyroot$  localf polbyroot$
 localf polbyvalue$  localf polbyvalue$
 localf pcoef$  localf pcoef$
   localf pmaj$
 localf prehombf$  localf prehombf$
 localf prehombfold$  localf prehombfold$
 localf sub3e$  localf sub3e$
Line 297  localf iscoef$
Line 303  localf iscoef$
 localf iscombox$  localf iscombox$
 localf sproot$  localf sproot$
 localf spgen$  localf spgen$
   localf spbasic$
 localf chkspt$  localf chkspt$
 localf cterm$  localf cterm$
 localf terms$  localf terms$
Line 337  localf divmattex$
Line 344  localf divmattex$
 localf dviout0$  localf dviout0$
 localf myhelp$  localf myhelp$
 localf isMs$  localf isMs$
   localf getline$
 localf showbyshell$  localf showbyshell$
 localf readcsv$  localf readcsv$
 localf tocsv$  localf tocsv$
Line 352  localf texsp$
Line 360  localf texsp$
 localf getbygrs$  localf getbygrs$
 localf mcop$  localf mcop$
 localf shiftop$  localf shiftop$
   localf shiftPfaff;
 localf conf1sp$  localf conf1sp$
 localf confexp$  localf confexp$
 localf confspt$  localf confspt$
   localf mcvm$
 localf s2csp$  localf s2csp$
 localf partspt$  localf partspt$
 localf pgen$  localf pgen$
Line 390  localf primroot$
Line 400  localf primroot$
 localf varargs$  localf varargs$
 localf ptype$  localf ptype$
 localf pfargs$  localf pfargs$
   localf regress$
 localf average$  localf average$
 localf tobig$  localf tobig$
 localf sint$  localf sint$
 localf frac2n$  localf frac2n$
   localf openGlib$
 localf xyproc$  localf xyproc$
 localf xypos$  localf xypos$
 localf xyput$  localf xyput$
Line 414  localf periodicf$
Line 426  localf periodicf$
 localf cmpf$  localf cmpf$
 localf areabezier$  localf areabezier$
 localf saveproc$  localf saveproc$
   localf xyplot$
 localf xygraph$  localf xygraph$
 localf xy2graph$  localf xy2graph$
 localf addIL$  localf addIL$
Line 469  extern Canvas$
Line 482  extern Canvas$
 extern ID_PLOT$  extern ID_PLOT$
 extern Rand$  extern Rand$
 extern LQS$  extern LQS$
 extern  SV=SVORG$  extern SV=SVORG$
 #endif  #endif
 static S_Fc,S_Dc,S_Ic,S_Ec,S_EC,S_Lc$  static S_Fc,S_Dc,S_Ic,S_Ec,S_EC,S_Lc$
 static S_FDot$  static S_FDot$
 extern AMSTeX$  extern AMSTeX$
 Muldif.rr="00190301"$  extern Glib_math_coordinate$
   extern Glib_canvas_x$
   extern Glib_canvas_y$
   Muldif.rr="00200302"$
 AMSTeX=1$  AMSTeX=1$
 TeXEq=5$  TeXEq=5$
 TeXLim=80$  TeXLim=80$
Line 782  def myediff(P,X)
Line 798  def myediff(P,X)
         return red(X*diff(P,X));          return red(X*diff(P,X));
 }  }
   
   def mypdiff(P,L)
   {
           if(type(P)>3) return map(os_md.mypdiff,P,L);
           for(Q=0;L!=[];L=cdr(L)){
                   Q+=mydiff(P,car(L))*L[1];
                   L=cdr(L);
           }
           return red(Q);
   }
   
   def pTaylor(S,X,N)
   {
           if(!isvar(T=getopt(time))) T=t;
           if(type(S)<4) S=[S];
           if(type(X)<4) X=[X];
           if(findin(T,varargs(S|all=2))>=0){
                   S=cons(z_z,S);X=cons(z_z,X);FT=1;
           }else FT=0;
           LS=length(S);
           FR=(getopt(raw)==1)?1:0;
           if(!FR) R=newvect(LS);
           else R=R1=[];
           for(L=[],I=0,TS=S,TX=X;I<LS;I++,TS=cdr(TS),TX=cdr(TX)){
                   if(!FR) R[I]=car(TX)+car(TS)*T;
                   else{
                           R=cons(car(TX),R);R1=cons(car(TS),R1);
                   }
                   L=cons(car(TS),cons(car(TX),L));
           }
           L=reverse(L);
           if(FR) R=[reverse(R1),reverse(R)];
           for(K=M=1;N>1;N--){
                   S=mypdiff(S,L);
                   K*=++M;
                   for(TS=S,I=0,R1=[];TS!=[];TS=cdr(TS),I++){
                           if(!FR) R[I]+=car(TS)*t^M/K;
                           else R1=cons(car(TS)/K,R1);
                   }
                   if(FR) R=cons(reverse(R1),R);
           }
           if(FT){
                   if(!FR){
                           S=newvect(LS-1);
                           for(I=1;I<LS;I++) S[I-1]=R[I];
                   }else{
                           for(S=[];R!=[];R=cdr(R)){
                                   S=cons(cdr(car(R)),S);
                           }
                           R=S;
                   }
                   R=subst(S,z_z,0);
           }
           return (FR&&!FT)?reverse(R):R;
   }
   
 def m2l(M)  def m2l(M)
 {  {
         if(type(M) < 4)          if(type(M) < 4)
Line 804  def m2l(M)
Line 875  def m2l(M)
   
 def mydeg(P,X)  def mydeg(P,X)
 {  {
         if(type(P) < 3)          if(type(P) < 3 && type(X)==2)
                 return deg(P,X);                  return deg(P,X);
         II = -1;          II=(type(X)==4)?-100000:-1;
         Opt = getopt(opt);          Opt = getopt(opt);
         if(type(P) >= 4){          if(type(P) >= 4){
                 S=(type(P) == 6)?size(P)[0]:0;                  S=(type(P) == 6)?size(P)[0]:0;
                 P = m2l(P);                  P = m2l(P);
                 for(I = 0, Deg = -3; P != []; P = cdr(P), I++){                  for(I = 0, Deg = -100000; P != []; P = cdr(P), I++){
                         if( (DT = mydeg(car(P),X)) == -2)                          if( (DT = mydeg(car(P),X)) == -2&&type(X)!=4)
                                 return -2;                                  return -2;
                         if(DT > Deg){                          if(DT > Deg){
                                 Deg = DT;                                  Deg = DT;
Line 822  def mydeg(P,X)
Line 893  def mydeg(P,X)
                 return (Opt==1)?([Deg,(S==0)?II:[idiv(II,S),irem(II,S)]]):Deg;                  return (Opt==1)?([Deg,(S==0)?II:[idiv(II,S),irem(II,S)]]):Deg;
         }          }
         P = red(P);          P = red(P);
         if(deg(dn(P),X) == 0)          if(type(X)==2){
                 return deg(nm(P),X);                  if(deg(dn(P),X) == 0)
                           return deg(nm(P),X);
           }else{
                   P=nm(red(P));
                   for(D=-100000,I=deg(P,X[1]);I>=0;I--){
                           if(TP=mycoef(P,I,X[1])){
                                   TD=mydeg(TP,X[0])-I;
                                   if(D<TD) D=TD;
                           }
                   }
                   return D;
           }
         return -2;          return -2;
 }  }
   
Line 1217  def vprod(V1,V2)
Line 1299  def vprod(V1,V2)
 def dnorm(V)  def dnorm(V)
 {  {
         if(type(V)<2) return dabs(V);          if(type(V)<2) return dabs(V);
           if((M=getopt(max))==1||M==2){
                   if(type(V)==5) V=vtol(V);
                   for(S=0;V!=[];V=cdr(V)){
                           if(M==2) S+=dabs(car(V));
                           else{
                                   if((T=dabs(car(V)))>S) S=T;
                           }
                   }
                   return S;
           }
         R=0;          R=0;
         if(type(V)!=4)          if(type(V)!=4)
                 for (I = length(V)-1; I >= 0; I--) R+= V[I]^2;                  for (I = length(V)-1; I >= 0; I--) R+= real(V[I])^2+imag(V[I])^2;
         else{          else{
                 if(type(V[0])>3){                  if(type(V[0])>3){
                         V=ltov(V[0])-ltov(V[1]);                          V=ltov(V[0])-ltov(V[1]);
                         return dnorm(V);                          return dnorm(V);
                 }                  }
                 for(;V!=[]; V=cdr(V))   R+=car(V)^2;                  for(;V!=[]; V=cdr(V))   R+=real(car(V))^2+imag(car(V))^2;
         }          }
         return dsqrt(R);          return dsqrt(R);
 }  }
Line 2249  def vgen(V,W,S)
Line 2341  def vgen(V,W,S)
 def mmc(M,X)  def mmc(M,X)
 {  {
         Mt=getopt(mult);          Mt=getopt(mult);
         if(type(M)==7) M=os_md.s2sp(M);          if(type(M)==7) M=s2sp(M);
         if(type(M)!=4||type(M[0])!=6) return 0;          if(type(M)!=4) return 0;
           if(type(M[0])<=3){
                   for(RR=[];M!=[];M=cdr(M)) RR=cons(mat([car(M)]),RR);
                   M=reverse(RR);
           }
         if(type(M[0])!=6){                      /* spectre type -> GRS */          if(type(M[0])!=6){                      /* spectre type -> GRS */
                 G=s2sp(M|std=1);                  G=s2sp(M|std=1);
                 L=length(G);                  L=length(G);
Line 3413  def lsort(L1,L2,T)
Line 3509  def lsort(L1,L2,T)
                                         }else{                                          }else{
                                                 for(I=0;LT!=[];I++,LT=cdr(LT))                                                  for(I=0;LT!=[];I++,LT=cdr(LT))
                                                         if(findin(I,C1)<0) RT=cons(car(LT),RT);                                                          if(findin(I,C1)<0) RT=cons(car(LT),RT);
                                                 RT=reverse(RT);  
                                         }                                          }
                                         R=cons(RT,R);                                          R=cons(reverse(RT),R);
                                 }                                  }
                                 return reverse(R);                                  return reverse(R);
                         }                          }
Line 3771  def lgcd(L)
Line 3866  def lgcd(L)
         return [];          return [];
 }  }
   
   #if0
 def llcm(L)  def llcm(L)
 {  {
           if(type(L)==5||type(L)==6) L=m2l(L);
           if(type(L)<4) L=[L];
         if(type(L)==4){          if(type(L)==4){
                 F=getopt(poly);                  F=getopt(poly);
                 V=car(L);                  V=car(L);
Line 3786  def llcm(L)
Line 3884  def llcm(L)
                 if(F!=1&&V<0) V=-V;                  if(F!=1&&V<0) V=-V;
                 return V;                  return V;
         }          }
         else if(type(L)==5||type(L)==6)  
                 return llcm(m2l(L)|option_list=getopt());  
         return [];          return [];
 }  }
   #else
   def llcm(R)
   {
           if(type(R)==4||type(R)==5) R=m2l(R);
           if(type(R)<4) R=[R];
           if(type(R)!=4) return 0;
           V=getopt(poly);
           if(type(V)<1){
                   for(L=R;L!=[];L=cdr(L)){
                           if(type(car(L))>1){
                                   V=1; break;
                           }
                   }
           }
           if(getopt(dn)!=1){
                   for(L=[];R!=[];R=cdr(R)) if(R!=0) L=cons(1/car(R),L);
                   R=L;
           }
           P=1;
           if(type(V)<1){
                   for(;R!=[];R=cdr(R)){
                           if(!(TL=car(R))) continue;
                           else P=ilcm(P,dn(TL));
                   }
                   return P;
           }
           for(;R!=[];R=cdr(R)){
                   if(!car(R)) continue;
                   D=dn(red(car(R)));
                   N=red(P/D);
                   if(type(V)<2){
                           if(type(N)!=3) continue;
                           P*=dn(N);
                           continue;
                   }
                   if(ptype(N,V)>2){
                           L=fctr(dn(N));
                           for(;L!=[];L=cdr(L)){
                                   if(ptype(car(L)[0],V)<2) continue;
                                   P*=car(L)[0]^car(L)[1];
                           }
                   }
           }
           return P;
   }
   #fi
   
 def ldev(L,S)  def ldev(L,S)
 {  {
Line 3865  def lnsol(VV,L)
Line 4007  def lnsol(VV,L)
   
 def ladd(X,Y,M)  def ladd(X,Y,M)
 {  {
           if(Y==0){
                   Y=X[1];X=X[0];
           }
         if(type(Y)==4) Y=ltov(Y);          if(type(Y)==4) Y=ltov(Y);
         if(type(X)==4) X=ltov(X);          if(type(X)==4) X=ltov(X);
         return vtol(X+M*Y);          return vtol(X+M*Y);
Line 4574  def mtransbys(FN,F,LL)
Line 4719  def mtransbys(FN,F,LL)
         return call(FN, cons(F,LL)|option_list=Opt);          return call(FN, cons(F,LL)|option_list=Opt);
 }  }
   
   def trcolor(S)
   {
           if(type(S)!=7) return S;
           return ((I=findin(S,LCOPT))>=0)?COLOPT[I]:0;
   }
   
 def drawopt(S,T)  def drawopt(S,T)
 {  {
         if(type(S)!=7) return -1;          if(type(S)!=7) return -1;
Line 4605  def drawopt(S,T)
Line 4756  def drawopt(S,T)
         return -1;          return -1;
 }  }
   
   def openGlib(W)
   {
           extern Glib_canvas_x;
           extern Glib_canvas_y;
           extern Glib_math_coordinate;
   
           if(W==0){
                   glib_clear();
                   return;
           }
           if(type(W)==4&&length(W)==2){
                   Glib_canvas_x=W[0];
                   Glib_canvax_y=W[1];
           }
           Glib_math_coordinate=1;
           if(getopt(null)!=1) return glib_open();
   }
   
 def execdraw(L,P)  def execdraw(L,P)
 {  {
         if((Proc=getopt(proc))!=1) Proc=0;          if((Proc=getopt(proc))!=1) Proc=0;
Line 4857  def execdraw(L,P)
Line 5026  def execdraw(L,P)
                                                 LOut=cons(T[2],Out);                                                  LOut=cons(T[2],Out);
                                         }                                          }
                                 }                                  }
                           }else if(T[0]==6){      /* plot */
                                   F++;
                                   if((T1=findin(T[1],LCOPT))>-1) T1=COLOPT(T1);
                                   else if(type(T1)!=1 && T1!=0) T1=0xffffff;
                                   for(T2=ptaffine(M,T[2]|option_list=Org);T2!=[];T2=cdr(T2))
                                           draw_obj(Id,Ind,[rint(car(T2)[0]),rint(car(T2)[1])],T1);
                         }else if(Proc==1&&type(T[0])==2){                          }else if(Proc==1&&type(T[0])==2){
                                 if(length(T)<3) call(T[0],T[1]);                                  if(length(T)<3) call(T[0],T[1]);
                                 else call(T[0],T[1]|option_list=T[2]);                                  else call(T[0],T[1]|option_list=T[2]);
Line 4936  def execdraw(L,P)
Line 5111  def execdraw(L,P)
                                         if(P[0]==2)     dviout(T[2]|option_list=T[1]);                                          if(P[0]==2)     dviout(T[2]|option_list=T[1]);
                                         else LOut=cons(T[2],Out);                                          else LOut=cons(T[2],Out);
                                 }                                  }
                           }else if(T[0]==6){      /* plot */
                                   F++;
                                   if(type(T[1])==7) T1=[T[1],"."];
                                   else T1=".";
                                   for(T2=ptaffine(M,T[2]|option_list=Org);T2!=[];T2=cdr(T2))
                                           str_tb(xypos([car(T2)[0],car(T2)[1],T1]),Out);
                         }else if(T[0]==-2)                          }else if(T[0]==-2)
                                 str_tb(["%",T[1],"\n"],Out);                                  str_tb(["%",T[1],"\n"],Out);
                         else if(Proc==1&&type(T[0])==2){                           else if(Proc==1&&type(T[0])==2){
                                 if(length(T)<3) call(T[0],T[1]);                                  if(length(T)<3) call(T[0],T[1]);
                                 else call(T[0],T[1]|option_list=T[2]);                                  else call(T[0],T[1]|option_list=T[2]);
                         }                          }
Line 5105  def mmulbys(FN,P,F,L)
Line 5286  def mmulbys(FN,P,F,L)
   
 def appldo(P,F,L)  def appldo(P,F,L)
 {  {
           if(getopt(Pfaff)==1){
                   L = vweyl(L);
                   X = L[0]; DX = L[1];
                   for(I=mydeg(P,DX);I>0;I--){
                           if(!(TP=mycoef(P,D,DX))) continue;
                           P=red(P+TP*(muldo(D^(I-1),F,L)-D^I));
                   }
                   return P;
           }
         if(type(F) <= 3){          if(type(F) <= 3){
                 if(type(L) == 4 && type(L[0]) == 4)                  if(type(L) == 4 && type(L[0]) == 4)
                         return applpdo(P,F,L);                          return applpdo(P,F,L);
Line 5269  def mce(P,L,V,R)
Line 5459  def mce(P,L,V,R)
 {  {
         L = vweyl(L);          L = vweyl(L);
         X = L[0]; DX = L[1];          X = L[0]; DX = L[1];
         P = sftexp(laplace1(P,L),L,V,R);          P = sftexp(laplace1(P,L),L,V,R|option_list=getopt());
         return laplace(P,L);          return laplace(P,L);
 }  }
   
 def mc(P,L,R)  def mc(P,L,R)
 {  {
         return mce(P,L,0,R);          return mce(P,L,0,R|option_list=getopt());
 }  }
   
 def rede(P,L)  def rede(P,L)
Line 5452  def mulpdo(P,Q,L);
Line 5642  def mulpdo(P,Q,L);
   
 def transpdosub(P,LL,K)  def transpdosub(P,LL,K)
 {  {
           if(type(P)>3) return
   #ifdef  USEMODULE
                   mtransbys(os_md.transpdosub,P,[LL,K]);
   #else
                   mtransbys(transpdosub,P,[LL,K]);
   #endif
         Len = length(K)-1;          Len = length(K)-1;
         if(Len < 0 || P == 0)          if(Len < 0 || P == 0)
                 return P;                  return P;
Line 5477  def transpdosub(P,LL,K)
Line 5673  def transpdosub(P,LL,K)
   
 def transpdo(P,LL,K)  def transpdo(P,LL,K)
 {  {
         if(type(K[0]) < 4)  
                 K = [K];  
         Len = length(K)-1;          Len = length(K)-1;
         K1=K2=[];          K1=K2=[];
         if(type(LL)!=4) LL=[LL];          if(type(LL)!=4) LL=[LL];
         if(type(LL[0])!=4) LL=[LL];          if(type(LL[0])!=4) LL=[LL];
           if(type(car(K)) < 4 && length(LL)!=length(K)) K = [K];
         if(getopt(ex)==1){          if(getopt(ex)==1){
                 for(LT=LL, KT=K; KT!=[]; LT=cdr(LT), KT=cdr(KT)){                  for(LT=LL, KT=K; KT!=[]; LT=cdr(LT), KT=cdr(KT)){
                         L = vweyl(LL[J]);                          L = vweyl(LL[J]);
Line 5491  def transpdo(P,LL,K)
Line 5686  def transpdo(P,LL,K)
                 }                  }
                 K2=append(K1,K2);                  K2=append(K1,K2);
         }else{          }else{
                   if(length(LL)==length(K) && type(car(K))!=4){
                           for(DV=V=TL=[],J=length(LL)-1;J>=0;J--){
                                   TL=cons(vweyl(LL[J]),TL);
                                   V=cons(car(TL)[0],V);
                                   DV=cons(car(TL)[1],DV);
                           }
                           LL=TL;
                           if(type(RK=solveEq(K,V|inv=1))!=4) return TK;
                           if(!isint(Inv=getopt(inv))) Inv=0;
                           if(iand(Inv,1)){J=K;K=RK;RK=J;}
                           M=jacobian(RK,V|mat=1);
                           M=mulsubst(M,[V,K]|lpair=1);
                           RK=vtol(M*ltov(DV));
                           if(Inv>1) return RK;
                           K=lpair(K,RK);
                   }
                 for(J = length(K)-1; J >= 0; J--){                  for(J = length(K)-1; J >= 0; J--){
                         L = vweyl(LL[J]);                          L = vweyl(LL[J]);
                         if(L[0] != K[J][0])                          if(L[0]!= K[J][0]) K1=cons([L[0],K[J][0]],K1);
                                 K1 = cons([L[0],K[J][0]],K1);  
                         K2 = cons(K[J][1],K2);                          K2 = cons(K[J][1],K2);
                 }                  }
                 P = mulsubst(P, K1);                  P = mulsubst(P, K1);
Line 6265  def baseODE(L)
Line 6475  def baseODE(L)
         }          }
         if(type(To=getopt(to))<2||type(To)>4) To=0;          if(type(To=getopt(to))<2||type(To)>4) To=0;
         else if(!isvar(To)){          else if(!isvar(To)){
                 if(type(To)!=4) To=cons(red(To),cdr(Var));                  if(type(To)!=4){
                           To=red(To);
                           for(K=0;K<length(Var);K++){
                                   I=mydeg(nm(To),Var[K]);J=mydeg(dn(To),Var[K]);
                                   if(I+J>0&&I<2&&J<2) break;
                           }
                           if(K==length(Var)) return -9;
                           J=To;
                           for(To=[],I=length(Var)-1;I>=0;I--)
                                   if(I!=K) To=cons(Var[I],To);
                           To=cons(J,To);
                   }
                 if(type(To)==4){                  if(type(To)==4){
                         if(type(car(To))==4){                          if(type(car(To))==4){
                                 R=1;To=car(To);                                  R=1;To=car(To);
                         }else R=0;                          }else R=0;
                         if(type(IL=solveEq(To,Var|inv=1))!=4) return -7;                          if(type(IL=solveEq(To,Var|inv=1))!=4) return IL;
                         if(R==1){                          if(R==1){
                                 R=To;To=IL;IL=R;                                  R=To;To=IL;IL=R;
                         }                          }
Line 6296  def baseODE(L)
Line 6517  def baseODE(L)
                 }                  }
         }          }
         if(F==-3) return [Var,L];          if(F==-3) return [Var,L];
         for(I=0;I<M;I++) L=subst(L,Var[I], makev([SV[I],0]));          for(I=0;I<M;I++) L=subst(L,Var[I],makev([SV[I],0]));
         if(TeX){          if(TeX){
                 for(TL=L,I=0;I<M;I++)                  for(TL=L,I=0;I<M;I++)
                         TL=subst(TL,makev([SV[I],0]),Var[I]);                          TL=subst(TL,makev([SV[I],0]),Var[I]);
Line 6320  def baseODE(L)
Line 6541  def baseODE(L)
                 TL=cons(nm(red(T)),TL);                  TL=cons(nm(red(T)),TL);
         }          }
         if(isvar(To)){          if(isvar(To)){
                 T=rtostr(T);                  T=rtostr(To);
                 IT=findin(T,SV);                  IT=findin(T,SV);
                 if(IT>=0 && IT<M){                  if(IT>=0 && IT<M){
                         R=[SV[IT]];                          R=[SV[IT]];
Line 6351  def baseODE(L)
Line 6572  def baseODE(L)
                 TL=RL;                  TL=RL;
         }          }
         L=append(TL,L);          L=append(TL,L);
         for(I=0;I<M;I++) L=subst(L,Var[I],makev([SV[I]]));          for(I=0;I<M;I++) L=subst(L,makev([SV[I],0]),Var[I]);
         for(V=VV=[],I=0;I<M;I++){          for(V=VV=[],I=0;I<M;I++){
                 for(J=0;J<M;J++) V=cons(J?makev([SV[I],J]):makev([SV[I]]),V);                  for(J=0;J<M;J++) V=cons(J?makev([SV[I],J]):makev([SV[I]]),V);
                 if(!I||In) V=cons(makev([SV[0],M]),V);                  if(!I||In) V=cons(makev([SV[0],M]),V);
Line 6485  def toeul(F,L,V)
Line 6706  def toeul(F,L,V)
         L = vweyl(L);          L = vweyl(L);
         X = L[0]; DX = L[1];          X = L[0]; DX = L[1];
         I = mydeg(F,DX);          I = mydeg(F,DX);
         if(V == "infty"){          if(getopt(raw)!=1){
                 for(II=I; II>=0; II--){                  for(II=I; II>=0; II--){
                         J = mydeg(P=mycoef(F,I,DX),X);                          J = mydeg(P=mycoef(F,II,DX),X);
                         if(II==I) S=II-J;                          if(II==I) S=II-J;
                         else if(P!=0 && II-J>S) S=II-J;                          else if(P!=0 && II-J>S) S=II-J;
                 }                  }
                 F *= X^S;                  F *= X^S;
                 R = 0;          }
                 for( ; I >= 0; I--)          if(V == "infty"){
                   for(R=0; I >= 0; I--)
                          R += red((mysubst(mycoef(F,I,DX),[X,1/X])*(x*DX)^I));                           R += red((mysubst(mycoef(F,I,DX),[X,1/X])*(x*DX)^I));
                 return(subst(pol2sft(R,DX),DX,-DX));                  return(subst(pol2sft(R,DX),DX,-DX));
         }          }
         F = subst(F,X,X+V);          for(R=0; I >= 0; I--)
         for(II=I; II>=0; II--){  
                 J = mymindeg(P=mycoef(F,II,DX),X);  
                 if(II==I) S=II-J;  
                 else if(P!=0 && II-J>S) S=II-J;  
         }  
         F *= X^S;  
         R = 0;  
         for( ; I >= 0; I--)  
                 R += (red(mycoef(F,I,DX)/X^I))*DX^I;                  R += (red(mycoef(F,I,DX)/X^I))*DX^I;
         return pol2sft(R,DX);          return pol2sft(R,DX);
 }  }
Line 6541  def fromeul(P,L,V)
Line 6755  def fromeul(P,L,V)
                 S = DX*(S*X + mydiff(S,DX));                  S = DX*(S*X + mydiff(S,DX));
                 R += mycoef(P,J,DX)*S;                  R += mycoef(P,J,DX)*S;
         }          }
         while(mycoef(R,0,X) == 0)          if(getopt(raw)!=1){
                 R = tdiv(R,X);                  while(mycoef(R,0,X) == 0)
                           R = tdiv(R,X);
           }
         if(V != "infty" && V != 0)          if(V != "infty" && V != 0)
                 R = mysubst(R,[X,X-V]);                  R = mysubst(R,[X,X-V]);
         return R;          return R;
Line 6551  def fromeul(P,L,V)
Line 6767  def fromeul(P,L,V)
 def sftexp(P,L,V,N)  def sftexp(P,L,V,N)
 {  {
         L = vweyl(L); DX = L[1];          L = vweyl(L); DX = L[1];
         P = mysubst(toeul(P,L,V),[DX,DX+N]);          P = mysubst(toeul(P,L,V|opt_list=getpt()),[DX,DX+N]);
         return fromeul(P,L,V);          return fromeul(P,L,V|opt_list=getopt());
 }  }
   
   
Line 6947  def expat(F,L,V)
Line 7163  def expat(F,L,V)
   
 def polbyroot(P,X)  def polbyroot(P,X)
 {  {
           if(isvar(V=getopt(var))&&length(P)>1&&isint(car(P))){
                   for(Q=[],I=car(P);I<=P[1];I++) Q=cons(makev([V,I]),Q);
                   P=Q;
           }
         R = 1;          R = 1;
         while(length(P)){          while(length(P)){
                 R *= X-car(P);                  R *= X-car(P);
Line 7135  def pcoef(P,L,Q)
Line 7355  def pcoef(P,L,Q)
         return Coef;          return Coef;
 }  }
   
   def pmaj(P)
   {
           if(type(P)==4){
                   Opt=delopt(getopt(),"var"|inv=1);
                   P=mtransbys(os_md.pmaj,P,[]|optilon_list=Opt);
                   if(Opt==[]) return P;
                   X=Opt[0][1];
                   D=mydeg(P,X);
                   for(S=0;D>=0;D--) S+=lmax(mycoef(P,D,X))*X^D;
                   return S;
           }
           V=vars(P);
           if(!(K=length(V))) return abs(P);
           for(R=0,D=deg(P,X=V[0]);D>=0;D--){
                   Q=coef(P,D,X);
                   if(Q!=0) R+=((type(Q)>1)?pmaj(Q):abs(Q))*X^D;
           }
           if(isvar(Y=getopt(var))) for(;V!=[];V=cdr(V)) R=subst(R,car(V),Y);
           return R;
   }
   
 def prehombf(P,Q)  def prehombf(P,Q)
 {  {
         if((Mem=getopt(mem))!=1 && Mem!=-1)          if((Mem=getopt(mem))!=1 && Mem!=-1)
Line 8272  def spgen(MO)
Line 8513  def spgen(MO)
         }else{          }else{
                 L0=0; L1=MO+1;                  L0=0; L1=MO+1;
         }          }
         if(MO<=0){          if(M0<=0){
                 MO=-MO;                  MO=-MO;
                 if(iand(MO,1)==1) return [];                  if(iand(MO,1)==1) return [];
                 if(MO>1){                  MO=MO/2;
                         if(isMs()==0) return [];                  B=spbasic(-2*MO,0|str=1);
                         Cmd="okubo "+rtostr(-MO);                  if(L1<3) L1=MO+4;
                         MO/=2;  
                         if(L1>0) Cmd=Cmd+"+"+rtostr(L0)+"-"+rtostr(L1);  
                         else L1=MO+4;  
                         Cmd=Cmd+" B";  
                         Id=getbyshell(Cmd);  
                         if(Id<0) return [];  
                         B=[];  
                         while((S=get_line(Id)) !=0){  
                                 P0=str_chr(S,1,":")+1;  
                                 if(P0>1){  
                                         P1=str_chr(S,P,"\n");  
                                         if(P1<0) P1=str_len(S);  
                                         B=cons(sub_str(S,P0,P1-1),B);  
                                 }  
                         }  
                         close_file(Id);  
                 }else{  
                         MO/=2;  
                         if(L1<=1) L1=MO+4;  
 BB=[  
 ["11,11,11,11","111,111,111","1^4,1^4,22","1^6,222,33"],  
 ["11,11,11,11,11","1^4,1^4,211","211,22,22,22","1^6,2211,33",  
 "2211,222,222","22211,2^4,44","2^511,444,66","1^4,22,22,31",  
 "2^5,3331,55","1^5,1^5,32","1^8,332,44","111,111,21,21","1^5,221,221"],  
 ["11,11,11,11,11,11","1^4,1^4,1^4","1^4,22,22,22","111,111,111,21",  
 "1^6,21^4,33","21^4,222,222","221^4,2^4,44","2^41^4,444,66",  
 "1^5,1^5,311","1^8,3311,44","1^6,222,321","321,33,33,33",  
 "3321,333,333","33321,3^4,66","3^721,666,99","2^5,3322,55",  
 "1^6,1^6,42","222,33,33,42","1^a,442,55","1^6,33,33,51",  
 "222,222,33,51","1^9,333,54","2^7,554,77","1^5,2111,221",  
 "2^41,333,441","1^7,2221,43","211,211,22,22","2211,2211,222",  
 "22211,22211,44","1^4,211,22,31","2^411,3331,55","1^4,1^4,31,31",  
 "22,22,22,31,31","1^7,331,331","2221,2221,331","111,21,21,21,21"],  
 ["11,11,11,11,11,11,11","111,111,111,111","1^6,1^6,33",  
 "1^6,222,222","222,33,33,33","1^5,1^5,221",  
 "1^4,211,22,22","1^4,1^4,22,31","22,22,22,22,31",  
 "111,111,21,21,21","21^6,2^4,44","2221^6,444,66",  
 "1^6,222,3111","3111,33,33,33","33111,333,333",  
 "333111,3^4,66","3^5111,666,99","2^5,33211,55",  
 "1^8,3221,44","3222,333,333","33222,3^4,66",  
 "3^4222,666,99","1^6,1^6,411","222,33,33,411",  
 "1^a,4411,55","2^4,2^4,431","431,44,44,44",  
 "2^6,4431,66","4431,444,444","44431,4^4,88",  
 "4^531,888,cc","1^a,433,55","1^7,1^7,52",  
 "1^c,552,66","3^4,444,552","1^8,2^4,53",  
 "1^8,44,44,71","3^5,555,771","21^4,2211,222",  
 "221^4,22211,44","2221^4,3331,55","1^6,2211,321",  
 "2^411,3322,55","1^7,322,331","2211,33,33,42",  
 "3^42,4442,77","2211,222,33,51","3^51,5551,88",  
 "2^611,554,77","2221,2221,322","2^41,2^41,54",  
 "1^5,2111,2111","222111,333,441","1^7,22111,43",  
 "1^5,1^5,41,41","1^9,441,441","22111,2221,331",  
 "1^5,221,32,41","221,221,221,41","211,211,211,22",  
 "2211,2211,2211","1^4,211,211,31","211,22,22,31,31",  
 "1^4,22,31,31,31","1^5,32,32,32","221,221,32,32","21,21,21,21,21,21"],  
 ["11,11,11,11,11,11,11,11","1^4,1^4,22,22","1^8,2^4,44",  
 "1^6,2211,222","2211,33,33,33","111,111,111,21,21",  
 "1^5,1^5,2111","1^4,211,211,22","1^4,1^4,211,31",  
 "211,22,22,22,31","1^4,22,22,31,31","111,21,21,21,21,21",  
 "221^8,444,66","2^5,331^4,55","1^8,32111,44",  
 "32211,333,333","332211,3^4,66","3^42211,666,99",  
 "2^5,32221,55","1^7,1^7,511","1^c,5511,66",  
 "3^4,444,5511","541,55,55,55","5541,555,555",  
 "55541,5^4,aa","5^541,aaa,ff","1^8,1^8,62",  
 "1^a1^4,662,77","1^a,55,55,91","2^71,555,87",  
 "21^6,22211,44","221^6,3331,55","1^6,2211,3111",  
 "2^411,33211,55","1^7,3211,331","2211,33,33,411",  
 "3^42,44411,77","22211,2^4,431","2^511,4431,66",  
 "1^8,332,431","3^42,4433,77","1^8,22211,53",  
 "2221,2221,3211","221^5,333,441","1^7,21^5,43",  
 "1^b,443,65","21^5,2221,331","2^51,3332,65",  
 "21^4,21^4,222","221^4,221^4,44","1^6,21^4,321",  
 "2221^4,3322,55","21^4,33,33,42","21^4,222,33,51",  
 "2^51^4,554,77","2^4,3311,3311","3^411,4442,77",  
 "321,321,33,33","3321,3321,333","33321,33321,66",  
 "222,321,33,42","1^6,321,33,51","222,222,321,51",  
 "1^9,3321,54","1^7,322,322","3^422,5551,88",  
 "1^6,33,42,42","1^6,222,42,51","33,33,33,42,51",  
 "1^6,1^6,51,51","222,33,33,51,51","1^b,551,551",  
 "1^5,221,311,41","2^41,3321,441","22111,2221,322",  
 "2^51,443,551","222111,2^41,54","21^4,2211,2211",  
 "1^5,311,32,32","3331,3331,442","2211,2211,33,51",  
 "221,221,311,32","22111,22111,331","1^5,2111,32,41",  
 "2111,221,221,41","2111,221,32,32","211,211,211,211",  
 "211,211,22,31,31","1^4,211,31,31,31","22,22,31,31,31,31"],  
 ["11,11,11,11,11,11,11,11,11","1^5,1^5,1^5","2^5,2^5,55",  
 "111,111,111,111,21","2^41,333,333","1^4,1^4,211,22",  
 "211,22,22,22,22","1^8,22211,44","1^4,1^4,1^4,31",  
 "1^4,22,22,22,31","1^7,1^7,43","1^7,2221,331",  
 "2221,2221,2221","1^6,21^4,222","21^4,33,33,33",  
 "1^6,1^6,321","222,321,33,33","1^6,33,33,42",  
 "222,222,33,42","1^6,222,33,51","222,222,222,51",  
 "33,33,33,33,51","1^6,2211,2211","111,111,21,21,21,21",  
 "1^5,1^5,32,41","1^5,221,221,41","1^5,221,32,32",  
 "221,221,221,32","1^4,211,211,211","211,211,22,22,31",  
 "1^4,211,22,31,31","1^4,1^4,31,31,31","22,22,22,31,31,31",  
 "21,21,21,21,21,21,21","21^a,444,66","1^8,31^5,44",  
 "321^4,333,333","3321^4,3^4,66","3^421^4,666,99",  
 "2^5,322111,55","32^41,3^4,66","3332^41,666,99",  
 "1^8,1^8,611","2^4,44,44,611","1^d,6611,77",  
 "4^5,66611,aa","2^6,444,651","3^4,3^4,651",  
 "651,66,66,66","3^6,6651,99","6651,666,666",  
 "66651,6^4,cc","6^551,ccc,ii","2^8,655,88",  
 "1^9,1^9,72","1^g,772,88","1^c,444,75",  
 "2^6,3^4,75","1^c,66,66,b1","3^4,444,66,b1",  
 "3^7,777,ba","1^7,2221,4111","2^41,333,4311",  
 "1^9,2^41,63","21^8,3331,55","2^411,331^4,55",  
 "1^7,31^4,331","2^411,32221,55","22211,2^4,422",  
 "2^511,4422,66","1^8,332,422","2^5,3331,541",  
 "22211,44,44,62","2^411,2^5,64","2^711,664,88",  
 "1^a,3331,64","2221,2221,31^4","21^7,333,441",  
 "333,333,441,81","2^6111,555,87","21^6,221^4,44",  
 "221^6,3322,55","2^41^6,554,77","1^6,21^4,3111",  
 "3111,321,33,33","33111,3321,333","333111,33321,66",  
 "222,3111,33,42","1^6,3111,33,51","222,222,3111,51",  
 "1^9,33111,54","2221^4,33211,55","1^7,3211,322",  
 "3^4211,5551,88","2^4,3221,3311","333221,4442,77",  
 "3222,3321,333","33222,33321,66","1^9,3222,54",  
 "21^4,33,33,411","3^411,44411,77","222,321,33,411",  
 "1^6,33,411,42","1^6,222,411,51","33,33,33,411,51",  
 "221^4,2^4,431","2^41^4,4431,66","1^8,3311,431",  
 "3^411,4433,77","33321,444,552","1^8,221^4,53",  
 "3311,44,44,53","4^42,5553,99","2^4,3311,44,71",  
 "3^421,555,771","4^52,7771,bb","3^611,776,aa",  
 "2^41,33111,441","22111,2221,3211","2^41,3222,441",  
 "2^61,4441,76","3331,3331,4411","22211,22211,431",  
 "3331,3331,433","3^41,3^41,76","1^7,1^7,61,61",  
 "1^d,661,661","21^5,2221,322","221^5,2^41,54",  
 "2^51,33311,65","21^5,22111,331","3^41,4441,661",  
 "1^7,331,43,61","2221,2221,43,61","2221,331,331,61",  
 "21^4,21^4,2211","21^4,2211,33,51","22211,3311,3311",  
 "1^5,311,311,32","2211,321,33,42","2211,222,321,51",  
 "3322,3331,442","2211,222,42,42","2^411,442,442",  
 "1^6,2211,42,51","2211,33,33,51,51","221,221,311,311",  
 "1^5,2111,311,41","222111,3321,441","22111,22111,322",  
 "222111,222111,54","2111,221,311,32","2111,2111,221,41",  
 "1^5,221,41,41,41","2221,43,43,43","1^5,32,32,41,41",  
 "331,331,43,43","221,221,32,41,41","221,32,32,32,41",  
 "211,211,211,31,31","211,22,31,31,31,31","1^4,31,31,31,31,31"]];  
                         B=BB[MO];  
                 }  
                 if(St!=1){                  if(St!=1){
                         for(R=[]; B!=[]; B=cdr(B)){                          for(R=[]; B!=[]; B=cdr(B)){
                                 RT=F?s2sp(car(B)|std=F):s2sp(car(B));                                  RT= F?s2sp(car(B)|std=F): s2sp(car(B));
                                 if(length(RT)<L0 || length(RT)>L1) continue;                                  if(length(RT)<L0 || length(RT)>L1) continue;
                                 R=cons(RT,R);                                  R=cons(RT,R);
                         }                          }
Line 8523  BB=[
Line 8623  BB=[
         return LL;          return LL;
 }  }
   
   def spbasic(Idx,D)
   {
   /*
     D<=3|Idx|+6,  D<=|Idx|+2 (p>3),  p<=|Idx|/2+4
     Idx=2*D^2-(D^2-\sum m_{j,\nu}^2); \sum(D-m_{j,1})>=2*D;
     \sum (m_{j,1)-m_{j,\nu})*m_{j,\nu)
     0<=(2*D-\sum(D-m_{j,1})})*D=\sum_(m_{j,1}-m_{j,\mu})*m_{j,\nu} -|Idx|
     (-2,0)                                    13ŒÂ (9+3+?)
     (-4,0)                                    37ŒÂ (25+9+?)
     (-6,0) :  8.5sec  ?sec          0.05sec   69ŒÂ (46+17+?)
     (-8,0) : 97  sec  1sec          0.13sec  113ŒÂ (73+29+?)   <- (-2,0)
     (-10,0):          4sec          0.27sec  198ŒÂ (127+50+?)
   @(-12,0)          28sec   4.2sec 0.64sec  291ŒÂ (182+76+?)
     (-14,0)          27sec  10.2sec 1.31sec  415ŒÂ (249+115+?)
     (-16,0)                 34.0sec 2.47sec  647ŒÂ (395+172+?) <- (-4,0)
     (-18,0)                         4.42sec  883ŒÂ (521+243+?) <- (-2,0)
     (-20,0)                         8.17sec 1186ŒÂ (680+345+?)
   */
           Idx=-Idx;
           if((Str=getopt(str))!=1) Str=0;
           if(!isint(Idx)||!isint(Idx/2)||Idx<0||!isint(D)||D<0||D==1||D>3*Idx+6) return [];
           if(D==0){
                   for(R=[],D=3*Idx+6;D>=2;D--) R=append(spbasic(-Idx,D|str=Str),R);
                   return R;
           }
           if(!Idx){
                   R=0;
                   if(D==2) R="11,11,11,11";
                   if(D==3) R="111,111,111";
                   if(D==4) R="22,1111,1111";
                   if(D==6) R="33,222,111111";
                   if(!R) return [];
                   return [(Str==1)?R:s2sp(R)];
           }
           if(D>Idx+2){
                   L=3;
                   if(D==3*Idx+6){
                           R=[[D/2,D/2],[D/3,D/3,D/3],[D/6,D/6,D/6,D/6,D/6,D/6-1,1]];
                           return [(Str==1)?s2sp(R):R];
                   }
                   if(iand(D,1)&&(D-3)/2>Idx) return [];
           }else L=Idx/2+4;
           V=newvect(L);SV=newvect(L);
           for(S1=[],I=0;I<D;I++) S1=cons(1,S1);
           for(T=D-1;T>1;T--){
                   K=D%T;
                   if((T-K)*K<=Idx) break;
           }
           J=(T-K)*K;SJ=K^2+(D-K)*T;
           TV=K?[K]:[];
           for(I=(D-K)/T;I>0;I--) TV=cons(T,TV);
           for(I=0;I<L;I++){
                   SV[I]=2*D^2-(I+1)*(D^2-J)-Idx;
                   V[I]=TV;
           }
           if(SV[2]>0) return [];
           if(D>Idx+2 && V[0][0]+V[1][0]>=D && V[1][0]>1){
                   T=V[1][0]-1;K=D%T;TV=K?[K]:[];
                   for(I=(D-K)/T;I>0;I--) TV=cons(T,TV);
                   V[1]=V[2]=TV;
           }
           for(R=[];;){
                   if(D>Idx+2){
                           if(3*V[0][0]<D) break;
                           if(V[0][0]+V[1][0]>=D && (T=D-V[0][0]-1)>0){
                                   K=D%T;TV=K?[K]:[];
                                   for(I=(D-K)/T;I>0;I--) TV=cons(T,TV);
                                   V[1]=V[2]=TV;
                           }
                           S2=V[0][0]+V[1][0]+V[2][0]-D;
                           if(V[0][0]+2*V[1][0]<D ||(S2<0&&V[1][0]==1) ){
                                   V[0]=V[1]=V[2]=nextpart(V[0]);
                                   T=V[0][0];
                                   T=D-2*T;
                                   if(T==0){
                                           V[1]=[D/2-1,1];
                                           V[2]=S1;
                                   }else if(T>0){
                                           J=D%T;
                                           K=J?[J]:[];
                                           for(J=(D-J)/T;J>0;J--) K=cons(T,K);
                                           V[2]=K;
                                   }
                                   continue;
                           }
                           if(S2<0||V[2][0]<=S2){
                                   V[1]=V[2]=nextpart(V[1]);
                                   continue;
                           }else if(S2>0){
                                   T=V[2][0]-S2;J=D%T;
                                   K=J?[J]:[];
                                   for(J=(D-J)/T;J>0;J--) K=cons(T,K);
                                   V[2]=K;
                           }
                   }
                   for(S=-2*D,IL=0;IL<L;IL++){
                           S+=D-car(V[IL]);
                           if(S>=0) break;
                   }
                   if((I=IL)==L){  /* reducible i.e. IL=L && S<0 */
                           for(LL=L-1;LL>=0;LL--){
                                   if((K=car(V[LL]))+S>0){
                                           K+=S;
                                           for(TV=[],TD=D;TD>=K;TD-=K) TV=cons(K,TV);
                                           if(TD>0) V[LL]=append(TV,[TD]);
                                           else V[LL]=TV;
                                           break;
                                   }else{
                                           S+=K-1;
                                           V[LL]=S1;
                                   }
                           }
                           if(LL<0) break;
                           continue;
                   }
                   for(S0=K=0;K<=IL;K++){
                           ST=car(V[K]);J=V[K][length(V[K])-1];S0+=(ST-J)*J;
                           if(S0>Idx) break;
                   }
                   if(S0>Idx && car(V[K])!=1){
                           ST=car(V[K]);
                           S0-=(ST-J)*J;
                           for(ST--;ST>0;ST--){
                                   J=D%ST;
                                   if(S0+(ST-J)*J <= Idx) break;
                           }
                           V[K]=J?[J]:[];
                           for(J=D-J;J>0;J-=ST) V[K]=cons(ST,V[K]);
                           for(J=K+1;J<L;J++) V[J]=V[K];
                           continue;
                   }
   
                   for(K=SS=0;K<L&&SS<=Idx;K++){
                           ST=car(V[K]);
                           for(S0=0,TV=cdr(V[K]);TV!=[];TV=cdr(TV)) S0+=(ST-car(TV))*car(TV);
                           SS+=S0;
                   }
                   if(SS>Idx && K<=IL && K!=L){
                           SS0=Idx-SS+S0;
                           for(TV=car(V[K]);TV>1;TV--){
                                   U=D%TV;
                                   if((D-U)*U<=SS0) break;
                           }
                           if(TV==car(V[K])){
                                   K=K-1;
                                   V[K]=nextpart(V[K]); /* to be improves */
                           }else{
                                   V[K]=U?[U]:[];  /* to be improved */
                                   for(J=D-U;J>0;J-=TV) V[K]=cons(TV,V[K]);
                           }
                           for(J=K+1;J<L;J++) V[J]=V[K];
                           continue;
                   }
   
                   for(Ix=2*D^2+Idx,J=0;J<L;J++){
                           IxF=Ix;
                           for(Ix-=D^2,TV=V[J];TV!=[];TV=cdr(TV)) Ix+=car(TV)^2;
                           if(Ix<=0) break;
                   }
                   if(Ix==0&&(J>=I||IL==2)){
                           for(TR=[],K=J;K>=0;K--) TR=cons(V[K],TR);
                           R=cons((Str==1)?s2sp(TR):TR,R);
                   }
                   if(J>=0 && J<L && Ix<=0){
                           I=V[J][0];K=D%I;S0=(D-K)*I+K^2;
                           if(I>1&& IxF-D^2+S0<0){
                                   for(V[J]=[],K=D-I;K>0;K--) V[J]=cons(1,V[J]);
                                   V[J]=cons(I,V[J]);
                                   V[J]=nextpart(V[J]);
                                   for(I=J+1;I<L;I++) V[I]=V[J];
                                   continue;
                           }
                   }
                   if(J>=0 && J<L && Ix<=0 && car(V[J])>(U=V[J][length(V[J])-1])+1){
                           TV=reverse(V[J]);
                           for(S0=0,K=[];TV!=[];TV=cdr(TV),S0++){
                                   if((I=car(TV))<U+2||(length(TV)>1&&S0<2)){
                                           while(I-->0) K=cons(1,K);
                                   }else K=cons(car(TV),K);
                           }
                           V[I=J]=K;
                   }else{
                           if(J>=L) J=L-1;
                           for(I=J;I>=0&&length(V[I])==D;I--);
                           if(I<0) break;
                   }
                   V[I]=nextpart(V[I]);                    /* to be improved */
                   for(J=I+1;J<L;J++) V[J]=V[I];
           }
           return R;
   }
   
 def spType2(L)  def spType2(L)
 {  {
         C=0;R=[];          C=0;R=[];
Line 10272  def str_str(S,T)
Line 10564  def str_str(S,T)
                 }else if(type(S)==4){                  }else if(type(S)==4){
                         for(; J<=LE; S=cdr(S),J++){                          for(; J<=LE; S=cdr(S),J++){
                                 if(car(S) != LP){                                  if(car(S) != LP){
                                         if(SJIS && (V=S[J])>128){                                          if(SJIS && (V=car(S))>128){
                                                 if(V<160 || (V>223 && V<240)) J++;                                                  if((V<160 || (V>223 && V<240))&&S!=[]) {
                                                           J++;S=cdr(S);
                                                   }
                                         }                                          }
                                         continue;                                          continue;
                                 }                                  }
Line 10776  def my_tex_form(S)
Line 11070  def my_tex_form(S)
                 }                  }
                 SS = cons(S[I], SS);                  SS = cons(S[I], SS);
         }          }
           SS=str_subst(SS,"\n\\\\\n\\end{pmatrix}","\n\\end{pmatrix}"|raw=1);
         SS=str_subst(SS,"\\\\\n\\end{pmatrix}","\n\\end{pmatrix}"|raw=1);          SS=str_subst(SS,"\\\\\n\\end{pmatrix}","\n\\end{pmatrix}"|raw=1);
         Subst=getopt(subst);          Subst=getopt(subst);
         Sub0=["{asin}","{acos}","{atan}"];          Sub0=["{asin}","{acos}","{atan}"];
Line 11190  def tocsv(L)
Line 11485  def tocsv(L)
                 if(type(LT)==5) LT=vtol(LT);                  if(type(LT)==5) LT=vtol(LT);
                 if(type(LT)<4) LT=[LT];                  if(type(LT)<4) LT=[LT];
                 for(N=0; LT!=[]; LT=cdr(LT),N++){                  for(N=0; LT!=[]; LT=cdr(LT),N++){
                         if(N) str_tb(", ",Tb);                          if(N) str_tb(",",Tb);
                         if((T=car(LT))==Null) continue;                          if((T=car(LT))==Null) continue;
                         if(type(T)==7){                          if(type(T)==7){
                                 K=str_len(T);                                  K=str_len(T);
Line 11280  def readcsv(F)
Line 11575  def readcsv(F)
         return L;          return L;
 }  }
   
   def getline(ID)
   {
           if(isint(Maxlen=getopt(Max))>0) Maxlen=1024;
           if(type(CR=getopt(CR))!=4) CR=[13];
           if(type(LF=getopt(LF))!=4) LF=[10];
           S=[];
           for(I=0; I<1023; I++){
                   C=get_byte(ID);
                   if(C<0) return 0;
                   if(findin(C,CR)>=0) continue;
                   if(findin(C,LF)>=0) break;
                   S=cons(C,S);
           }
           return asciitostr(reverse(S));
   }
   
 def showbyshell(S)  def showbyshell(S)
 {  {
         Id = getbyshell(S);          Id = getbyshell(S);
Line 12396  def rungeKutta(F,N,Lx,Y,IY)
Line 12707  def rungeKutta(F,N,Lx,Y,IY)
         if((Pr=getopt(prec))==1){          if((Pr=getopt(prec))==1){
                 One=eval(exp(0));                  One=eval(exp(0));
         }else{          }else{
                 One=1;Pr=0;                  One=deval(exp(0));Pr=0;
         }          }
         if((FL=getopt(last))!=1) FL=0;          if(!isint(FL=getopt(mul))||!FL) FL=1;
         if(length(Lx)>2){          if(length(Lx)>2){
                 V=car(Lx);Lx=cdr(Lx);                  V=car(Lx);Lx=cdr(Lx);
         }else V=x;          }else V=x;
         if(Pr==0) Lx=[deval(Lx[0]),deval(Lx[1])];          if(Pr==1) Lx=[eval(Lx[0]),eval(Lx[1])];
         else Lx=[eval(Lx[0]),eval(Lx[1])];          else Lx=[deval(Lx[0]),deval(Lx[1])];
         if(type(Y)==4){          if(type(Y)==4){
                 if((Sing=getopt(single))==1||type(F)!=4)                  if((Sing=getopt(single))==1||type(F)!=4)
                         F=append(cdr(Y),[F]);                          F=append(cdr(Y),[F]);
Line 12417  def rungeKutta(F,N,Lx,Y,IY)
Line 12728  def rungeKutta(F,N,Lx,Y,IY)
         }          }
         if(getopt(val)==1) V1=1;          if(getopt(val)==1) V1=1;
         else V1=0;          else V1=0;
         H=(Lx[1]-Lx[0])/N;H2=H/2;          if(FL>0) N*=FL;
           H=(Lx[1]-Lx[0])/N*One;H2=H/2;
         FV=findin(V,vars(F));          FV=findin(V,vars(F));
         K=newvect(4);          K=newvect(4);
         if(L==1){          if(L==1){
                 R=[[T=Lx[0],S=IY]];                  R=[[T=Lx[0],S=IY]];
                 if(!H) return R;                  if(!H) return R;
                 for(;;){                  for(C=0;C<N;C++){
                         for(I=0;I<4;I++){                          for(I=0;I<4;I++){
                                 if(I==0)      W=[[V,T],[Y,S]];                                  if(I==0)      W=[[V,T],[Y,S]];
                                 else if(I==3) W=[[V,T+H],[Y,S+H*K[2]]];                                  else if(I==3) W=[[V,T+H],[Y,S+H*K[2]]];
Line 12432  def rungeKutta(F,N,Lx,Y,IY)
Line 12744  def rungeKutta(F,N,Lx,Y,IY)
                                 K[I]=Pr?myfeval(F,W)*One:myfdeval(F,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;                          S+=(K[0]+2*K[1]+2*K[2]+K[3])*H/6;T+=H;
                         if(!FL) R=cons([deval(T),S],R);                          if(FL>0&&!((C+1)%FL)) R=cons([deval(T),S],R);
                         if((T+H-Lx[1])*H>0) break;  
                 }                  }
         }else{          }else{
                 T=Lx[0];                  T=Lx[0];
                 R=[cons(T,V1?[car(IY)]:IY)];                  R=[cons(T,V1?[car(IY)]:IY)];
                 S=ltov(IY);                  S=ltov(IY);
                 if(!H) return R;                  if(!H) return R;
                 for(;;){                  for(C=0;C<N;C++){
                         for(I=0;I<4;I++){                          for(I=0;I<4;I++){
                                 if(I==0)      W=cons([V,T   ],lpair(Y,vtol(S)));                                  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 if(I==3) W=cons([V,T+H ],lpair(Y,vtol(S+H*K[2])));
Line 12453  def rungeKutta(F,N,Lx,Y,IY)
Line 12764  def rungeKutta(F,N,Lx,Y,IY)
                         }                          }
                         S+=(K[0]+2*K[1]+2*K[2]+K[3])*H/6;T+=H;                          S+=(K[0]+2*K[1]+2*K[2]+K[3])*H/6;T+=H;
                         TS=vtol(S);                          TS=vtol(S);
                           if(FL<0||(C+1)%FL) continue;
                         if(V1) TS=[car(TS)];                          if(V1) TS=[car(TS)];
                         if(!FL) R=cons(cons(deval(T),TS),R);                          R=cons(cons(deval(T),TS),R);
                         if((T+H-Lx[1])*H>0) break;  
                 }                  }
         }          }
         return FL?(V1?S[0]:S):reverse(R);          L=(FL<0)?(V1?S[0]:S):reverse(R);
           return L;
 }  }
   
   #if 1
   def pwTaylor(F,N,Lx,Y,Ly,M)
   {
           if(!isint(FL=getopt(mul))||!FL) FL=1;
           if(getopt(val)==1) V1=1;
           else V1=0;
           if(length(Lx)>2){
                   V=car(Lx);Lx=cdr(Lx);
           }else V=t;
           if(!isvar(T=getopt(var))) V=t;
           if(isint(Pr=getopt(prec))&&Pr>0){
                   One=eval(exp(0));
                   if(Pr>9){
                           setprec(Pr);
                           ctrl("bigfloat",1);
                   }
                   Pr=1;
           }else{
                   One=deval(exp(0));Pr=0;
           }
           if(Pr==1) Lx=[eval(Lx[0]),eval(Lx[1])];
           else Lx=[deval(Lx[0]),deval(Lx[1])];
           if(type(Y)==4){
                   if(type(F)!=4)  F=append(cdr(Y),[F]);
           }else Y=[Y];
           if(type(Ly)!=4) Ly=[Ly];
           if(type(Er=getopt(err))==4){
                   Er=car(Er);ErE=1;               /* Šî€‰ð‚ð•Ô‚· */
           }else ErE=0;
           if(isint(Er)){
                   LyE=Ly;
                   if(Er<0){                               /* Œë·‚̃mƒ‹ƒ€ */
                           Er=-Er;ErR=1;
                   }else ErR=0;                    /* Œë·‚̃mƒ‹ƒ€‚̑Δ */
           }else Er=0;
           if(FL>0) N*=FL;
           S=vtol(pTaylor(F,Y,M|time=V));
           FM=pmaj(F|var=x);
           LS=length(S);
   
           if(type(Vw=getopt(view))==4){   /* Dislay on Canvas */
                   glib_window(car(Vw)[0], car(Vw)[1],car(Vw)[2],car(Vw)[3]);
                   if(length(Vw)>1 && (C=trcolor(Vw[1]))!=0) Opt=[["color",C]];
                   else Opt=[];
                   if(length(Vw)>2){
                           Mt=Vw[2];
                           if(LS==1){
                                   if(type(Mt)>1) Mt=0;
                           }else{
                                   if(type(Mt)!=6||size(Mt)!=[]) Mt=0;
                           }
                   }else Mt=0;
                   if(!Mt){
                           if(LS>1){
                                   Mt=newmat(2,LS);Mt[0][0]=Mt[1][1]=1;
                           }else Mt=1;
                           if(LS==1) glib_putpixel(Lx[0],Mt*Ly[0]|option_list=Opt);
                           else{
                                   YT=ptaffine(Mt,Ly);
                                   glib_putpixel(YT[0],YT[1]|option_list=Opt);
                           }
                   }
           }else Vw=0;
   
           RE=R=[cons(T=Lx[0],Ly)];
           H=(Lx[1]-Lx[0])/N*One;
           Ck=N+1;
           if(type(Inf=getopt(Inf))==4&&length(Inf)>1){    /* explosion */
                   Ck=Inf[0];Ckm=Inf[1];
                   MM=2;
                   if(length(Inf)>2) MM=Inf[2];
                   if(!isint(Mt)||MM<1) MM=2;
                   if(length(Inf)>3) C1=Inf[3];
                   if(!isint(C1)||C1<0) C1=0;
           }else if(Inf==1||Ck<5){
                   Ck=100;Ckm=4;MM=2;C1=0;
           }else Inf=0;
           SS=subst(S,V,H);
           if(Er>0){
                   HE=H/(Er+1);SSE=subst(S,V,HE);
           }
           for(C=CC=0;C<N;C++,CC++){
                   if(CC>=Ck){                                     /* check explosion */
                           CC=0;
                           D0=dnorm(Ly|max=1);
                           for(Dy=F,TY=Y,TL=Ly;TY!=[];TY=cdr(TY),TL=cdr(TL))
                                   Dy=subst(Dy,car(TY),One*car(TL));
                           D1=dnorm(Dy|max=1);D2=subst(FM,x,2*D0+C1);D3=D1+D2;
                           HH=2*(D0+C1)/Ck/Ckm;
                           if(HH<H*D3){
                                   H=HH/D3;
                                   SS=subst(S,V,H);
                                   if(Er){
                                           HE=H/(Er+1);
                                           SSE=subst(S,V,HE);
                                   }
                                   if(MM>1) N*=MM;
                                   MM=0;
                           }
                           CC=0;
                   }
   
                   T+=H;
                   for(Dy=SS,TY=Y,TL=Ly;TY!=[];TY=cdr(TY),TL=cdr(TL))
                           Dy=subst(Dy,car(TY),One*car(TL));
                   Ly=Dy;
   
                   if(Er>0){               /* estimate error */
                           for(CE=-1;CE<Er;CE++){
                                   for(Dy=SSE,TY=Y,TL=LyE;TY!=[];TY=cdr(TY),TL=cdr(TL))
                                           Dy=subst(Dy,car(TY),One*car(TL));
                                   LyE=Dy;
                           }
                   }
   
                   if(FL<0||(C+1)%FL) continue;
                   if(Vw){
                           if(LS==1) glib_putpixel(deval(T),Mt*Ly[0]|option_list=Opt);
                           else{
                                   YT=ptaffine(Mt,Ly);
                                   glib_putpixel(YT[0],YT[1]|option_list=Opt);
                           }
                           continue;
                   }
                   TR=cons((Inf)?eval(T):deval(T),(V1)?[car(Ly)]:Ly);
                   R=cons(TR,R);
                   if(Er){
                           TRE=cons((Inf)?eval(T):deval(T),(V1)?[car(LyE)]:LyE);
                           RE=cons(TRE,RE);
                   }
           }
           if(Vw) return 1;
           L=(FL<0)?((V1)?car(Ly):Ly):reverse(R);
           if(Er){                                                                 /* Estimate error */
                   LE=(FL<0)?((V1)?car(LyE):LyE):reverse(RE);
                   if(FL>0){
                           for(S=L,T=LE,D=[];S!=[];S=cdr(S),T=cdr(T)) D=cons(os_md.ladd(car(S),car(T),-1),D);
                           F=map(os_md.dnorm,reverse(D));
                           if(!ErR) F=map(os_md.nlog,F);
                           if(ErE) return [L,F,LE];
                   }else if(V1){
                           D=ladd(L,LE,-1);F=dnorm(D);
                           if(!ErR) F=nlog(dnorm(D));
                   }else{
                            D=abs(L-LE);
                           if(!ErR) F=nlog(D);
                   }
                   return [L,F];
           }
           return L;
   }
   
   #else
   def pwTaylor(F,N,Lx,Y,Ly,M)
   {
           if(!isint(FL=getopt(mul))||!FL) FL=1;
           if(getopt(val)==1) V1=1;
           else V1=0;
           if(isint(Er=getopt(er))&&Er>0&&type(getopt(Inf))<0){
                   Opt=delopt(getopt(),["er","mul","Inf"]);
                   L0=pwTaylor(F,N,Lx,Y,Ly,M|option_list=cons(["mul",FL*(Er+1)],Opt));
           }else Er=0;
           if(length(Lx)>2){
                   V=car(Lx);Lx=cdr(Lx);
           }else V=t;
           if(!isvar(T=getopt(var))) V=t;
           if((Pr=getopt(prec))==1){
                   One=eval(exp(0));
           }else{
                   One=deval(exp(0));Pr=0;
           }
           if(Pr==1) Lx=[eval(Lx[0]),eval(Lx[1])];
           else Lx=[deval(Lx[0]),deval(Lx[1])];
           if(type(Y)==4){
                   if(type(F)!=4)  F=append(cdr(Y),[F]);
           }else Y=[Y];
           if(type(Ly)!=4) Ly=[Ly];
           if(FL>0) N*=FL;
           S=vtol(pTaylor(F,Y,M|time=V));
           FM=pmaj(F|var=x);
           LS=length(S);
   
           if(type(Vw=getopt(view))==4){   /* Dislay on Canvas */
                   glib_window(car(Vw)[0], car(Vw)[1],car(Vw)[2],car(Vw)[3]);
                   if(length(Vw)>1 && (C=trcolor(Vw[1]))!=0) Opt=[["color",C]];
                   else Opt=[];
                   if(length(Vw)>2){
                           Mt=Vw[2];
                           if(LS==1){
                                   if(type(Mt)>1) Mt=0;
                           }else{
                                   if(type(Mt)!=6||size(Mt)!=[]) Mt=0;
                           }
                   }else Mt=0;
                   if(!Mt){
                           if(LS>1){
                                   Mt=newmat(2,LS);Mt[0][0]=Mt[1][1]=1;
                           }else Mt=1;
                           if(LS==1) glib_putpixel(Lx[0],Mt*Ly[0]|option_list=Opt);
                           else{
                                   YT=ptaffine(Mt,Ly);
                                   glib_putpixel(YT[0],YT[1]|option_list=Opt);
                           }
                   }
           }else Vw=0;
   
           R=[cons(T=Lx[0],Ly)];
           H=(Lx[1]-Lx[0])/N*One;
           Ck=N+1;
           if(type(Inf=getopt(Inf))==4&&length(Inf)>1){    /* explosion */
                   Ck=Inf[0];Ckm=Inf[1];
                   MM=2;
                   if(length(Inf)>2) MM=Inf[2];
                   if(!isint(Mt)||MM<1) MM=2;
                   if(length(Inf)>3) C1=Inf[3];
                   if(!isint(C1)||C1<0) C1=0;
           }else if(Inf==1||Ck<5){
                   Ck=100;Ckm=4;MM=2;C1=0;
           }else Inf=0;
           for(C=CC=0;C<N;C++,CC++){
                   if(CC>=Ck){
                           CC=0;
                           D0=dnorm(Ly|max=1);
                           for(Dy=F,TY=Y,TL=Ly;TY!=[];TY=cdr(TY),TL=cdr(TL))
                                   Dy=subst(Dy,car(TY),One*car(TL));
                           D1=dnorm(Dy|max=1);D2=subst(FM,x,2*D0+C1);D3=D1+D2;
                           HH=2*(D0+C1)/Ck/Ckm;
                           if(HH<H*D3){
                                   H=HH/D3;
                                   SS=subst(S,V,H);
                                   if(MM>1) N*=MM;
                                   MM=0;
                           }
                           CC=0;
                   }
                   if(!C){
                           SS=subst(S,V,H);
                   }
                   T+=H;
                   for(Dy=SS,TY=Y,TL=Ly;TY!=[];TY=cdr(TY),TL=cdr(TL))
                           Dy=subst(Dy,car(TY),One*car(TL));
   /*              Ly=subst(Dy,V,H); */
                   Ly=Dy;
                   if(FL<0||(C+1)%FL) continue;
                   if(Vw){
                           if(LS==1) glib_putpixel(deval(T),Mt*Ly[0]|option_list=Opt);
                           else{
                                   YT=ptaffine(Mt,Ly);
                                   glib_putpixel(YT[0],YT[1]|option_list=Opt);
                           }
                           continue;
                   }
                   TR=cons((Inf)?eval(T):deval(T),(V1)?[car(Ly)]:Ly);
                   R=cons(TR,R);
           }
           if(Vw) return 1;
           L=(FL<0)?((V1)?car(Ly):Ly):reverse(R);
           if(Er){
                   if(FL>0){
                           for(S=L,T=L0,D=[];S!=[];S=cdr(S),T=cdr(T)) D=cons(os_md.ladd(car(S),car(T),-1),D);
                           E=map(os_md.dnorm,reverse(D));F=map(os_md.nlog,E);
                   }else if(V1){
                           D=ladd(L,L0,-1);F=nlog(dnorm(D));
                   }else F=nlog(abs(L-L0));
                   return [L,F];
           }
           return L;
   }
   #endif
   
 def xy2graph(F0,N,Lx,Ly,Lz,A,B)  def xy2graph(F0,N,Lx,Ly,Lz,A,B)
 {  {
         /* (x,y,z) -> (z sin B + x cos A cos B + y sin A cos B,          /* (x,y,z) -> (z sin B + x cos A cos B + y sin A cos B,
Line 13044  def mytan(Z)
Line 13626  def mytan(Z)
 def mylog(Z)  def mylog(Z)
 {  {
         if(type(Z=eval(Z))>1) return todf(os_md.mylog,[Z]);          if(type(Z=eval(Z))>1) return todf(os_md.mylog,[Z]);
         if((Im=imag(Z))==0) return dlog(Z);          if(imag(Z)==0&&Z>=0) return dlog(Z);
         return dlog(dabs(Z))+@i*myarg(Z);          return dlog(dabs(Z))+@i*myarg(Z);
 }  }
   
   def nlog(X)
   {
           return mylog(X)/dlog(10);
   }
   
 def mypow(Z,R)  def mypow(Z,R)
 {  {
         if(type(Z=eval(Z))>1||type(R=eval(R))>1) return todf(os_md.mypow,[Z,R]);          if(type(Z=eval(Z))>1||type(R=eval(R))>1) return todf(os_md.mypow,[Z,R]);
Line 13814  def fcont(F,LX)
Line 14401  def fcont(F,LX)
         return reverse(L);          return reverse(L);
 }  }
   
   def xyplot(L,LX,LY)
   {
           LX=map(deval,LX);LY=map(deval,LY);
           Opt=getopt();Opt0=delopt(Opt,["dviout","proc"]);
           for(S="",L0=[],TL=L;TL!=[];TL=cdr(TL)){
                   TTL=map(deval,car(TL));
                   if(TTL[0]<LX[0]||TTL[0]>LX[1]||TTL[1]<LY[0]||TTL[1]>LY[1]){
                           S+=xylines(reverse(L0)|option_list=Opt0);
                           L0=[];
                   }else{
                           L0=cons(TTL,L0);
                   }
           }
           if(length(L0)>1) S+=xylines(reverse(L0)|option_list=Opt0);
           if(type(AX=getopt(ax))==4) S+=xygraph([0,0],0,LX,LX,LY|option_list=delopt(Opt0,"opt"));
           if(getopt(dviout)!=1) return S;
           xyproc(S|dviout=1);
   }
   
 def xygraph(F,N,LT,LX,LY)  def xygraph(F,N,LT,LX,LY)
 {  {
         if((Proc=getopt(proc))!=1&&Proc!=2&&Proc!=3) Proc=0;          if((Proc=getopt(proc))!=1&&Proc!=2&&Proc!=3) Proc=0;
Line 14126  def xygraph(F,N,LT,LX,LY)
Line 14732  def xygraph(F,N,LT,LX,LY)
                         if(length(Ax)>3){                          if(length(Ax)>3){
                                 D=Ax[3];                                  D=Ax[3];
                                 if(type(D)==1 && D>0){                                  if(type(D)==1 && D>0){
                                         I0=ceil((LY[0]-Ax[1])/D); I1=floor((LY[1]-Ax[0])/D);                                          I0=ceil((LY[0]-Ax[1])/D); I1=floor((LY[1]-Ax[1])/D);
                                         for(DD=[],I=I0; I<=I1; I++){                                          for(DD=[],I=I0; I<=I1; I++){
                                                 if(length(Ax)<5) DD=cons(I*D,DD);                                                  if(length(Ax)<5) DD=cons(I*D,DD);
                                                 else if(I!=0){                                                  else if(I!=0){
Line 15717  def ptcopy(L,V)
Line 16323  def ptcopy(L,V)
         }          }
 }  }
   
   def regress(L)
   {
           E=deval(exp(0));
           for(S0=T0=0,S=L;S!=[];S=cdr(S)){
                   S0+=car(S)[0]*E;T0+=car(S)[1]*E;
           }
           K=length(L);S0/=K;T0/=K;
           for(SS=TT=0,S=L;S!=[];S=cdr(S)){
                   SS+=(car(S)[0]-S0)^2*E;TT+=(car(S)[1]-T0)^2*E;
                   ST+=(car(S)[0]-S0)*(car(S)[1]-T0)*E;
           }
           if(!SS||!TT) return [];
           A=ST/SS;
           L=[A,A*S0-T0,ST/dsqrt(SS*TT),S0,dsqrt(SS/K),T0,dsqrt(TT/K)];
           if(isint(N=getopt(sint))){
                   R=reverse(L);
                   for(L=[];R!=[];R=cdr(R)) L=cons(sint(car(R),N|str=0),L);
           }
           return L;
   }
   
 def     average(L)  def     average(L)
 {  {
         if(getopt(opt)=="co"){          if(getopt(opt)=="co"){
Line 17387  def shiftop(M,S)
Line 18014  def shiftop(M,S)
         return [QQ,P,RS];          return [QQ,P,RS];
 }  }
   
   
   def shiftPfaff(A,B,G,X,M)
   {
           if(type(G)==4){
                   G0=G[1];G1=G[0];
           }
           if(type(G)==6){
                   G=map(red,G);
                   G0=llcm(G);G1=map(red,G0*G);
           }
           if(type(G)==3){
                   G=red(G);G0=dn(G);G1=nm(G);
           }
           if(type(M)==4){
                   M0=M[0];M1=M[1];
           }else{
                   M0=M;M1=0;
           }
           X=vweyl(X);
           D0=mydeg(G0,X[0]);D1=mydeg(G1,X[0]);
           if(M1>=0){
                   D=(D1-M1>D0)?D1-M1:D0;
                   G0=muldo(X[1]^D,G0,X);G1=muldo(X[1]^(D+M1),G1,X);
           }else{
                   D=(D0+M1>D1)?D0+M1:D1;
                   G0=muldo(X[1]^(D-M1),G0,X);G1==muldo(X[1]^D,G1,X);
           }
           G0=map(mc,G0,X,M0);G1=map(mc,G1,X,M0+M1);
           G0=appldo(G0,A,X|Pfaff=1);
           G1=sppldo(G1,B,X|Pfaff=1);
           return rmul(myinv(G0),G1);
   }
   
 def conf1sp(M)  def conf1sp(M)
 {  {
         if(type(M)==7) M=s2sp(M);          if(type(M)==7) M=s2sp(M);
Line 17680  def confspt(S,T)
Line 18340  def confspt(S,T)
 }  }
 #endif  #endif
   
   def mcvm(N)
   {
     X=getopt(var);
     if((Z=getopt(z))!=1) Z=0;
     if(type(N)==4){
       if((K=length(N))==1&&isvar(X)) X=[X];
       if(type(X)!=4){
         for(X=[],I=0;I<K;I++) X=cons(asciitostr([97+I]),X);
         X=reverse(X);
       }
           if(getopt(e)==1){
             if(length(N)==4){
                   N=ltov(N);
                   if(N[1]<N[3]){
                           I=N[1];N[1]=N[3];N[3]=I;
                   }
                   if(N[2]<N[3]||N[2]>=N[1]+N[3]) return 0;
                   X=X[0];
                   for(R=[],I=1;I<N[3];I++) R=cons(makev([X[0],I]),R);
                   for(L=[],I=N[1];I<=N[2];I++) L=cons(makev([X[0],I]),L);
                   for(S=0,I=N[1];I<=N[2];I++){
                     V=makev([X[0],I]);
             S+=polbyroot(R,V)/polbyroot(lsort(L,V,1),V);
                     S=red(S);
                   }
                   return S;
         }
           }
           for(M=[],I=S=0;I<K;Z=0,I++){
                   M=cons(mcvm(N[I]|var=X[I],z=Z),M);
                   S+=N[I];
           }
           M=newbmat(K,K,reverse(M));
       NR=N;
           N=S;
     }else{
           if(type(X)==7) X=strtov(X);
           if(!isvar(X)) X=a;
       M=newmat(N,N);
       for(I=0;I<N;I++){
         V=makev([X,I+1]);
         for(J=0;J<=I;J++){
           R=polbyroot([1,J],V|var=X);
           if(Z==1) R*=V;
           M[I][J]=R;
         }
       }
     }
     if((Get=getopt(get))==1){
       for(R=[],I=0;I<N;I++){
         U=newmat(N,N);
         for(J=0;J<N;J++) U[J][J]=M[J][I];
         R=cons(rmul(rmul(myinv(M),U),M),R);
       }
       return reverse(R);
     }else if(Get==2||Get==3||Get==4){
           for(V=[],I=N;I>0;I--) V=cons(makev(["a0",I]),V);
       MI=myinv(M);
           V=ltov(V)*MI;
           for(R=[],I=0;I<N;I++){
         for(J=I+1;J<N;J++){
           K=newmat(N,N);
           K[I][I]=V[J];K[I][J]=-V[J];K[J][J]=V[I];K[J][I]=-V[I];
              R=cons(rmul(rmul(MI,K),M),R);
          }
           }
       R=reverse(R);
           if(Get==2||length(NR)!=2||Z==1) return R;
       for(V1=[],I=NR[0];I>0;I--) V1=cons(os_md.makev([X[0],I]),V1);
       for(V2=[],I=NR[1];I>0;I--) V2=cons(os_md.makev([X[1],I]),V2);
       R=subst(R,car(V1),0,car(V2),0);
       V1=subst(V1,car(V1),0);
       V2=subst(V2,car(V2),0);
       for(V=[],S=V1;S!=[];S=cdr(S)) for(T=V2;T!=[];T=cdr(T)) V=cons(car(T)-car(S),V);
       V=reverse(V);
       Mx=length(V);
       for(A0=[],I=J=NR[0]-1;J>=0;I+=--J) for(K=0;K<NR[1];K++,I++) A0=cons(R[I],A0);
       A0=reverse(A0);
       for(F0=[],T=1,I=Mx-1;I>=0;I--) F0=cons(1/(x-V[I]), F0);
       MV=confexp([F0,V]|sym=3);
       RR=newvect(Mx);
       for(K=0;K<Mx;K++) for(RR[K]=0,I=0;I<Mx;I++) RR[K]=map(red,RR[K]+MV[I][K]*A0[I]);
       RR0=mysubst(RR,[append(cdr(V1),cdr(V2)),vtol(newvect(Mx-2))]|lpair=1);
       RR0=vtol(RR0);
       return (Get==3)?[RR,RR0]:RR0;
     }
     return M;
   }
   
 def confexp(S)  def confexp(S)
 {  {
           if((Sym=getopt(sym))==1||Sym==2||Sym==3){
                   D=polbyroot(S[1],x);
                   for(R=[],T=S[0];T!=[];T=cdr(T)){
                           M=D*car(T);
                           if(type(M)>3) M=map(red,M);
                           else M=red(M);
                           R=cons(M,R);
                   }
                   R=reverse(R);
                   if(Sym==2) return R;
                   M=length(R);N=length(S[1]);
                   E=newmat(M,N);
                   for(I=0;I<M;I++){
                           for(J=0;J<N;J++) E[I][J]=mycoef(R[I],N-J-1,x);
                   }
                   if(Sym==3){
                           for(R=[],P=1,T=S[1];T!=[];T=cdr(T)) R=cons(P/=(x-car(T)),R);
                           R=confexp([reverse(R),S[1]]|sym=1);
                           return E*myinv(R);
                   }
                   return E;
           }
         if(type(S[0])==4){          if(type(S[0])==4){
                 for(E=[];S!=[];S=cdr(S))                  for(E=[];S!=[];S=cdr(S)) E=cons(confexp(car(S),E));
                         E=cons(confexp(car(S),E));  
                 return reverse(E);                  return reverse(E);
         }          }
         V=x;E=[];          V=x;E=[];
Line 17825  def newbmat(M,N,R)
Line 18594  def newbmat(M,N,R)
         S  = newvect(M);          S  = newvect(M);
         T  = newvect(N);          T  = newvect(N);
         IM = length(R);          IM = length(R);
           if(type(car(R))!=4 && M==N && M==IM){
                   for(RR=TR=[],I=0;I<M;I++){
                           for(TR=[R[I]],J=0;J<I;J++) TR=cons(0,TR);
                           RR=cons(TR,RR);
                   }
                   R=reverse(RR);
           }
         for(I = 0; I < IM; I++){          for(I = 0; I < IM; I++){
                 RI = R[I];                  RI = R[I];
                 JM = length(RI);                  JM = length(RI);
Legend:
Removed from v.1.47  
changed lines
  Added in v.1.59
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>