OpenXM/src/asir-contrib/packages/src/os_muldif.rr - diff

Return to os_muldif.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

Diff for /OpenXM/src/asir-contrib/packages/src/os_muldif.rr between version 1.75 and 1.81

-version 1.75, 2020/10/05 03:42:01
+version 1.81, 2021/07/04 22:31:52
 Line 1
 Line 1
 Line 1
- /* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.74 2020/10/02 07:19:20 takayama Exp $ */
+ /* $OpenXM: OpenXM/src/asir-contrib/packages/src/os_muldif.rr,v 1.80 2020/11/23 10:23:51 takayama Exp $ */
  /* The latest version will be at https://www.ms.u-tokyo.ac.jp/~oshima/index-j.html
   scp os_muldif.[dp]* ${USER}@lemon.math.kobe-u.ac.jp:/home/web/OpenXM/Current/doc/other-docs
  */
 Line 6
 Line 6
 Line 6
  /* #undef USEMODULE */
  /*             os_muldif.rr (Library for Risa/Asir)
-  *          Toshio Oshima (Nov. 2007 - Oct. 2020)
+  *          Toshio Oshima (Nov. 2007 - July 2021)
   *
   *   For polynomials and differential operators with coefficients
   *   in rational funtions (See os_muldif.pdf)
 Line 85  localf trpos$
 Line 85  localf trpos$
 Line 85  localf trpos$
  localf sprod$
  localf sinv$
  localf slen$
+ localf sexps$
  localf sord$
  localf vprod$
  localf dvangle$
-Line 369  localf getbyshell$
+Line 370  localf getbyshell$
 Line 369  localf getbyshell$
 Line 370  localf getbyshell$
  localf show$
  localf dviout$
  localf rtotex$
+ localf togreek$
  localf mtotex$
  localf ltotex$
  localf texbegin$
-Line 381  localf shiftPfaff;
+Line 383  localf shiftPfaff;
 Line 381  localf shiftPfaff;
 Line 383  localf shiftPfaff;
  localf conf1sp$
  localf confexp$
  localf confspt$
+ localf vConv$
  localf mcvm$
  localf s2csp$
  localf partspt$
-Line 464  localf ptaffine$
+Line 467  localf ptaffine$
 Line 464  localf ptaffine$
 Line 467  localf ptaffine$
  localf ptlattice$
  localf ptpolygon$
  localf ptwindow$
+ localf pt5center$
  localf ptconvex$
  localf ptbbox$
  localf darg$
-Line 516  extern AMSTeX$
+Line 520  extern AMSTeX$
 Line 516  extern AMSTeX$
 Line 520  extern AMSTeX$
  extern Glib_math_coordinate$
  extern Glib_canvas_x$
  extern Glib_canvas_y$
- Muldif.rr="00201003"$
+ Muldif.rr="00210702"$
  AMSTeX=1$
  TeXEq=5$
  TeXLim=80$
-Line 1054  def cmpsimple(P,Q)
+Line 1058  def cmpsimple(P,Q)
 Line 1054  def cmpsimple(P,Q)
 Line 1058  def cmpsimple(P,Q)
  def simplify(P,L,T)
  {
-         if(type(P) > 3)
+         if(type(P) > 3){
  #ifdef USEMODULE
                  return map(os_md.simplify,P,L,T);
  #else
                  return map(simplify,P,L,T);
  #endif
+         }
          if(type(L[0]) == 4){
                  if(length(L[0]) > 1)
  #if USEMODULE
-Line 1284  def slen(S)
+Line 1289  def slen(S)
 Line 1284  def slen(S)
 Line 1289  def slen(S)
          return V;
  }
+ def sexps(S)
+ {
+         K=length(S);S=ltov(S);
+         for(R=[],I=0;I<K-1;I++){
+                 for(J=I;J>=0&&S[J]>S[J+1];J--){
+                         T=S[J];S[J]=S[J+1];S[J+1]=T;
+                         R=cons(J,R);
+                 }
+         }
+         return R;
+ }
  def sord(W,V)
  {
          L = length(W);
-Line 2400  def vgen(V,W,S)
+Line 2417  def vgen(V,W,S)
 Line 2400  def vgen(V,W,S)
 Line 2417  def vgen(V,W,S)
  def mmc(M,X)
  {
+         if(getopt(full)==1){
+                 M=mmc(M,X|option_list=delopt(getopt(),"full"));
+                 if(type(M)<4) return -1;
+                 L=length(M);
+                 Mt=getopt(mult);
+                 if((L>=6 && Mt!=0)||(L==3&&Mt==1)){
+                         for(SS=2,I=3; I<L; I+=(++SS));
+                         if(I==L) Mt=1;
+                         else Mt=0;
+                 }
+                 if(Mt!=1){
+                         for(R=[],I=S=0;I<L;I++){
+                                 S=radd(S,M[I]);
+                                 R=cons([[0,I+1],M[I]],R);
+                         }
+                         R=cons([[0,I+1],-S],R);
+                         return reverse(R);
+                 }
+                 for(R=[],I=S=0;I<SS;I++)
+                         for(J=I+1;J<=SS;J++,S++) R=cons([[I,J],M[S]],R);
+                 for(I=0;I<=SS;I++){
+                         for(J=S=0;J<=SS;J++){
+                                 if(I==J) continue;
+                                 S=radd(S,delopt(R,(I<J)?[I,J]:[J,I]|get=1));
+                         }
+                         R=cons([[I,SS+1],-S],R);
+                 }
+                 return qsort(R);
+         }
          Mt=getopt(mult);
          if(type(M)==7) M=s2sp(M);
-         if(type(M)!=4) return 0;
+         if(type(M)!=4&&type(M)!=5) return 0;
          if(type(M[0])<=3){
+                 if(type(M)==5) M=vtol(M);
                  for(RR=[];M!=[];M=cdr(M)) RR=cons(mat([car(M)]),RR);
                  M=reverse(RR);
          }
          if(type(M[0])!=6){                      /* spectre type -> GRS */
-                 G=s2sp(M|std=1);
+                 G=M;
                  L=length(G);
                  for(V=[],I=L-2;I>=0;I--) V=cons(makev([I+10]),V);
                  V=cons(makev([L+9]),V);
-Line 2421  def mmc(M,X)
+Line 2469  def mmc(M,X)
 Line 2421  def mmc(M,X)
 Line 2469  def mmc(M,X)
                  if(Mt!=1) Mt=0;
                  if(R[2]!=2 || R[3]!=0 || !(R=getbygrs(G,1|mat=1))) return 0;
                  MZ=newmat(1,1);
-                 SS=length(G);
+                 SS=length(G)-1;
-                 if(Mt==1) SS=SS*(SS-1)/2;
+                 if(Mt==1) SS=SS*(SS+1)/2;
                  for(M=[],I=0;I<SS;I++) M=cons(MZ,M);
                  for(RR=R; RR!=[]; RR=cdr(RR)){
                          RT=car(RR)[0];
                          if(type(RT)==4){
                                  if(RT[0]!=0) M=mmc(M,[RT[0]]|simplify=Simp);
-                                 M=mmc(M,[cdr(RT)]);
+                                 for(TT=cdr(RT);TT!=[];TT=cdr(TT)){
+                                         if(car(TT)!=0){
+                                                 M=mmc(cdr(M),cdr(RT));
+                                                 break;
+                                         }
+                                 }
                          }
                  }
- /*              for(R=cdr(R);R!=[];R=cdr(R)) M=mmc(M,[car(R)[0]]|mult=Mt); */
          }
          if(X==0) return M;
          L=length(M);
-Line 2442  def mmc(M,X)
+Line 2494  def mmc(M,X)
 Line 2442  def mmc(M,X)
 Line 2494  def mmc(M,X)
          }else{
                  SS=L;Mt=0;
          }
+         if(type(X[0])==4){
+                 for(;X!=[];X=cdr(X)) M=mmc(M,car(X));
+                 return M;
+         }
          if(length(X)==SS+1){
-                 if(car(X)!=0&&(M=mmc(M,[car(X)]|mult=Mt))==0) return M;
+                 if(car(X)!=0) M=mmc(M,[car(X)]|simplify=Simp);
-                 return mmc(M,cdr(X)|mult=Mt);
+                 return mmc(M,cdr(X));
          }
          for(I=X;I!=[];I=cdr(I)) if(I[0]!=0) break;
          if(I==[]) return M;
-Line 2453  def mmc(M,X)
+Line 2509  def mmc(M,X)
 Line 2453  def mmc(M,X)
 Line 2509  def mmc(M,X)
          N=newvect(L);
          for(I=0;I<L;I++) N[I]=dupmat(M[I]);
          S=size(N[0])[0];
-         if(type(X)==4&&length(X)>SS){   /* addition */
+         if(type(X)==4&&length(X)>=SS){  /* addition */
-                 for(I=0;I<SS;I++,X=cdr(X)) if(X[I] != 0) N[I] = radd(N[I],car(X));
+                 for(I=0;I<SS;I++,X=cdr(X)) if(car(X) != 0) N[I] = radd(N[I],diagm(S,[car(X)]));
          }
-         if(length(X)!=1) return 0;
+         if(length(X)!=1||!X[0]) return N;
          X=X[0];
          MZ = newmat(S,S);
          MM = newvect(L);
-Line 2501  def mmc(M,X)
+Line 2557  def mmc(M,X)
 Line 2501  def mmc(M,X)
 Line 2557  def mmc(M,X)
          if(length(KE) == 0) return MM;
          KK = mtoupper(lv2m(KE),0);
          for(I=0;I<L;I++) MM[I] = mmod(MM[I],KK);
-         if(Simp!=0) MM = mdsimplify(MM|type=Simp);
+         if(Simp!=0){
+                 MM = mdsimplify(MM|type=Simp,show=1);
+                 if(getopt(verb)) show([size(MM[0][0]),MM[1]]);
+                 MM=MM[0];
+         }
          return MM;
  }
-Line 3068  def mdsimplify(L)
+Line 3128  def mdsimplify(L)
 Line 3068  def mdsimplify(L)
 Line 3128  def mdsimplify(L)
          return L;
  }
+ #if 1
  def m2mc(M,X)
  {
          if(type(M)<2){
-Line 3158  def m2mc(M,X)
+Line 3219  def m2mc(M,X)
 Line 3158  def m2mc(M,X)
 Line 3219  def m2mc(M,X)
                          if(X[1]=="dviout") Show=2;
                          if(X[1]=="TeX") Show=1;
                  }
-                 if(X[0]=="GRS"||X[0]=="GRSC"||X[0]=="sp"){
+                 if(X[0]=="GRS"||X[0]=="GRSC"||X[0]=="sp"||X[0]=="extend"){
                          Y=radd(-M[0],-M[1]-M[2]);
+                         if(X[0]=="extend")
+                                 return [M[1],M[0],M[2],Y, M[3],M[4],radd(-M[1],-M[3]-M[4]),
+                                         radd(Y,-M[3]-M[4]),radd(M[1],M[2]+M[4]), radd(M[0],M[1]+M[3])];
                          if(X[0]!="GRSC"){
                                  L=meigen([M[0],M[1],M[2],M[3],M[4],Y,radd(-M[1],-M[3]-M[4]),radd(Y,-M[3]-M[4])]|mult=1);
                                  if(X[0]=="sp"){
-Line 3255  def m2mc(M,X)
+Line 3319  def m2mc(M,X)
 Line 3255  def m2mc(M,X)
 Line 3319  def m2mc(M,X)
          KK = mtoupper(lv2m(KE),0);
          for(I=0;I<5;I++)
                  MM[I] = mmod(MM[I],KK);
+         if(Simp!=0){
+                 MM = mdsimplify(MM|type=Simp);
+                 if(getopt(verb)) show([size(MM[0][0]),MM[1]]);
+                 MM=MM[0];
+         }
+         return MM;
+ }
+ #else
+ def m2mc(M,X)
+ {
+         if(type(M)<2){
+         mycat([
+ "m2mc(m,t) or m2mc(m,[t,s])\t Calculation of Pfaff system of two variables\n",
+ " m : list of 5 residue mat. or GRS/spc for rigid 4 singular points\n",
+ " t : [a0,ay,a1,c], swap, GRS, GRSC, sp, irreducible, pair, pairs, Pfaff, All\n",
+ " s : TeX, dviout, GRSC\n",
+ " option : swap, small, simplify, operator, int\n",
+ " Ex: m2mc(\"21,21,21,21\",\"All\")\n"
+ ]);
+                 return 0;
+         }
+         if(type(M)==7) M=s2sp(M);
+         if(type(X)==7) X=[X];
+         Simp=getopt(simplify);
+         if(Simp!=0 && type(Simp)!=1) Simp=2;
+         Small=(getopt(small)==1)?1:0;
+         if(type(M[0])==4){
+                 if(type(M[0][0])==1){ /* spectral type */
+                         XX=getopt(dep);
+                         if(type(XX)!=4 || type(XX[0])>1) XX=[1,length(M[0])];
+                         M=sp2grs(M,[d,a,b,c],[XX[0],XX[1],-2]|mat=1);
+                         if(XX[0]>1 && XX[1]<2) XX=[XX[0],2];
+                         if(getopt(int)!=0){
+                                 T=M[XX[0]-1][XX[1]-1][1];
+                                 for(V=vars(T);V!=[];V=cdr(V)){
+                                         F=coef(T,1,car(V));
+                                         if(type(F)==1 && dn(F)>1)
+                                          M = subst(M,car(V),dn(F)*car(V));
+                                 }
+                         }
+                         V=vars(M);
+                         if(findin(d1,V)>=0 && findin(d2,V)<0 && findin(d3,V)<0)
+                                 M=subst(M,d1,d);
+                 }
+                 RC=chkspt(M|mat=1);
+                 if(RC[2] != 2 || RC[3] != 0){ /* rigidity idx and Fuchs cond */
+                         erno(0);return 0;
+                 }
+                 R=getbygrs(M,1|mat=1);
+                 if(getopt(anal)==1) return R;   /* called by mc2grs() */
+                 Z=newmat(1,1,[[0]]);
+                 N=[Z,Z,Z,Z,Z,Z];
+                 for(RR=R; RR!=[]; RR=cdr(RR)){
+                         RT=car(RR)[0];
+                         if(type(RT)==4){
+                                 if(RT[0]!=0) N=m2mc(N,RT[0]|simplify=Simp);
+                                 N=m2mc(N,[RT[1],RT[2],RT[3]]|simplify=Simp);
+                         }
+                 }
+                 if(type(X)==4 && type(X[0])==7)
+                         return m2mc(N,X|keep=Keep,small=Small);
+                 return N;
+         }
+         if(type(X)==4 && type(X[0])==7){
+                 Keep=(getopt(keep)==1)?1:0;
+                 if(X[0]=="All"){
+                         dviout("Riemann scheme"|keep=1);
+                         m2mc(M,[(findin("GRSC",X)>=0)?"GRSC":"GRS","dviout"]|keep=1);
+                         dviout("Spectral types : "|keep=1);
+                         m2mc(M,["sp","dviout"]|keep=1);
+                         dviout("\\\\\nBy the decompositions"|keep=1);
+                         R=m2mc(M,["pairs","dviout"]|keep=1);
+                         for(R0=R1=[],I=1; R!=[]; I++, R=cdr(R)){
+                                 for(S=0,RR=car(R)[1][0];RR!=[]; RR=cdr(RR)) S+=RR[0];
+                                 if(S==0) R0=cons(I,R0);
+                                 else if(S<0) R1=cons(I,R1);
+                         }
+                         S="irreducibility\\ $"+((length(R0)==0)?"\\Leftrightarrow":"\\Leftarrow")
+                                 +"\\ \\emptyset=\\mathbb Z\\cap$";
+                         dviout(S|keep=1);
+                         m2mc(M,["irreducible","dviout"]|keep=1);
+                         if(R0!=[])
+                                 dviout(ltotex(reverse(R0))|eq=0,keep=1,
+                                  title="The following conditions may not be necessary for the irreducibility.");
+                         if(R1!=[])
+                                 dviout(ltotex(reverse(R1))|eq=0,keep=1,title="The following conditions can be omitted.");
+                         if(getopt(operator)!=0){
+                                 dviout("The equation in a Pfaff form is"|keep=1);
+                                 m2mc(M,["Pfaff","dviout"]|keep=Keep,small=Small);
+                         }
+                         else if(Keep!=1) dviout(" ");
+                         return M;
+                 }
+                 Show=0;
+                 if(length(X)>1){
+                         if(X[1]=="dviout") Show=2;
+                         if(X[1]=="TeX") Show=1;
+                 }
+                 if(X[0]=="GRS"||X[0]=="GRSC"||X[0]=="sp"||X[0]=="extend"){
+                         Y=radd(-M[0],-M[1]-M[2]);
+                         if(X[0]=="extend")
+                                 return [M[1],M[0],M[2],Y, M[3],M[4],radd(-M[1],-M[3]-M[4]),
+                                         radd(Y,-M[3]-M[4]),radd(M[1],M[2]+M[4]), radd(M[0],M[1]+M[3])];
+                         if(X[0]!="GRSC"){
+                                 L=meigen([M[0],M[1],M[2],M[3],M[4],Y,radd(-M[1],-M[3]-M[4]),radd(Y,-M[3]-M[4])]|mult=1);
+                                 if(X[0]=="sp"){
+                                         L=chkspt(L|opt="sp");
+                                         V=[L[1],L[0],L[2],L[5]]; W=[L[1],L[3],L[4],L[6]];
+                                         if(Show==2) dviout(s2sp(V)+" : "+s2sp(W)|keep=Keep);
+                                         return [V,W];
+                                 }
+                                 S="x=0&x=y&x=1&y=0&y=1&x=\\infty&y=\\infty&x=y=\\infty\\\\\n";
+                         }else{
+                                 L=meigen([M[0],M[1],M[2],M[3],M[4],Y,radd(-M[1],-M[3]-M[4]),radd(Y,-M[3]-M[4]),
+                                         radd(M[0],M[1]+M[3]),radd(M[1],M[2]+M[4])]|mult=1);
+                                 S="x=0&x=y&x=1&y=0&y=1&x=\\infty&y=\\infty&x=y=\\infty&x=y=0&x=y=1\\\\\n";
+                         }
+                         T=ltotex(L|opt="GRS",pre=S,small=Small);
+                         if(Show==2) dviout(T|eq=0,keep=Keep);
+                         if(Show==1) L=T;
+                         return L;
+                 }
+                 if(X[0]=="Pfaff"){
+                         S=ltotex(M|opt=["Pfaff",u,x,x-y,x-1,y,y-1],small=Small);
+                         if(Show==2) dviout(S|eq=0,keep=Keep);
+                         return S;
+                 }
+                 if(X[0]=="irreducible"){
+                         L=meigen([M[0],M[1],M[2],radd(-M[0],-M[1]-M[2])]|mult=1);
+                         S=getbygrs(L,10|mat=1);
+                         if(Show==2) dviout(ltotex(S)|eq=0,keep=Keep);
+                         return S;
+                 }
+                 if(X[0]=="pairs"||X[0]=="pair"){
+                         L=meigen([M[0],M[1],M[2],radd(-M[0],-M[1]-M[2])]|mult=1);
+                         S=chkspt(L|opt=0);
+                         V=(Show==2)?1:0;
+                         S=sproot(L,X[0]|dviout=V,keep=Keep);
+                         return S;
+                 }
+                 if(X[0]=="swap"){
+                         Swap=getopt(swap);
+                         if(type(Swap)<1 || Swap==1)
+                                 return newvect(6,[M[3],M[1],M[4],M[0],M[2],M[5]]);
+                         if(Swap==2)
+                                 return newvect(5,[radd(M[0],M[1]+M[3]),M[4],M[2],radd(-M[1],-M[3]-M[4]),M[1]]);
+                         if(type(Swap)==4 && length(Swap)==3){
+                                 MX=radd(-M[0],-M[1]-M[2]); MY=radd(-M[3],-M[1]-M[4]);
+                                 if(Swap[0]==1){
+                                         MX0=M[2];MY0=M[4];
+                                 }
+                                 else if(Swap[0]==2){
+                                         MX0=MX;MY0=MY;
+                                 }else{
+                                         MX0=M[0];MY0=M[3];
+                                 }
+                                 if(Swap[1]==1){
+                                         MX1=M[2];MY1=M[4];
+                                 }
+                                 else if(Swap[1]==2){
+                                         MX1=MX;MY1=MY;
+                                 }else{
+                                         MX1=M[0];MY1=M[3];
+                                 }
+                                 return newvect(5,MX0,M[1],MX1,MY0,MY1);
+                         }
+                 }
+                 return 0;
+         }
+         if(getopt(swap)==1)
+                  return m2mc(m2mc(m2mc(M,"swap"),X),"swap");
+         N=newvect(6);
+         for(I=0;I<6;I++)
+                 N[I]=M[I];
+         S=size(N[0])[0];
+         if(type(X)==4){
+                  for(I=0;I<3;I++){
+                          if(X[I] != 0)
+                                         N[I] = radd(N[I],X[I]);
+                  }
+                  if(length(X)==3) return N;
+                  X=X[3];
+         }
+         MZ = newmat(S,S);
+         ME = mgen(S,0,[X],0);
+         MM = newvect(6);
+         MM[0] = newbmat(3,3, [[N[0]+ME,N[1],N[2]], [MZ], [MZ]]);        /* A01 */
+         MM[1] = newbmat(3,3, [[MZ], [N[0],N[1]+ME,N[2]], [MZ]]);        /* A02 */
+         MM[2] = newbmat(3,3, [[MZ], [MZ], [N[0],N[1],N[2]+ME]]);        /* A03 */
+         MM[3] = newbmat(3,3, [[N[3]+N[1],-N[1]], [-N[0],radd(N[0],N[3])], [MZ,MZ,N[3]]]);       /* A12 */
+         MM[4] = newbmat(3,3, [[N[4]], [MZ,N[4]+N[2],-N[2]], [MZ,-N[1],radd(N[4],N[1])]]);       /* A23 */
+         MM[5] = newbmat(3,3, [[MZ,N[5]+N[2],-N[2]], [N[5]], [MZ,-N[0],radd(N[5],N[0])]]);       /* A13 */
+         M0 = newbmat(3,3, [[N[0]], [MZ,N[1]], [MZ,MZ,N[2]]]);
+         M1 = radd(MM[0],MM[1]+MM[2]);
+         KE = append(mykernel(M0|opt=1),mykernel(M1|opt=1));
+         if(length(KE) == 0) return MM;
+         KK = mtoupper(lv2m(KE),0);
+         for(I=0;I<6;I++)
+                 MM[I] = mmod(MM[I],KK);
          if(Simp!=0) MM = mdsimplify(MM|type=Simp);
          return MM;
  }
+ #endif
  def easierpol(P,X)
  {
-Line 5421  def appldo(P,F,L)
+Line 5685  def appldo(P,F,L)
 Line 5421  def appldo(P,F,L)
 Line 5685  def appldo(P,F,L)
                  L = vweyl(L);
                  X = L[0]; DX = L[1];
                  for(I=mydeg(P,DX);I>0;I--){
-                         if(!(TP=mycoef(P,D,DX))) continue;
+                         if(!(TP=mycoef(P,I,DX))) continue;
-                         P=red(P+TP*(muldo(D^(I-1),F,L)-D^I));
+                         P=red(P-TP*DX^I+TP*muldo(DX^(I-1),F,L));
                  }
                  return P;
          }
-Line 8216  def stoe(M,L,N)
+Line 8480  def stoe(M,L,N)
 Line 8216  def stoe(M,L,N)
 Line 8480  def stoe(M,L,N)
          L = vweyl(L);
          Size = size(M);
          S = Size[0];
-         NN = 0;
+         NN = -1;
          if(type(N) == 4){
                  NN=N[0]; N=N[1];
+                 if(N==NN) return 1;
          }else if(N < 0){
                  NN=-N; N=0;
          }
-Line 8228  def stoe(M,L,N)
+Line 8493  def stoe(M,L,N)
 Line 8228  def stoe(M,L,N)
 Line 8493  def stoe(M,L,N)
          MN = dupmat(M);
          MD = newmat(S,S);
          DD = D[0];
-         DD[N] = 1; DD[S] = 1;
+         DD[N]=1; DD[S] = 1;
          for(Lcm = I = 1; ; ){
                  DD = D[I];
                  MM = MN[N];
-Line 8241  def stoe(M,L,N)
+Line 8506  def stoe(M,L,N)
 Line 8241  def stoe(M,L,N)
 Line 8506  def stoe(M,L,N)
                           DD[J] = red(DD[J]*Lcm);
                  if(I++ >= S)
                          break;
-                 if(I==S && NN>0){
+                 if(I==S && NN>=0){
                          DD = D[I];
-                         DD[0]=-z_zz; DD[NN]=1;
+                         DD[S]=z_zz; DD[NN]=1;
                          break;
                  }
                  Mm = dupmat(MN*M);
-Line 8261  def stoe(M,L,N)
+Line 8526  def stoe(M,L,N)
 Line 8261  def stoe(M,L,N)
 Line 8526  def stoe(M,L,N)
                  if(mydeg(P[I][0],L[1]) > 0)
                           R *= P[I][0]^P[I][1];
          }
-         if(NN > 0)
+         if(NN >= 0)
                  R = -red(coef(R,0,z_zz)/coef(R,1,z_zz));
          return R;
  }
-Line 10066  def mc2grs(G,P)
+Line 10331  def mc2grs(G,P)
 Line 10066  def mc2grs(G,P)
 Line 10331  def mc2grs(G,P)
                          }
                  }
                  if(F=="rest"||F=="eigen"||F=="rest0"||F=="rest1"){
-                         if(F!="eigen") G=mc2grs(G,"homog");
+                         if((Hg=getopt(homog))!=0) Hg=1;
+                         if(F!="eigen"&&Hg) G=mc2grs(G,"homog");
                          if(length(P)==1){
                                  for(R=[],I=0;I<4;I++){
                                          for(J=I+1;J<5;J++){
-                                                 S=mc2grs(G,[F,[I,J]]);
+                                                 S=mc2grs(G,[F,[I,J]]|homog=Hg);
                                                  if(S!=[]) R=cons(cons([I,J],S),R);
                                          }
                                  }
-Line 10855  def mcmgrs(G,P)
+Line 11121  def mcmgrs(G,P)
 Line 10855  def mcmgrs(G,P)
 Line 11121  def mcmgrs(G,P)
  def delopt(L,S)
  {
+         if(getopt(get)==1){
+                 for(;L!=[];L=cdr(L)) if(car(L)[0]==S) return car(L)[1];
+                 return [];
+         }
          if((Inv=getopt(inv))!=1&&Inv!=2) Inv=0;
          if(Inv&&type(S)==4&&type(car(S))==4){
                  for(R=[];L!=[];L=cdr(L)){
-Line 12572  def rtotex(P)
+Line 12842  def rtotex(P)
 Line 12572  def rtotex(P)
 Line 12842  def rtotex(P)
          return (str_len(S) == 1)?S:"{"+S+"}";
  }
+ def togreek(P,T)
+ {
+         R0=[a,b,c,d,e,i,k,l,m,n,o,p,r,s,t,u,x,z];
+         R1=[alpha,beta,gamma,delta,epsilon,iota,kappa,lambda,
+                 mu,nu,omega,pi,rho,sigma,theta,tau,xi,zeta];
+         if(T==0||T==[]) T=[a,b,c];
+         for(S=[],TR=T;TR!=[];TR=cdr(TR)){
+                 if(type(TR[0])!=4){
+                         if((I=findin(car(TR),R0))>=0) S=cons([car(TR),R1[I]],S);
+                 }else if((I=findin(car(TR)[0],R0))>=0){
+                         for(U=car(TR)[1];U!=[];U=cdr(U))
+                                 S=cons([makev([R0[I],car(U)]),makev([R1[I],car(U)])],S);
+                 }
+         }
+         if(getopt(raw)==1) return S;
+         if(getopt(inv)==1) return mysubst(P,S|inv=1);
+         else return mysubst(P,S);
+ }
  def mtotex(M)
  {
          /* extern TexLim;       */
-Line 12959  def xypos(P)
+Line 13248  def xypos(P)
 Line 12959  def xypos(P)
 Line 13248  def xypos(P)
  def xyput(P)
  {
+         if(type(T=car(P))==4||type(car(P)==5)){
+                 P=cdr(P);P=cons(T[1],P);P=cons(T[0],P);
+         }
          if((type(Sc=getopt(scale))==1 && Sc!=1) || type(Sc)==4){
                  if(type(Sc)==1) Sc=[Sc,Sc];
                  Sx=Sc[0];Sy=Sc[1];
-Line 16837  def xybezier(L)
+Line 17129  def xybezier(L)
 Line 16837  def xybezier(L)
 Line 17129  def xybezier(L)
  def xybox(L)
  {
          K=length(L);
+         P=L[0];Q=L[1];
- P=L[0];Q=L[1];
          if(K==2)
                  LL=[ P, [P[0],Q[1]], Q, [Q[0],P[1]] ];
          else{
-Line 24470  P=L[0];Q=L[1];
+Line 17157  P=L[0];Q=L[1];
 Line 24470  P=L[0];Q=L[1];
 Line 17157  P=L[0];Q=L[1];
  def xyang(S,P,Q,R)
  {
          Opt=delopt(getopt(),"ar");
+         if(type(S)>2) S=dnorm([S,P]);
          if(type(Prec=getopt(prec))!=1) Prec=0;
          if(type(Q)>2){
                  if(type(Ar=getopt(ar))!=1) Ar=0;
-Line 25071  def ptwindow(L,X,Y)
+Line 17759  def ptwindow(L,X,Y)
 Line 25071  def ptwindow(L,X,Y)
 Line 17759  def ptwindow(L,X,Y)
          return reverse(R);
  }
+ def pt5center(P,Q,R)
+ {
+         L=newvect(7);
+         L[2]=ptcommon([P,Q],[P,R]|in=-1);
+         Q1=ptcommon([P,R],[Q,0]);R1=ptcommon([P,Q],[R,0]);
+         L[3]=ptcommon([Q,Q1],[R,R1]);
+         P=ltov(P);Q=ltov(Q);R=ltov(R);
+     A=dnorm([Q,R]);B=dnorm([P,R]);C=dnorm([P,Q]);
+         L[0]=vtol((P+Q+R)/3);
+         L[1]=vtol((A*P+B*Q+C*R)/(A+B+C));
+         L[4]=vtol((-A*P+B*Q+C*R)/(-A+B+C));
+         L[5]=vtol((A*P-B*Q+C*R)/(A-B+C));
+         L[6]=vtol((A*P+B*Q-C*R)/(A+B-C));
+         return vtol(L);
+ }
  def lninbox(L,W)
  {
          if(L[0]==L[1]) return 0;
-Line 27158  def confspt(S,T)
+Line 19862  def confspt(S,T)
 Line 27158  def confspt(S,T)
 Line 19862  def confspt(S,T)
  }
  #endif
+ def vConv(K,I,J)
+ {
+         if(type(X=getopt(var))!=7) X="a";
+         if(getopt(e)==2) return subst(vConv(K,I+1,J+1),makev([X,1]),0);
+         if(J>K){L=J;J=K;K=L;}
+         if(K>I||J<1||K+J<I+1) return 0;
+         if(K+J==I+1) return 1;
+         else
+ #if 1
+         L=I-K<J-2?I-K+1:J;
+         for(S=0,M=0;M<L;M++) S+=(makev([X,K+M])-makev([X,J-M-1]))*vConv(K+M,I,J-M-1|var=X);
+         return S;
+ #else
+         return  vConv(K+1,I,J-1|var=X)+(makev([X,K])-makev([X,J-1]))*vConv(K,I,J-1|var=X);
+ #endif
+ }
  def mcvm(N)
  {
    X=getopt(var);
-Line 27165  def mcvm(N)
+Line 19886  def mcvm(N)
 Line 27165  def mcvm(N)
 Line 19886  def mcvm(N)
    if(type(N)==4){
      if((K=length(N))==1&&isvar(X)) X=[X];
      if(type(X)!=4){
-       for(X=[],I=0;I<K;I++) X=cons(asciitostr([97+I]),X);
+       for(X=[],I=0;I<K;I++) X=cons(asciitostr([97+I]),X);       /* a,b,... */
        X=reverse(X);
      }
-         if(getopt(e)==1){
+         if((E=getopt(e))==1||E==2){
-           if(length(N)==4){
+           if(length(N)==4) N=cdr(N);
-                 N=ltov(N);
+           if(length(N)==3) return vConv(N[0],N[1],N[2]|var=X,e=E);
-                 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);
-Line 27240  def mcvm(N)
+Line 19947  def mcvm(N)
 Line 27240  def mcvm(N)
 Line 19947  def mcvm(N)
      MV=confexp([F0,V]|sym=3);
      RR=newvect(Mx);
      for(K=0;K<Mx;K++) for(RR[K]=0,I=0;I<Mx;I++) RR[K]=map(red,RR[K]+MV[I][K]*A0[I]);
-     RR0=mysubst(RR,[append(cdr(V1),cdr(V2)),vtol(newvect(Mx-2))]|lpair=1);
+         for(RR0=RR,VV=append(cdr(V1),cdr(V2));VV!=[];VV=cdr(VV)) RR0=subst(RR0,car(VV),0);
      RR0=vtol(RR0);
      return (Get==3)?[RR,RR0]:RR0;
    }
-Line 29412  def bernoulli(N)
+Line 22119  def bernoulli(N)
 Line 29412  def bernoulli(N)
 Line 22119  def bernoulli(N)
  }
  /* linfrac01([x,y]) */
- /* linfrac01(newvect(10,[0,1,2,3,4,5,6,7,8,9]) */
+ /* (x_0,x_1,x_2,x_3,...,x_{q+3})=(x,0,1,y_1,...,y_q,\infty)
- /* 0:x=0, 1:x=y, 2:x=1, 3:y=0, 4:y=1, 5:x=\infty, 6:y=\infty, 7:x=y=0, 8:x=y=1, 9:x=y=\infty
-:y_2=0, 11:y_2=x, 12:y_2=y, 13: y_2=1,   14: y_2=\infty
-:y_3=0, 16:y_3=x, 17:y_3=y, 18: y_3=y_2, 19: y_3=1, 20:y_3=\infty
-          X[0],X[11],X[2],X[10],X[13],X[5],X[14],X[7],X[8],X[9],
-          X[3],X[1],X[12],X[4],X[6]
          T=0   (x_2,x_1,x_3,x_4,...)
          T=-j  (x_1,x_2,..,x_{j-1},x_{j+1},x_j,x_{j+2},...)
          T=1   (1-x_1,1-x_2,1-x_3,1-x_4,...)
          T=2   (1/x_1,1/x_2,1/x_3,1/x_4,...)
-         T=3   (x_1,x_1/x_2,x_1/x_3,x_1/x_4,...)
+         T=3   (1/x_1,x_2/x_1,x_3/x_1,x_4/x_1,...)
+     ...
  */
  def lft01(X,T)
  {
-         MX=getopt();
+         S=0;
          if(type(X)==4){
+                 if(type(car(X))==4){
+                         S=X[1];X=car(X);
+                 }
                  K=length(X);
                  if(K>=1) D=1;
          }
-         if(type(X)==5){
-                 K=length(X);
-                 for(J=5, F=K-10; F>0; F-=J++);
-                 if(F==0) D=2;
-         }
          if(D==0) return 0;
-         if(T==0){  /* x <-> y */
+         if(type(T)==4&&(length(T)==K+3||length(T)==2)){
-                 if(D==1){
+                 for(U=[],I=K+2;I>=0;I--) U=cons(I,U);
-                         R=cdr(X); R=cdr(R);
+                 if(length(T)==2) T=mperm(U,[T],0);
-                         R=cons(X[0],R);
+                 L=sexps(T);
-                         return cons(X[1],R);
+                 for(R=[X,S];L!=[];L=cdr(L)){
+                         if(!(I=car(L))) I=4;
+                         /* else if(I==1) I=1; */
+                         else if(I==2) I=5;
+                         else if(I==K+1) I=6;
+                         else if(I>2) I=2-I;
+                         R=lft01(R,I);
                  }
-                 R=newvect(K,[X[3],X[1],X[4],X[0],X[2],X[6],X[5]]);
-                 for(I=7;I<K;I++) R[I]=X[I];
-                 for(I=11,J=5; I<K; I+=J++){
-                         R[I]=X[I+1]; R[I+1]=X[I];
-                 }
                  return R;
          }
-         if(T==1){
+         if(!S) S=getopt(tr);
-                 if(D==1){
+         if(type(S)==4&&length(S)==K+3){
-                         for(R=[];X!=[];X=cdr(X)) R=cons(1-car(X),R);
+                 D=2;
-                         return reverse(R);
+         }else if(S==1) for(S=[],I=K+2;I>=0;I--) S=cons(I,S);
+         else S=0;
+         if(T<=0){  /* y_i <-> y_{i+1}, y_0=x=x_0, y_i=x_{i+2} */
+                 R=mperm(X,[[-T,1-T]],0);
+                 if(S){
+                         if(!T) S=mperm(S,[[0,3]],0);
+                         else   S=mperm(S,[[2-T,3-T]],0); /* : J J=3,...,K; */
+                         R=[R,S];
                  }
-                 R=newvect(K,[X[2],X[1],X[0],X[4],X[3],X[5],X[6],X[8],X[7],X[9]]);
-                 for(I=11;I<K;I++) R[I]=X[I];
-                 for(I=10, J=5; I<K; I+=J++){
-                         R[I]=X[I+J-2]; R[I+J-2]=X[I];
-                 }
                  return R;
-         }
+         }else if(T==1){ /* (x_1=0, x_2=1) : 1 */
-         if(T==2){
+                 for(R=[];X!=[];X=cdr(X)) R=cons(1-car(X),R);
-                 if(D==1){
+                 if(S) S=mperm(S,[[1,2]],0);
-                         for(R=[]; X!=[]; X=cdr(X)) R=cons(red(1/car(X)),R);
+         }else if(T==2){ /* (x_1=0, x_{K+2}=infty) */
-                         return reverse(R);
+                 for(R=[]; X!=[]; X=cdr(X)) R=cons(red(1/car(X)),R);
-                 }
+                 if(S) S=mperm(S,[[1,K+2]],0);
-                 R=newvect(K,[X[5],X[1],X[2],X[6],X[4],X[0],X[3],X[9],X[8],X[7]]);
+         }else if(T==3){ /* (x_0=x, x_2=1) */
-                 for(I=11;I<K;I++) R[I]=X[I];
+                 T=car(X);
-                 for(I=10,J=5; I<K; I+=J++){
+                 for(R=[red(1/T)],X=cdr(X); X!=[]; X=cdr(X)) R=cons(red(car(X)/T),R);
-                         R[I]=X[I+J-1]; R[I+J-1]=X[I];
+                 if(S) S=mperm(S,[[0,2]],0);
-                 }
+         }else if(T==4){ /* (x_0=x,x_1=0) : 0 */
-                 return R;
+                 T=car(X);
-         }
+                 for(R=[red(T/(T-1))],X=cdr(X); X!=[]; X=cdr(X)) R=cons(red((T-car(X))/(T-1)),R);
-         if(T==3){
+                 if(S) S=mperm(S,[[0,1]],0);
-                 if(D==1){
+         }else if(T==5){ /* (x_2=1,x_3=y) : 2 */
-                         T=car(X);
+                 T=X[1];
-                         for(R=[T],X=cdr(X); X!=[]; X=cdr(X))
+                 for(R=[1/T,red(X[0]/T)],X=cdr(cdr(X));X!=[]; X=cdr(X)) R=cons(red(car(X)/T),R);
-                                 R=cons(red(T/car(X)),R);
+                 if(S) S=mperm(S,[[2,3]],0);
-                         return reverse(R);
+         }else if(T==6){ /* (x_{K+1}=y_{K-1}, x_{K+2}=infty) : K+1 */
-                 }
+                 T=X[K-1];
-                 R=newvect(K,[X[7],X[4],X[2],X[6],X[1],X[9],X[3],X[0],X[8],X[5]]);
+                 for(R=[];length(X)>1;X=cdr(X)) R=cons(red(car(X)*(1-T)/(car(X)-T)),R);
-                 for(I=10,J=5; I<K; I+=J++){
+                 R=cons(1-T,R);
-                         R[I]=X[I+J-1]; R[I+1]=X[I+J-2]; R[I+J-2]=X[I+1]; R[I+J-1]=X[I];
+                 if(S) S=mperm(S,[[K+1,K+2]],0);
-                 }
+         }else if(T==7){ /* x_2=1 <-> x_{K+2}=infty */
-                 return R;
+                 for(R=[];X!=[];X=cdr(X)) R=cons(red(car(X)/(car(X)-1)),R);
-         }
+                 if(S) S=mperm(S,[[2,K+2]],0);
-         if(T==-1){
+         }else return 0;
-                 if(D==1){
+         R=reverse(R);
-                         return append([X[1],X[2],X[0]],cdr(cdr(cdr(X))));
+         return S?[R,S]:R;
-                 }
-                 R=newvect(K,[X[0],X[11],X[2],X[10],X[13],X[5],X[14],X[7],X[8],X[9],
-                          X[3],X[1],X[12],X[4],X[6]]);
-                 for(I=11;I<K;I++) R[I]=X[I];
-                 for(I=17,J=5; I<K; I+=J++){
-                         R[I]=X[I+1]; R[I+1]=X[I];
-                 }
-                 return R;
-         }
-         if(T<0){
-                 if(D==1){
-                         for(R=[],I=0; X!=[]; X=cdr(X),I--){
-                                 if(I==T){
-                                         R=cons(X[1],R);
-                                         R=cons(X[0],R);
-                                         X=cdr(X);
-                                 }
-                                 else R=cons(car(X),R);
-                         }
-                         return reverse(R);
-                 }
-                 T=3-T;
-                 R=newvect(K);
-                 for(I=0;I<K;I++) R[I]=X[I];
-                 for(I=10,J=5;J<T;I+=J++);
-                 for(II=0; II<J-2; II++){
-                         R[I]=X[I+J]; R[I+J]=R[I];
-                 }
-                 for( ; II<J; II++){
-                         R[I]=X[I+J+1]; R[I+J+1]=X[I];
-                 }
-                 return R;
-         }
-         return 0;
  }
  def linfrac01(X)
  {
-         if(type(X)==4) K=length(X)-2;
+         if(type(X)==4){
-         else if(type(X)==5){
+                 K=length(X)-2;
-                 L=length(X);
+                 if(type(car(X))==4){
-                 for(K=0,I=10,J=5; I<L; K++,I+=J++);
+                         for(U=[],I=K+4;I>=0;I--) U=cons(I,U);
-                 if(I!=L) return 0;
+                         X=[car(X),U];
+                 }else U=0;
          }
          if(K>3 && getopt(over)!=1) return(-1);
          II=(K==-1)?3:4;
-Line 29549  def linfrac01(X)
+Line 22219  def linfrac01(X)
 Line 29549  def linfrac01(X)
 Line 22219  def linfrac01(X)
                          }
                  }
          }
-         return L;
+         return reverse(L);
  }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>