###########################################################################
#
# Module:  Reduction							   
# 						                                 					   
###########################################################################

Reduction := module ()
  option package;
  description "algorithms for reductions";
  export 
     GCanonicalForm,
     GKernelReduction,
     GKernelAndShellReduction, 
     HermiteReduce,       
     Indicator,
     GShellReduction1, 
     GShellReduction,
      ResidualForm0,
      LaurentPolynomialReduction1,
      LaurentPolynomialReduction2,
      LaurentPolynomialReduction;
   local
     GKernelAndShellReduction0,
     GKernelReduction1,
     RelativeFactor,
     GKernelReduction0,
     SIntegrable,
     GFindSepcialPolyWildCase,
     PolynomialReduction2,
     GMyRatSum,
      PolynomialReduction; 

 #########################################################################################################
 # 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 monomial extension K_n over C(x), maybe empty, tn is hyperexpoenetial or no special;
 #          dT = [t1', t2', .., tn'],   a list of derivatives of elements in T, maybe empty;
 #          f,   an element 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, l, d1, p1, p2, m, a, p, g, i, j, u, v, B, h,ord,s,r;

     # find the main "variable"

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

    # compute the normally proper and reduced parts with the monic denominator

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

    if T<>[] and degree (dT[-1],y)=1 then
	ord:=0;
	while rem(d,y,y)=0 do
		ord:=ord+1;
		d:=normal(d/y);
	end do;
	gcdex(y^(ord),d,a,y,'s','r');
	p1:=p1+r/y^(ord);
	a:=s;
    end if;
    # 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 := Rational:-ExtendedEuclidean(y, -a/j, u*DField:-Derivative(x, T, dT, v), v);
            g := g + B[1]/v^j;
            a := normal(-j*B[2] - u*DField:-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, p];
end proc:



#########################################################
# Name: Indicator
# Calling Sequence：Indicator(T, M)
# Input:  T, a list of indeterminates in ascending order;
#         M, a monomial in T.
# Output: i, the indicator of M.
#########################################################

Indicator  :=  proc(T, M)
    local i, n, ind;
    n := nops(T);
    if depends(M, T) then
        for i from 1 to n do 
            if depends(M, T[i]) then 
                return i;
            end if;
        end do;
    else 
        return n;
    end if;
end proc:




################################################################################
# Name: SIntegrable
# Calling sequence: SIntegrable(x, T, A, f, ind)
# Input: x, a variable; 
#           T, a list of indeterminates;
#           A, the associated matrix of the differential field generated by T;
#           f,  an element in the differential field generated by T in Matryoshka form;
#           ind, an integer between 0 and ind;
# Ouput:  C := [c_1, ..., c_[ind]] such that
#                     (1)  c_1, ..., c_[ind] are constant, not all  zero;
#                     (2)  f = c_1*T[1]' + ... + c_[ind]*T[ind]';
#             [], if no such constants exist.
################################################################################

SIntegrable := proc(x, T, A, f, ind)
   local S, Y, Z, X,  i,  L,  j,  M,  Tp, p, u; 
      S := {};  Y := [];  Tp := [x, op(T)]; 
      for i from 1 to ind do
            Y  := [op(Y), X[i]]; 
            L  := 0; 
            for j from 1 to ind do
                 L := L + A[i,j]*X[j];
            end do;
            L := L - f[i];
            S := S union {coeffs(collect(numer(L), Tp, 'distributed'), Tp)};
      end do;
      M := SolveTools:-Linear(S, {op(Y)}); 
      if M = NULL then
         return [];
      end if;
      for i from 1 to ind do
           u := M[i];
           member(lhs(u), Y, 'p');
           Z[p] := rhs(u);
      end do;
      return[seq(Z[j], j=1..ind)];
end proc;

#################################################################################
#Name: GCanonicalForm
#Calling sequence: GCanonicalForm(x, T, dT, f)
#Input:  x,   a variable of the base field C(x);
#          T = [t1, t2, .., tn],   a list of generators of a hyperexponential 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:  [k, s] such that 
#               f=k+s'/s;
#      where k is diffrential reduced;
#      neither num(s) nor den(s) has the factor tn
#       gcd(den(K), den(S))=1.
###################################################################################                  

GCanonicalForm := proc(x, T, dT, f)
local t, a, b, F, n, i, m, p, c, q,  r, u,  s, k;

   # find the main "variable" 

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

     a := numer(f);
     b := denom(f);
     F := factors(primpart(b, t));
     n := nops(F[2]);

     s := 1;
     k := f;
     for i from 1 to n do
         (p, m) := F[2][i][1], F[2][i][2];
         if m = 1 then
            if  t = x or rem(p, t,  t)<> 0 then
               q := normal(f-c*DField:-Derivative(x,T, dT, p)/p);
               r := rem(numer(q), p, t);
               u := normal(-subs(c=0, r)/lcoeff(r, c));
               if type(u, integer) then
                  s := normal(s*p^u);
                  k := normal(k-u*DField:-Derivative(x,T, dT, p)/p);
               end if;
            end if;
         end if;
      end do;

      return [k, s];
end proc:

############################################################################### 
# Name: GKernelReduction1
# Calling sequence:
#
#                   GKernelReduction1(x, T, dT, p, k1, k2, m)
#
# 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 [];
#            p, k1, k2, three polynomials with K := k1/k2 differential-reduced ；
#             m,         a positive integer.
# Output: [p1, p2], two polynomials such that
#              p/k2^m = (p1/k2^(m-1))' + p1/k2^(m-1)*K + p2/k2.
#####################################################################

GKernelReduction1  :=  proc(x, T, dT, p, k1, k2, m)
   local t, res, p1s, p2s, p1, p2, rds;
   if T = [] then
      t := x;
  else
     t := op(-1, T);
  end if;
  if m = 1 then
     return([0, p]);
  end if;

  res := Rational:-ExtendedEuclidean(t, p, k1-(m-1)*DField:-Derivative(x, T, dT ,k2), k2);
  p1s := res[1];
  p2s := normal(res[2] - DField:-Derivative(x, T, dT ,res[1]));
  rds := GKernelReduction1(x, T,  dT,  p2s,  k1,  k2,  m-1);
  p1 := normal(rds[1]*k2+p1s);
  p2 := rds[2];
  return [p1, p2];

end proc:

############################################################################ 
#Name: GKernelReduction(x, T, dT, f, k1, k2)
#Calling sequence:
#                       GKernelReduction(x, T, dT, f, k1, k2)
#
#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 w.r.t tn  in C(x)(T) whose denominator is a divisor of some power of k2,
#          k1, k2 two coprime polynomials in tn such that K := k1/k2 is diffrential reduced;
#Output: [u, p] such that
#             f = u'+K*u+p/k2,
#             where u is a rational fuction, p is a polynomial w.r.t tn .
###########################################################################
  
GKernelReduction := proc(x, T, dT, f, k1, k2)
   local  b, t, l, h, W; 

   #------------------------------------------------------
   #  Find the main variable
   #------------------------------------------------------
 
   if T = [] then
      t := x;
  else
     t := op(-1, T);
  end if;

  b := denom(f);

  #------------------------------------------------------
  # Find the correct multiplicity of k2
  #------------------------------------------------------

  l := 1;
  h := k2;
  while  rem(h, b, t)<>0 do
     l  := l+1;
     h := h*k2;
  end do;

  #----------------------------------------------------------------
  # Kernel reduction with multiplicity
  #----------------------------------------------------------------

  W := GKernelReduction1(x, T, dT, normal(h*f), k1, k2, l);
  return ([normal(W[1]*k2/h), W[2]]);

end proc:  



##################################################################################### 
# Name: RelativeFactor
# Calling sequence:
#            RelativeFactor(p, q, t)
# Input:
#        p, q, two nonzero polynomials in t;
#        t, a name
# Output,  [r, s] such that
#                p = r*s, 
#                no irreducilbe factors of r are factors of q,
#                all irreducible factors of s are factors of q;
################################################################################## 

RelativeFactor := proc(p, q, t)
         local r, s, g, gs;
         (r, s) := (p, 1);
         while true do
           g := gcd(numer(r), numer(q));
           if degree(g, t) = 0 then
              return [r, s];
          else
             gs := g/lcoeff(g, t);
             r := quo(r, gs, t);
             s  := s*gs;
          end if;
        end do;
end proc: 

################################################################################# 
#Name: GShellReduction1
#Calling sequence:
#                GShellReduction1(x, T, dT, a, d, i, k1, k2)
#Input : x,   a variable of the base field C(x);
#          T = [t1, t2, .., tn],   a list of generators of a hyperexponential tower K_n over C(x), may be empty;
#          dT=[t1', t2', ...,tn']   a list of derivatives of T;
#          a, d,   elements of K_{n-1}[tn] and d is normal if T <>[],
#                    elements of C[x] and d is squarefree if T =[] ;
#          i, an positive integer;
#          k1, k2,  two polynomials with K:=k1/k2 diffrential reduced and gcd(d, den(K))=1.
#Output: [p1, p2, p3]such that
#             a/d^i = (p1)' + p1*K + p2+p3;
#            where p1, p2, p3 are in C(x)(T)  with p2 =a2/den(K) and p3= a3/d if T <>[],
#                  p1, p2, p3 are in C(x) with p2 =a2/den(K) and p3= a3/d if T=[].
############################################################################### 

GShellReduction1:=proc(x, T, dT, a, d, i, k1, k2) 
    local t, j, p1, p2, p3, M1, M2, N, s ; 
    if T = [] then
       t :=x;
    else
       t := op(-1, T);
    end if; 

    if i=1 then
      return [0, 0, a/d];
    end if; 

    p1 := 0;
    p2 := 0;
    p3 := a; 
     
    gcdex(DField:-Derivative(x, T, dT, d), d, t, 's');
     for j from i by -1 to 2 do 

        #get a/(1-i)=M1*d'+M2*d

        M1 := rem( (s*p3)/(1-j), d, t);
        M2 := quo(p3/(1-j)-M1*DField:-Derivative(x, T, dT, d), d, t);
       #M := Rational:-ExtendedEuclidean(t, p3/(1-j), DField:-Derivative(x, T, dT, d), d );  #a/(1-i)=M[1]d'+M[2]d

       p1 := normal(p1+M1/d^(j-1));
       N  := Rational:-ExtendedEuclidean(t, -k1*M1, k2, d^(j-1));
       p2 := normal(p2+N[2]/k2);
       p3 := normal(-DField:-Derivative(x, T, dT, M1)+(1-j)*M2+N[1]); 

     end do; 

    return ([p1, p2, p3/d]);

 end proc: 

################################################################################## 

#Name: GShellReduction
#Calling sequence:
#                  GShellReduction(x, T, dT, f, k1, k2)
#Input : x,   a variable of the base field C(x);
#          T = [t1, t2, .., tn],   a list of generators of a hyperexponential tower K_n over C(x), may be empty;
#          dT=[t1', t2', ...,tn']   a list of derivatives of T;
#          f, a rational function in C(x)(T) whose any irreducible factor p  of den(f) satisfies  gcd(p, k2)=1 and p<>t_n  if T<>[];
#          f, a rational function in C(x) if T=[];
#          k1,k2,  K:=k1/k2 is diffrential reduced.
#Output: [u, h, p] such that
#             f= u' + u*K + h+p/den(K);
#             where u is in C(x)(T),  h is t_n-simple and gcd(den(K), den(h))=1, p is in C(x)(t1, ...,t_(n-1))[t_n, t_n^{-1}]  if T<>[];        
#                        u is in C(x), h is x-simple and gcd(den(K), den(h))=1, p is in C[x] if  T=[].

############################################################################### 

GShellReduction := proc(x, T, dT, f, k1, k2) 
    local t, a, b, u, p, h, S, j, d, m, M, N, i; 

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

    (a, b) :=numer(f), denom(f);
    u := 0;
    p := 0;
    h := 0;
  
  #reduce irreducible factors d of den(f]) satisfy gcd(d, k2)=1, 

  S := factors(b);
  for j from 1 to nops(S[2]) do 
      (d, m) := S[2][j][1], S[2][j][2];
      #M := Rational:-ExtendedEuclidean(t, a, normal(b/(d^m)), d^m);  #M[1]/d^m + M[2]/(b/d^m) =a/b
      M :=Rational:-EEA(t, a, [S[1], seq(S[2][i][1]^S[2][i][2], i=j+1..nops(S[2]))], d^m); #M[1]/d^m + M[2]/(b/d^m) =a/b
      N := GShellReduction1(x, T, dT, M[1], d, m, k1, k2);
      u := normal(u+N[1]);
      p := normal(p+N[2]);
      h := normal(h+N[3]);
      b := normal(b/d^m);
      a := M[2]; 
  end do;
  return ([u, h, normal(p*k2)]);
end proc:

############################################################################################################################ 
#Name: GKernelAndShellReduction (can only be applied to primitive and hyperexponential)
#Calling sequence:
#                    GKernelAndShellReduction(x, T, dT, f, k1, k2)
# Input:   x,   a variable of the base field C(x);
#          T = [t1, t2, .., tn],   a list of generators of a Trans. Liouvillian tower K_n over C(x), maybe empty;
#          dT = [t1', t2', .., tn'],   a list of derivatives of elements in T, maybe empty;
#          f,   an element in C(x)(T);
#          k1,k2,  K:=k1/k2 is diffrential reduced.
# Output:  [u, h,  p]
#          such that f = u' + K*u + h+p/den(K),
#          where u is in K_n, p is in C(x)(t1,...,t_{n-1})[tn, tn^{-1}] and h is t_n simple  with gcd(den(h), den(K))=1 if T<>[];
#                    u is in C(x), p is in C[x] and h is x-simple  with gcd(den(h), den(K))=1 if T=[];
############################################################################################################################ 

#Name: GKernelAndShellReduction
#Calling sequence:
#                    GKernelAndShellReduction(x, T, dT, f, k1, k2)
# Input:   x,   a variable of the base field C(x);
#          T = [t1, t2, .., tn],   a list of generators of a Liouvillian tower K_n over C(x), maybe empty;
#          dT = [t1', t2', .., tn'],   a list of derivatives of elements in T, maybe empty;
#          f,   an element in C(x)(T);
#          k1,k2,  K:=k1/k2 is diffrential reduced.
# Output:  [u, h,  p]
#          such that f = u' + K*u + h+p/den(K),
#          where u is in K_n, p is in C(x)(t1,...,t_{n-1})<t_n> and h is t_n simple  with gcd(den(h), den(K))=1 if T<>[];
#                    u is in C(x), p is in C[x] and h is x-simple  with gcd(den(h), den(K))=1 if T=[];
############################################################################################################################ 

GKernelAndShellReduction := proc(x, T, dT, f, k1, k2)
   local t, F, num, den, W, A, h1, S, R;

   if T = [] then
      t :=x;
      F:=[f,0];
   else
      t := op(-1, T);
      if degree(dT[-1],t)=1 then
        F := Rational:-NormallyProperAndLaurentPolynomialParts (f, t);
      else 
        F:=[Rational:-ProperAndPolynomialParts(f,t)];
      end if;
   end if;

   
      if F[1] = 0 then
         return ( [0, 0 ,normal(k2*F[2])] );
      end if; 

     #decompose  proper part of f

     F[1] :=normal(F[1]);
     (num, den) := numer(F[1]), denom(F[1]);
     W := RelativeFactor(den, k2, t);
     A  := Rational:-ExtendedEuclidean(t, num, W[2], W[1]);   #A[1]/W[1]+A[2]/W[2]=F[1]
     h1 := normal(A[2]/W[2]); 
     
     #Shell reduction

     S := GShellReduction(x, T, dT, normal(A[1]/W[1]), k1, k2);
     if h1 = 0 then
        return ([S[1], S[2], normal(S[3]+F[2]*k2)]);
     end if;

     #kernel reduction

     R := GKernelReduction(x, T, dT, h1, k1, k2);
     return([normal(S[1]+R[1]), S[2], normal(R[2]+S[3]+F[2]*k2)]);

end proc: 

###########################################################
GKernelAndShellReduction0 := proc(x, T, dT, f, k1, k2)
   local t, F, num, den, W, A, h1, S, R;

   if T = [] then
     t :=x;
     F:=[f,0];
   else
     t := op(-1, T);
     F := Rational:-NormallyProperAndLaurentPolynomialParts (f, t);
   end if;

   
      if F[1] = 0 then
         return ( [0, 0 ,normal(k2*F[2]),0] );
      end if; 
     (num, den) := numer(F[1]), denom(F[1]);
     W := RelativeFactor(den, k2, t);
   if W[1] = [] then
     return ([0, 0,normal(k2*F[2]), F[1]]);
  end if;

  A  := Rational:-ExtendedEuclidean(t, num, W[2], W[1]);   #A[1]/W[1]+A[2]/W[2]=F[1]
  h1 := normal(A[2]/W[2]); 

  S := GShellReduction(x, T, dT, normal(A[1]/W[1]), k1, k2);
  if h1 = 0 then
     return ([S[1], S[2], normal(S[3]+F[2]*k2)]);
  end if;

  R := GKernelReduction(x, T, dT, h1, k1, k2);
  return([normal(S[1]+R[1]), S[2], normal(R[2]+S[3]+F[2]*k2)]);

end proc: 

###########################################################

# GFindSepcialPolyWildCase(a, b, y)
# Input: a, b, two polynomials 
#        y, a variable name.
# Output: p, a polynomial in M_{a/b} and deg(p, y) < deg(a, y),
#         q, a polynomial such that
#            b*Dy(y^t)+a*y^t = b*Dy(q)+a*q + p (t=-la/lb).  
###########################################################

GFindSepcialPolyWildCase := proc(a, b, y)
  local t, la, lb, da, db, dd, p, q, lm, i;

    la := lcoeff(a, y);
    lb := lcoeff(b, y);
    da := degree(a, y);
    db := degree(b, y);
    t  := normal(-la/lb);
    if da=db-1 and type(t, integer) and t>0 then  
     q := y^t;
     p := collect(normal(t*b*y^(t-1)+a*y^t), y);
     dd := degree(p, y)-da;
     while dd >=0 do
         lm := coeff(p, y, da+dd)*y^dd/(lb*dd+la);
          p := normal(normal(p - b*diff(lm, y))-a*lm);
          q := normal(q - lm);
          dd := degree(p, y)-da; 
     end do;
     return([p, q]);
    end if;
    return([0, 0]);

end proc:

##########################################################

# Name: PolynomialReduction2
# Calling sequence:
#                  PolynomialReduction2(p, a, b, sp, y)
# Input: p, a, b,  three polynomials;
#        sp, a pair of polynomials; (used in the wild case);
#        y,        a variable name.
# Output: [p1, p2], where p1, p2 are polynomials satifying
#
#                   p = b* Dy(p1) + a* p1 + p2, and   
#                     
#                 p2 is in the complement of the reduction 
#                 subspace M_k := {b*Dy(f)+a*f | f in E[y]}.
#
###########################################################

PolynomialReduction2 := proc(p, a, b, sp, y)
    local la, lb, da, db, t, p1, p2, i, dd, tm, lm, dsp;
 
    ### Trivial case K=0 or p=0
    if normal(p)=0 then
       return([0, 0]);
    end if;
    if a=0 and b=1 then
       return([int(p, y), 0]);
    end if;

    la := lcoeff(a, y);
    da := degree(a, y);
    lb := lcoeff(b, y);
    db := degree(b, y);
     t := - normal(la/lb);

    p1 := 0;
    p2 := collect(p, y);
    ### print([da, db]);
     
   ### Case 1. da >= db
    
    if da >= db then
       ## print(case_1); 
       dd := degree(p2, y) - da;
       while dd >= 0 do
             lm := lcoeff(p2, y)*y^dd/la; 
             p1 := GMyRatSum(p1, lm);
             p2 := GMyRatSum(GMyRatSum(p2, -b*diff(lm, y)), -a*lm);
             dd := degree(p2, y) - da;
       end do;  
       return([p1, p2]);   
    end if;

   ### Case 2. da < db -1
    
    if da < db-1 then
       ### print(case_2);
       dd := degree(p2, y) - db;
       while dd >= 0 do
             lm := lcoeff(p2, y)*y^(dd+1)/(lb*(dd+1));
             p1 := GMyRatSum(p1, lm);
             p2 := GMyRatSum(GMyRatSum(p2, -b*diff(lm, y)), -a*lm);
             dd := degree(p2, y) - db;
       end do;  
       lm := coeff(collect(p2, y), y, da)/la;
       p1 := GMyRatSum(p1, lm);
       p2 := GMyRatSum(p2, -a*lm); 
       return([p1, p2]);  
    end if;

    ### Case 3. da = db -1 and t is a positive integer
     if type(t, integer) and t > 0 then
        ### print(case_3_wild); 
        tm := coeff(collect(p2, y), y, da+t)*y^(da+t);
        p2 := normal(p2 - tm);
        dd := degree(p2, y) - da;
       while dd >= 0 do
             lm := lcoeff(p2, y)*y^dd/(lb*dd+la);
             p1 := GMyRatSum(p1 , lm);
             p2 := GMyRatSum(GMyRatSum(p2,  -b*diff(lm, y)), -a*lm);
             tm := tm + coeff(collect(p2, y), y, da+t)*y^(da+t);
             p2 := normal(p2 - coeff(collect(p2, y), y, da+t)*y^(da+t));
             dd := degree(p2, y) - da;
       end do;
       dsp := degree(sp[1], y);
       p2 := collect(p2, y);
       p1 := normal(p1 + coeff(p2, y, dsp)*sp[2]/lcoeff(sp[1], y));
       p2 := normal(tm + p2 - coeff(p2, y, dsp)*sp[1]/lcoeff(sp[1], y));
       return([p1, p2]);
     end if;

    ### Case 4. da = db -1 and t is not a positive integer
       ### print(case_4);
       dd := degree(p2, y) - da;
       while dd >= 0 do
             lm := lcoeff(p2, y)*y^dd/(lb*dd+la);
             p1 := GMyRatSum(p1,  lm);
             p2 := GMyRatSum(GMyRatSum(p2, - b*diff(lm, y)), -a*lm);
             dd := degree(p2, y) - da;
       end do;
       return [p1, p2];

end proc:

##################
## GMyRatSum(f, g)
##
##################

GMyRatSum :=proc(f, g)
  local df, dg, nf, ng, gd, cdf, cdg;

  nf, ng := numer(f), numer(g);
  df, dg := denom(f), denom(g);
  
  gd := gcd(df, dg, cdf, cdg);
  return normal(nf*cdg+ng*cdf)/(gd*cdf*cdg);
 
end proc:
###########################################################
# Name: PolynomialReduction
# Calling sequence:
#                  PolynomialReduction(p, a, b, y)
# Input: p, a, b,  three polynomials;
#        y,        a variable name.
# Output: [p1, p2], where p1, p2 are polynomials satifying
#
#                   p = b* Dy(p1) + a* p1 + p2, and   
#                     
#                 p2 is in the complement of the reduction 
#                 subspace M_k := {b*Dy(f)+a*f | f in E[y]}.
#
###########################################################
PolynomialReduction := proc(p, a, b, y)
  local sp;
  sp := GFindSepcialPolyWildCase(a, b, y);
  return PolynomialReduction2(p, a, b, sp, y);

end proc:
##########################################################

#Name: ResidualForm0
#Calling sequence:
#                ResidualForm0(f, k1, k2, x)
#Input: f, a rational function in C(x);
#         k1, k2 are polynomial in C[x] with K:=k1/k2 differential reduced;
#         x, a variable name.
#Output: [u, h, p] such that
#           f=u'+K*u+h+p/k2;
#           u is in C(x),  h is x-simple  with gcd(den(h), den(K))=1, and
#           p is in the complement of the reduction subspace 
#              M_k := {k2*v'+k1*v | v in C[x]}.

###############################################################

ResidualForm0 :=proc(f, k1, k2, x)
local A, B;
  A := GKernelAndShellReduction(x, [], [], f, k1, k2);
  if A[3] = 0 then
     return A;
  end if;
  B := PolynomialReduction(A[3], k1, k2, x);
  return ([normal(A[1]+B[1]), A[2], B[2]]);
 end proc:



#############################################################
#Name: LaurentPolynomialReduction1
#Calling sequence:
#                LaurentPolynomialReduction1(x, t, dt, p, k1,k2)
#Input: x, a variable of the base field C(x);
#          t, an exponential generator over C(x);
#          dt, the derivative of t;
#          p, a polynomial in C(x)[t];
#          k1, k2,  polynimials of C(x)[t] satisfies K:=k1/k2 differential reduced
#Output: [u, p1] such that
#            p=k2*u'+k1*u+p1, where u is in C(x)(t), p1 is in C_k,
#            C_k := span_C(x){t^i | 0<=i<m}+b_m*span_C{ r_i*t^i | i>=m}
#            where r_i is a ((i-m)t'/t+a_m/b_m)-residual form in C(x) for i >=m if b_m<>0 
#           write  k2 =b_m*t^m+.......+b_0 and k1= a_m*t^m+.......+a_0
#####################################################################

LaurentPolynomialReduction1:=proc (x, t,dt, p, k1,k2)
 local m, a, b, n1, n2, pp,hc, u,  A, g1, g2, r, q;
    m := max(degree(k1, t), degree(k2, t));
     a := coeff(k1, t, m);
     b := coeff(k2, t, m);
    # e  := coeff(k1, t, 0);
     #f := coeff(k2, t, 0);
    #(num, den) := numer(p), denom(p);
    # pr := rem(num, den, t, 'pp');
     u := 0;
     q := 0;
      pp:=p;
      n1 := degree(pp, t);

       while n1>=m do
          hc := coeff(pp, t, n1);
          
          #reduce the leading term of p

            if a<>0 then
              u :=normal(u+(hc/a)*t^(n1-m));
              pp:=normal(pp-k2*DField:-Derivative(x, [t],[dt], u)-k1*u);
            else
               (g1, g2) :=numer((n1-m)*(dt/t)+a/b), denom((n1-m)*(dt/t)+a/b);
               A := ResidualForm0(hc/b, g1, g2, x);  #hc/b=A[1]'+((n1-m)*(dt/t)+a/b)*A[1]+A[2]+A[3]/g2;
               r := normal(A[2]+A[3]/g2);
               u := normal(u+A[1]*t^(n1-m));
               q := normal(q+r*b*t^n1);
               pp := normal(pp-k2*DField:-Derivative(x, [t],[dt], A[1]*t^(n1-m))-k1*(A[1]*t^(n1-m))-r*b*t^n1);
              end if;
              n1 := degree(pp, t);
         end do; 
             return([u, normal(pp+q)]);
   end proc:

#############################################################
#Name: LaurentPolynomialReduction2
#Calling sequence:
#                LaurentPolynomialReduction2(x, t, dt, p, k1,k2)
#Input: x, a variable of the base field C(x);
#          t, an exponential generator over C(x);
#          dt, the derivative of t;
#          p, a laurent polynomial in C(x)[t^{-1}];
#          k1, k2,  polynimials of C(x)[t] satisfies K:=k1/k2 differential reduced
#Output: [u, p1] such that
#            p=k2*u'+k1*u+p1, where u is in C(x)(t), p1 is in C_k,
#            C_k := span_C(x){t^i | 0<=i<m}+b_0*span_C{ r_i*t^i | i<0};
#            where  r_i is a (it'/t +a_0/b_0)-residual form  for i<0 if b_0<>0.
#           write  k2 =b_m*t^m+.......+b_0 and k1= a_m*t^m+.......+a_0
#####################################################################
LaurentPolynomialReduction2:=proc(x, t, dt, p, k1,k2)
  local m, e, f, u, q, pr, n2, tc, g1, g2, A, r;
    m := max(degree(k1, t),degree(k2, t));
     e  := coeff(k1, t, 0);
     f  := coeff(k2, t, 0);
     u := 0;
     q := 0;
     pr := p;
      n2 := ldegree(pr, t);
      while n2<0 do
            tc := coeff(pr, t, n2);
            
            #reduce the tail term of pr
            
            if e<>0 then
               u :=normal(u+(tc/e)*t^(n2));
              pr :=normal(pr-k2*DField:-Derivative(x, [t],[dt], u)-k1*u);
            else
               (g1, g2) :=numer((n2)*(dt/t)+e/f), denom((n2)*(dt/t)+e/f);
               A := ResidualForm0(tc/f, g1, g2, x);  #tc/f=A[1]'+((n2)*(dt/t)+e/f)*A[1]+A[2]+A[3]/g2;
               r := normal(A[2]+A[3]/g2);
               u := normal(u+A[1]*t^(n2));
               q := normal(q+r*f*t^(n2));
               pr := normal(pr-k2*DField:-Derivative(x, [t],[dt], A[1]*t^(n2))-k1*A[1]*t^(n2)-r*f*t^(n2));
              end if;
              n2 := ldegree(pr, t);
        end do;
       return([u, normal(q+pr)]);
end proc:

#############################################################
#Name: LaurentPolynomialReduction
#Calling sequence:
#                LaurentPolynomialReduction(x, t, dt, p, k1,k2)
#Input: x, a variable of the base field C(x);
#          t, an exponential generator over C(x);
#          dt, the derivative of t;
#          p, a laurent polynomial in C(x)[t, t^{-1}];
#          k1, k2,  polynimials of C(x)[t] satisfies K:=k1/k2 differential reduced
#Output: [u, p1] such that
#            p=k2*u'+k1*u+p1, where u is in C(x)(t), p1 is in C_k,
#            C_k := span_C(x){t^i | 0<=i<m}+b_m*span_C{ r_i*t^i | i>=m}+b_0*span_C{ r_i*t^i | i<0};
#            where r_i is a ((i-m)t'/t+a_m/b_m)-residual form in C(x) for i >=m if b_m<>0 or r_i is a 
#            (it'/t +a_0/b_0)-residual form  for i<0 if b_0<>0.
#           note thatC_k is the complement space of the reduction 
#                 subspace  I_k:={k2*v'+k1*v|v in C(x)[t, t^{-1}]} 
#           write  k2 =b_m*t^m+.......+b_0 and k1= a_m*t^m+.......+a_0
#####################################################################
LaurentPolynomialReduction := proc(x, t, dt, p, k1,k2)
local num, den, pp, pr, A, B;
   (num, den) := numer(p),denom(p);
    pr := rem(num, den, t, 'pp');
    pr := normal(pr/den);
    A := LaurentPolynomialReduction1(x, t, dt, pp, k1,k2);
    B :=[];
    if pr = 0 then
       return(A);
    end if;
    B :=LaurentPolynomialReduction2(x, t, dt, pr, k1,k2);
    return([normal(A[1]+B[1]), normal(A[2]+B[2])]);
 end proc:




       
  
end module:
         



