Return to taka_diffop.rr CVS log

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

Utility functions for diff ops.
poly_dact, poly_dmul, poly_diff2euler, poly_euler2diff.

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_diffop.rr,v 1.7 2022/02/01 00:36:22 takayama Exp $ */
#define DEBUG

#define ANK_SIZE 20

#define XX   x
#define LOG  log(x)
/* #define LOG  log(-x) */

/* module taka_dr1 is defined later. */

module tk_diff;
localf euler2diff$
localf diff2euler1$
localf diff2euler$
localf dact$
localf dmul$
def euler2diff(F,XL,TL) {
  DXL=poly_dvar(XL);
  return taka_dr1.euler2diff(F,XL,TL,DXL);
}
def diff2euler1(F,X,T) {
  /* T is theta_X=X*Dx */
  DN = dn(F);
  F = nm(F);
  Dx = poly_dvar(X);
  F = dp_ptod(F,[X,Dx]);
  R = 0;
  while (F != 0) {
    P=dp_hm(F);
    F=dp_rest(F);
    E = dp_etov(P);
    R += dp_hc(P)*dp_vtoe(newvect(2,[E[0]-E[1],0]))*taka_dr1.dp_poch(E[1]);
  }
  R = dp_dtop(R,[X,T]);
  R = eval_str(rtostr(R));
  DN2= dn(R);
  R = nm(R);
  /* factor out X */
  M = taka_dr1.root_multiplicity(R,X,0);
  return [red(R/X^M),X^M/(DN*DN2)];
}

def diff2euler(F,XL,TL) {
  L=[F,1];
  for (I=0; I<length(XL); I++) {
    L2=diff2euler1(L[0],XL[I],TL[I]);
    L=[L2[0],red(L[1]*L2[1])];
  }
  return(L);
}
/*&usage
  Example:
   ED=diff2euler((x-y)*dx*dy+a*dx+b*dy,[x,y],[tx,ty]);
   red(euler2diff(ED[0],[x,y],[tx,ty])*ED[1]);
*/
/*&usage
dact(L,F,XL) acts the operator L to F. XL is a list of X variables.
*/
def dact(L,F,XL) {
  if (length(XL)==0) return(L*F);
  D = dn(L);
  L = nm(L);
  X=XL[0]; Dx=poly_dvar(X);
  N = deg(L,Dx);
  R = 0;
  for (I=0; I<=N; I++) {
   G=dact(coef(L,I,Dx),taka_dr1.diff2(F,X,I),cdr(XL));
   R += G;
  }
  return red(R/D);
}
/*&usage
  Example:
  load("ok_diff.rr");
  L0=(x*dx+(1-x-y)*dy)^3; F0=1/(x^3-y^2); V0=[x,y];
  printf("actdiff=%a\n",red(dact(L0,F0,V0)-odiff_act(L0,F0,V0)));
*/
def dmul(F,G,XL) {
  if ((F==0) || (G==0)) return 0;
  if ((dn(F) != 1) || (dn(G) != 1)) error("dmul works only for polynomial differential operators.");
  DXL=poly_dvar(XL);
  V=append(XL,DXL);
  F=dp_ptod(F,V); G=dp_ptod(G,V);
  R=dp_weyl_mul(F,G);
  return dp_dtop(R,V);
}
endmodule;


module taka_dr1;

localf diff2 $
localf mult$
localf in$
localf in_coef$
localf division$
localf gcd$
localf act$
localf initANK$
localf ank$
localf thetapower2diff$
localf euler2diff$
localf diff2euler$
localf dp_poch $
localf int_roots $
localf root_multiplicity $
localf indicial $
localf kcoef0 $
localf test1 $
localf ord0 $
localf  pfq_infty $

extern ANK $
ANK = newmat(ANK_SIZE+1,ANK_SIZE+1) $

/* --------------------------------------------------------------------
   Ring of differential operators with rational functions coefficients
   in one variables.
   R_1 = K < x, dx >
   ---------------------------------------------------------------------
*/

/* cf. d_gcd.rr */
def mult(F,G){
	N=deg(nm(F),dx);
	A=0;
	for(K=0;K<=N;K++){
		A=A+(1/fac(K))*diff2(F,dx,K)*diff2(G,x,K);
	}
	return red(A);
}
def diff2(F,G,K){
    for(I=0;I<K;I++){
      F=diff(F,G);
	}
    return F;
}


def in(F) {
  Dn = dn(F); F=nm(F);
  D = deg(F,dx);
  C = coef(F,D,dx);
  return( red(C*dx^D/Dn) );
}
def in_coef(F) {
  Dn = dn(F); F=nm(F);
  D = deg(F,dx);
  C = coef(F,D,dx);
  return( red(C/Dn) );
}

def division(F,G) {
   Q = 0; R = F;
   DegG = deg(nm(G),dx);
   while ((R != 0) && (deg(nm(R),dx) >= DegG)) {
      /* D = red(in(R)/in(G));  It does not work for complex coef */
      DegR = deg(nm(R),dx);
      D = red(in_coef(R)/in_coef(G))*dx^(DegR-DegG);
      Q = red(Q+D);
      R = red(R-mult(D,G));
   }
   return([Q,R]);
}


def gcd(F,G) {
   S = 1; T = 0; SS=0; TT=1;
   if (deg(nm(F),dx) > deg(nm(G),dx)) {
     Changed = 0;
   }else {
     T=F; F=G; G=T; Changed = 1;
   }
   do {
     D=division(F,G);
     Q = D[0]; R = D[1];
     SSS = S-Q*SS; TTT = T-Q*TT;
     F = G; G = R;
     S = SS; T = TT; SS = SSS; TT = TTT;
   } while ( R != 0);
   if (Changed) {
     return [F,T,S];
   }else {
     return [F,S,T];
   }
}

def act(L,F) {
  D = dn(L);
  L = nm(L);
  N = deg(L,dx);
  R = 0;
  for (I=0; I<=N; I++) {
   R += coef(L,I,dx)*diff2(F,x,I);
  }
  return red(R/D);
}

/* --------------------------------------------------------
   euler operators
   --------------------------------------------------------
*/

def initANK() {  /* Note, 2/22, 2003 */
  extern ANK ;
  ANK[1][1] = 1;
  for (NN=2; NN<ANK_SIZE; NN++) {
    for (KK=1; KK<=NN; KK++) {
      if ((NN-1 < 1) || (KK > NN-1)) A = 0;
      else A = (ANK[NN-1][KK])*KK;
      if (KK-1 < 1) B = 0;
      else B = ANK[NN-1][KK-1];
      ANK[NN][KK] = A+B;
    }
  }
}

def ank(N,K) {
  extern ANK ;
  if (!ANK[1][1]) {
    initANK();
  }
  if ( N < 1 ) error("ank() out of index.");
  if ( K > N ) return 0;
  if ( K <= 0 ) return 0;
  if (N < ANK_SIZE) return ANK[N][K];
  return ank(N-1,K)*K+ank(N-1,K-1);
}

/*  Expand (X D)^J in X and D */
def thetapower2diff(J,X,D) {
  if (J < 0) error("thetapower2diff() out of index.");
  if (J == 0) return 1;
  R = 0;
  for (I = 1; I <= J; I++) {
    R += ank(J,I)*X^I*D^I;
  }
  return R;
}

def euler2diff(F,Xlist,Tlist,Dlist) {
  DN =dn(F);
  F = nm(F);
  for (I=0; I < length(Xlist); I++) {
    T = Tlist[I]; X = Xlist[I]; D = Dlist[I];
    Deg = deg(F,T); R = 0;
    for (J=0; J<=Deg; J++) {
      C = coef(F,J,T);
      R += C*thetapower2diff(J,X,D);
    }
    F = R;
  }
  return red(F/DN);
}

/* -----------------------------------
   diff2euler(x*dx+c,t) ==> [t+c,1]
   diff2euler(dx+c,t) ==> [t+cx,1/x]
   ----------------------------------
*/
def diff2euler(F,T) {
  DN = dn(F);
  F = nm(F);
  F = dp_ptod(F,[x,dx]);
  R = 0;
  while (F != 0) {
    P=dp_hm(F);
    F=dp_rest(F);
    E = dp_etov(P);
    R += dp_hc(P)*dp_vtoe(newvect(2,[E[0]-E[1],0]))*dp_poch(E[1]);
  }
  R = dp_dtop(R,[x,T]);
  R = eval_str(rtostr(R));
  DN2= dn(R);
  R = nm(R);
  /* factor out x */
  M = root_multiplicity(R,x,0);
  return [red(R/x^M),x^M/(DN*DN2)];
}

def dp_poch(K) {
  One = <<0,0>>;
  T = <<0,1>>;
  R = One;
  for (I=0; I<K; I++) {
    R = R*(T-I*One);
  }
  return R;
}

/* int_roots(x^2*(x-1)^3*(x^2+1),x);   --> [[1,3], [0,2]]
   non-negative integral roots.
*/
def int_roots(F,S) {
  Root = [];
  F = fctr(F);
  F = cdr(F);
  N = length(F);
  for (I=0; I<N; I++) {
    P  = F[I][0];  /* factor */
    Pm = F[I][1];  /* Multiplicity */
    if (deg(P,S) == 1) {
      P = -subst(P,S,0); /* Find root. */
      if (type(P) <=1) {
        if (P >= 0) Root = cons([P,Pm],Root);
      }
    }
  }
  return Root;
}

/*  root_multiplicity(x^2*(x-1),x,0); --> 2
    multiplicity at x = 0.
*/
def root_multiplicity(F,S,P) {
  F = fctr(F);
  F = cdr(F);
  N = length(F);
  for (I=0; I<N; I++) {
    Pf  = F[I][0];  /* factor */
    Pm = F[I][1];  /* Multiplicity */
    if (deg(Pf,S) == 1) {
      Pf = -subst(Pf,S,P); /* Find root. */
      if (Pf == 0) {
        return Pm;
      }
    }
  }
  return 0;
}

/* L is a differential operator in x and dx */
def indicial(L) {
  /* indicial polynomial is a polynomial in x and s */
  L = diff2euler(L,s)[0];
  L0 = subst(L,x,0);
  return [L0,-(L-L0)];  /* inidicial and the rest */
}


/* ---------------------------------------------------
  Determine a_{k,i}
  a_{k,i} x^k (log(x))^i,  i=0, 1,..., m
  L = [L0,LR] (euler op in s, output of indicial).  L0-LR
  Prev : series solution in x upto degree k-1 
  0 ==> slow version.
  DONT USE THE VARIABLE atmp_ 
----------------------------------------------------- */
def kcoef0(K,M,L,Prev) {
  L0 = L[0];  D_L0 = euler2diff(L0,[x],[s],[dx]);
  LR = L[1];  D_LR = euler2diff(LR,[x],[s],[dx]);
  V = [ ];
  F = 0;
  M0 = root_multiplicity(L0,s,K);  /* exclude starting terms */
  for (I=M0; I<=M; I++) {
    V = cons(util_v("atmp",[I]),V);
    F = F + util_v("atmp",[I])*(XX)^K*(LOG)^I;
  }
  V = reverse(V);
  /* print([M0,M,V,F]); */
  P0 = act(D_L0,F);
  P1 = act(D_LR,Prev);
  P0=red(subst(P0,LOG,y)/x^K);
/*
  P1=coef(subst(P1,LOG,y),K,x);
*/
  Ptmp = subst(P1,LOG,y);
  P1=red(coef(nm(Ptmp),K,x)/dn(Ptmp));

  #ifdef DEBUG
  print([P0,P1]);
  #endif

  /* compare coefficients of y=LOG */
  Eq = [ ];
  PPP = nm(P0-P1);
  for (I=0; I<=M; I++) {
    /* G = coef(P0,I,y)-coef(P1,I,y); */
    G = coef(PPP,I,y);
    if (G != 0) {
       Eq = cons(G,Eq); 
    }
  }
  EqAns = poly_solve_linear(Eq,V);
  #ifdef DEBUG
  print(EqAns);
  #endif
  return base_replace(F,EqAns);  
}

/* testing  kcoef0 */
def test1() {
  L = euler2diff(s^2-x*(s+1/2)^2,[x],[s],[dx]);
  LL = indicial(L);
  F0 = 1;
  print(sprintf("Indicial = %a",LL));
  for (I=1; I<4; I++) {
    Fnew = kcoef0(I,1,LL,F0);
    print(sprintf("I=%a, %a",I,Fnew));
    F0 = F0 + Fnew;
  }
  print("------------------------------");
  F1 = LOG;
  for (I=1; I<4; I++) {
    Fnew = kcoef0(I,1,LL,F1);
    print(sprintf("I=%a, %a",I,Fnew));
    F1 = F1 + Fnew;
  }
  print("------------------------------");

  return [F0,F1];
  /*  cf. SST, p.23
    [1+(1/4)*x+(9/64)*x,
     log(x)+(-2*(1/2)^2+1 + 1/4*log(x))*x]
  */
}

/* order at X=0 of F*/
def ord0(F,X) {
  F = nm(F);
  N = deg(F,X);
  for (I=0; I<=N; I++) {
    if (coef(F,I,X) != 0) return I;
  }
  return -1;
}

/*
   Construct series solutions of p F q (A,B) up to O(1/x^N)
  Use "CHECK" to check if a given series is an approximate sol.
  order must be equal to the last arg (5 or 3).
   taka_dr1.pfq_infty([1/2,1/2],[0],5);  
   taka_dr1.pfq_infty([a,b],[c],5);  
   taka_dr1.pfq_infty([a,a],[c],5);
   taka_dr1.pfq_infty([a,a,b],[c1,c2],3); 
   taka_dr1.pfq_infty([a,a,b],[2,3],3); 
   taka_dr1.pfq_infty([1/2,1/2,-1/3],[2,3],5);
   taka_dr1.pfq_infty([1/2,1/2,1/2,-3/2],[2,3,4],5); 
*/
/*&usage begin: taka_dr1.pfq_infty(A,B,N)
 It computes the set of the canonical series solutions up to degree {N}
 at z = infty for the generalized hypergeometric function p F q ({A},{B};z).
 The variable 1/z is put x.
 example: taka_dr1.pfq_infty([1/2,1/2,-3/2],[2,3],5); 
*/
def pfq_infty(A,B,N) {
  P = length(A);  
  Q = length(B);
  L0 = -s;
  for (I=0; I<Q; I++) {
    L0 = (-s+B[I]-1)*L0;
  }
  L1 = 1;
  for (I=0; I<P; I++) {
    L1 = (-s+A[I])*L1;
  }
  L = L1-x*L0;  /* operator at infty */

  L1 = cdr(fctr(L1));
  M = length(L1);
  Roots = [ ];
  LogMult = 0;  /* Max exponent of log(x) */
  for (I=0; I<M; I++) {
    Roots = cons(-subst(L1[I][0],s,0)/coef(L1[I][0],1,s),Roots);
    if (L1[I][1]-1 > LogMult) LogMult = L1[I][1];
  }
  Roots = reverse(Roots);
  M = length(Roots);
  Ans = [ ]; 
  for (Ri=0; Ri<M; Ri++) {
    /* x^E F(x). Le is the operator for F(x) */
    E = Roots[Ri];
    Le = subst(L,s,s+E);
    LL = indicial(Le);
    MM = root_multiplicity(LL[0],s,0);

    for (K=0; K<MM; K++) {
      F0 = (LOG)^K;

      #ifdef DEBUG
      print(sprintf("Indicial = %a",LL));
      #endif

      for (I=1; I<N; I++) {
        Fnew = kcoef0(I,LogMult,LL,F0);

        #ifdef DEBUG
        print(sprintf("I=%a, %a",I,Fnew));
        #endif

        F0 = F0 + Fnew;
      }

      #ifdef DEBUG
      LLL = euler2diff(Le,[x],[s],[dx]);
      FFF = act(LLL,F0);
      print(sprintf("order in x is %a",ord0(FFF,x)));
      #endif      

      print("------------------------------");

      #ifdef DEBUG
/*      debug; */
      #endif

      Ans = cons([x^E,F0],Ans);
    }
  }
  return reverse(Ans);
} 


endmodule;

end$