############################################################################
# Module: ParametricLogDer													                                 					
###########################################################################

ParametricLogDer := module ()
  option package;
  description "The module includes functions for determining parametric logarithmic derivatives";
  export
     ExtendedEuclidean,
     Derivative,
     HermiteReduce,
     RationalCombination,
     LiouvillianPLD,
     LiouvillianLD;
  local
     BaseCase,
     Recursion;

 ##################################################################################
 # Name: ExtendedEuclidean
 # Calling sequence: ExtendedEuclidean(x, p, a, b)
 # Input:   x,  a variable;
 #          p, a, b,  three polynomials of x with gcd(a, b) = 1;
 # Output:  [u, v],  two polynomials such that p = u*a + v*b with deg(u) < deg(b).
 ##################################################################################

 ExtendedEuclidean := proc(x, p, a, b)
     local h, u, v, pp1, pp2, t, c1, c2, q, s;
     h := gcdex(a, b, x, 's');
     if degree(h, x) > 0 then
         error "a nontrivial gcd is found";
     end if;
     u := rem(s*p, b, x);
     pp1 := primpart(normal(p - u*a), x, 'c1');
     pp2 := primpart(b, x, 'c2');
     t := evala(Divide(pp1, pp2, 'qu'));
     if t = false then
        error "the division is not exact";
     end if;
     v := normal(c1/c2)*qu;
     [u, v];
 end proc:

#############################################################################################
# Name: Derivative
# Calling sequence: Derivative(x, T, dT, f)
# Input:   x,    a variable of the base field C(x);
#          T = [t1, t2, .., tn],    a list of indeterminates over the base field, maybe [];
#          dT = [t1', t2', .., tn'],   a list of derivatives of elements in T, maybe [];
#          f,     a rational function in multivariate differential field C(x)(T).
# Output:  the derivative of f.
#############################################################################################

Derivative := proc(x, T, dT, f)
   local i, n, df;

   if T = [] then
       return diff(f, x);
   end if;
   n := nops(T);
   df := diff(f, x);
   for i from 1 to n do
       df := df + diff(f, T[i])*dT[i];
   end do:
   df := normal(df);
end proc:

 #########################################################################################################
 # Name: HermiteReduce
 # Calling sequence: HermiteReduce(x, T, dT, f)
 # Input:   x,   a variable of the base field C(x);
 #          T = [t1, t2, .., tn],   a list of generators of a transcendental and elementary tower K_n over C(x), maybe empty;
 #          dT = [t1', t2', .., tn'],   a list of derivatives of elements in T, maybe empty;
 #          f,   a rational function in multivariate differential field C(x)(T).
 # Output:  [g, h, p]
 #          such that f = g' + h + p,
 #          where g, h are in K_n and h is t_n -simple, p is in K_{n-1}[t_n], if T <> [];
 #                g, h are in C(x) and h is x-simple, p = 0, if T = [].
 #########################################################################################################

 HermiteReduce := proc(x, T, dT, f)
     local y, d, deg, l, d1, p1, p2, m, a, p, g, i, j, u, v, B, h, fp, up, vp, k, M, N, pp;

     # find the main "variable"

     if T = [] then
         y := x;
     else
         y := op(-1, T);
     end if;

     # test hyperexponential or primitive

     if T = [] then
        deg := 0;
     else
        deg := degree(op(-1, dT), y);
     end if;

    # compute the proper and polynomial parts with the monic denominator

    pp := 0;
    if deg = 0 then
       # primitive case
       fp := f;
    else
       # hyperexponential case
       (M, N) := denom(f), numer(f);
       if normal(subs(y=0, M)) <> 0 then
          # the denominator has no special factor
          fp := f;
       else
          k := 1;
          M := normal(M/y);
          while normal(subs(y=0, M)) = 0 do
                k := k + 1;
                M := normal(M/y);
          end do;
          gcdex(M, y^k, y, 'up', 'vp');
          fp := normal((N*vp)/M);   # normal denominator
          pp := normal((N*up)/y^k); # the nontrivial special denominator
        end if;
     end if;


    d := denom(fp);
    l := lcoeff(d, y);
    a := rem(numer(fp), d, y, 'p1');
    d := normal(d/l);
    a := normal(a/l);

    # square-free decomposition

    d1 := sqrfree(d, y)[2];

    # Hermite Reduction

    g := 0;
    for i from 1 to nops(d1) do
        (v, m) := d1[i][1], d1[i][2];
        u := normal(d/v^m);

        for j from m - 1 by -1 to 1 do
            B := ExtendedEuclidean(y, -a/j, u*Derivative(x, T, dT, v), v);
            g := g + B[1]/v^j;
            a := normal(-j*B[2] - u*Derivative(x, T, dT, B[1]));
        end do;
        d := normal(u*v);
    end do;

    a := rem(a, d, y, 'p2');
    p := normal(p1 + p2);
    h := normal(a/d);

    if T = [] then
        g := g + int(p, x);
        p := 0;
    end if;

    return [g, h, normal(p+pp)];
end proc:

BaseCase := proc(x, f, paralist)
      local ff, k, U, S, h, a, d, F, G, j, m, p, q, r, P, u;

      ff := normal(f);
      k := nops(paralist);

      #------------------------
      # trivial case
      #------------------------

      if ff = 0 then
         return {}, paralist;
      end if;

      #-------------------------
      # Rational reduction
      #-------------------------

      U := HermiteReduce(x, [], [], ff);

      #---------------------------
      # x-proper part
      #---------------------------

      S := {coeffs(expand(numer(U[1])), x)} minus {0};



      #---------------------------
      # x-simple part
      #---------------------------

      h := U[2];

      if h = 0 then
         return S, paralist;
      end if;

      (a, d) := numer(h), denom(h);
      F := factors(d);
      G := F[2];
      m := nops(G);
      P := paralist;

      for j from 1 to m do
         p := G[j][1];
         q := quo(d, p, x);
         gcdex(q, p, x, 'u');
         r := rem(a*u, p, x);
         S := S union {coeffs(collect(r-c[k+j]*diff(p,x),x),x)} minus {0};
         P := P union {[c[k+j], p]};
      end do;
      return S, P;
end proc:

Recursion := proc(x, T, dT, f, paralist)
   local ff, t, ind, U, a, d, r, g, h, p, S, W, b, L, k, m, j, w, deg, s, n, Sp, Lp;

   if T = [] then
      return BaseCase(x, f, paralist);
   end if;

   ff := normal(f);

   #------------------------
   # trivial case
   #------------------------


   if ff = 0 then
      return {}, paralist;
   end if;

   #-----------------------
   # initialization
   #-----------------------

   t := T[-1];
   ind := [x, op(T)];


   #-------------------------
   # Hermite reduction
   #-------------------------

   U := HermiteReduce(x, T, dT, ff);
   (a, d) := numer(U[1]), denom(U[1]);
   r := rem(a, d, t, 'q');
   g := r/d;
   h := U[2];
   p := U[3]+Derivative(x, T, dT, q);


   #---------------------------
   # tn-proper part
   #---------------------------

   S := {coeffs(collect(numer(g), ind, distributed), ind)} minus {0};

   #----------------------------
   # t-simple part
   #----------------------------

   d := denom(U[2]);
   W := factors(d)[2];
   b := 0;
   if degree(dT[-1],t) = 1 then
      b := normal(dT[-1]/t);
   end if;
   L := paralist;
   k := nops(L);
   m := nops(W);
   for j from 1 to m do
         w := W[j][1];
         if degree(w, t)>0 then
            w := w/lcoeff(w,t);
            k := k+1;
            h := h - c[k]*Derivative(x, T, dT, w)/w;
            if b <> 0 then
               h := h + c[k]*degree(w, t)*b;
               p := p - c[k]*degree(w, t)*b;
            end if;
            L := L union {[c[k],w]};
         end if;
   end do;
   S := S union {coeffs(collect(numer(h), ind, distributed), ind)} minus {0};

   #------------------------------
   # t-reduced part
   #------------------------------

   (a, d) := numer(p), denom(p);
   deg := degree(d, t);
   s := coeff(a, t, deg)/lcoeff(d, t);
   p := p - s;
   S := S union {coeffs(collect(numer(p), ind, distributed), ind)} minus {0};

   #-------------------------------
   # recursion
   #-------------------------------

   n := nops(T);
   (Sp, Lp) := Recursion(x, [op(1..n-1,T)], [op(1..n-1, dT)], s, L);
   return S union Sp, L union Lp;
end proc:

###########################################################################
# Name: RationalCombination
# Calling sequence:
#       RationalCombination(x, T, dT, f)
# Input: x, a variable,
#        T, a list of indeterminates, maybe empty,
#        dT, a list of derivatives, maybe empty,
#        f, an element of C(x, T)
# Output: either a list [ [r1, a1], ... [r_k, a_k]], where r1, ..., r_k are rational
#         numbers, and a1, ..., a_k are elements of C(x, T) such that
#
#                 f = r1*a1'/a1 + ... + r_k*a_k'/a_k
#         or false if there does not exist such a list.
##############################################################################

RationalCombination := proc(x, T, dT, f)
        local ff, P, n, a, i, t, tp, S, sol,c;
        ff := normal(f);
        if ff = 0 then
           return([[1, 1]]);
        end if;

        P := {};
        n := nops(T);
        a := 0;
        for i from 1 to n do
            (t, tp) := T[i], dT[i];
            if degree(tp, t) = 1 then
               a := a + c[i]*lcoeff(tp, t);
               P := P union {[c[i], t]};
            end if;
        end do;

        (S, P) := Recursion(x, T, dT, f - a, P);
        sol := solve(S, indets(S));
        if sol = NULL then
           return false;
        end if;
        return subs(sol, P);
end proc:

###########################################################################
# Name: LiouvillianPLD
# Calling sequence: LiouvillianPLD(x, T, dT, f)
# Input: x, a variable;
#          T, a list of regular and Liouvillian monomials; not empty; T[-1] hyperexponential
#          dT, derivatives of T;
#           f, an element in C(x)(T[1..n-1]);
#Output: [n, u] such that
#                  f=n*T[-1]'/T[-1]+u'/u, where n is an integer and u is in C(x)(T[1..n-1]) and nonzero,
#             Otherwise, 
#             [], if  no such[n,u] exists.
##############################################################################

LiouvillianPLD := proc(x, T, dT, f)
        local    ff, P, m, u, n, i;

       # trivial case

        ff := normal(f);
        if ff = 0 then
           return([0, 1]);
        end if;
        
       P :=RationalCombination(x, T, dT, f);
       if P=false then
          return [];
       end if;
       m := nops(P);
        u  := 1;
        n  :=1;
       for i from 1 to m do
             if  type( P[i][1], 'integer') then
                 if P[i][2]=T[-1] then
                     n := P[i][1];
                  else
                     u :=u*P[i][2]^P[i][1]
                  end if;
              else
                 return [];
                 break;
             end if;
      end do;
      return [n,u];
      end proc:


###########################################################################
# Name: LiouvillianLD
# Calling sequence: LiouvillianLD(x, T, dT, f)
# Input: x, a variable;
#          T, a list of regular and Liouvillian monomials; not empty; 
#          dT, derivatives of T;
#           f, an element in C(x)(T);
#Output: u such that
#                  f=u'/u, where u is in C(x)(T) and nonzero,
#             Otherwise, 
#             [], if  no such u exists.
##############################################################################

LiouvillianLD := proc(x, T, dT, f)
        local    ff, P, m, u, n, i;

       # trivial case

        ff := normal(f);
        if ff = 0 then
           return(1);
        end if;
        
       P :=RationalCombination(x, T, dT, ff);
       if P=false then
          return [];
       end if;
       m := nops(P);
        u  := 1;
        n  :=1;
       for i from 1 to m do
             if  type( P[i][1], 'integer') then
                 u :=u*P[i][2]^P[i][1]
              else
                 return [];
                 break;
             end if;
      end do;
      return u;
      end proc:





end module: 