## Copyright (c) January 2013, NCSU. All rights reserved.
## Author: Shaoshi Chen <schen@amss.ac.cn> 
## Load the modules: 
##  Decomposition.mm, Division.mm, HyperexpReduction.mm, List.mm, Hyperexp.mm
##
## Modified on April 12, 2013 (+ a code to add rational functions)
################################################################################

HermiteCT := module()
    description "Creative Telescoping for Hyperexponential Functions via Hermite-Like Reduction":

    export
                    
          HermiteReduction,
          HermiteTelescoping,
          KernelReduction, 
          PolynomialReduction,  
          ShellReduction;                          
    local
          PolynomialReduction2,
          SimpleReduction,
          MyRatSum,
          OrderBound,
          FindSepcialPolyWildCase,
          ExtendedEuclidean;       
    option package:

#######################  Data  Description ####################################
# List representation:
# A hyperexponential function H(x1, ..., xn) is representated 
# by the following list
#
#              [[c1, x1], ..., [cn, xn]],  ----------------------  (*)
#
# where ck's are rational functions and xk's are variable names s.t.
#
#           ck = diff(H, xk)/H, for k = 1, ..., n, 
#
# and ck's satisfy compatible conditions
#
#           diff(ci, xj) = diff(cj, xi) for all i<>j in {1, ..., n}.
# Remark: if H = 0, we set [] its list representation.
#
###############################################################################


#######################  Task  Description ####################################
# Telescoping problem for hyperexponential functions:
# Given H(x, y) a hyperexponential function with list form [[f, x], [g, y]], find
# a nonzero linear differential operator L(x, Dx) in F(x)<Dx> such that
# 
#                   L(x, Dx)(H) = Dy(T),
#
# where T is another hyperexponential function.
#
###############################################################################


##################################
#                                #
#     Export Functions           #
#                                #
##################################

#############################################################################
# Name: HermiteReduction
# Calling sequence:
#                             HermiteReduction(H, y)
# Input: H,       a hyperexponential function with Dy(H)/H = Dy(S)/S+K;
#        y,       a variable name.
#
# Output:  u, r,  two rational functions;
#          T,     a hyperexponential function with certificate K such that
#
#             1.  r is the residual form of S w.r.t. K
#
#             2.  H = Dy(u*T) + r*T.
#
##############################################################################

HermiteReduction := proc(H, y)

 local g, rcf, S, K, s1, s2, sp, gex, v, k1, k2, T, u, r, a, b, 
       i, as, ps, sred, q2, q1, pred; 

#### Find the certificate   ######
  
     g := Hyperexp:-FromFunctionToList(H, [y])[1][1]:
 

#### Trivial case: H is free of y ###
  if g = 0 then
     return([0, 1, 0, 1, 0, H]);
  end if;

#### Compute the canonical DRNF of g ###
  rcf := HyperexpReduction:-RationalCanonicalForm(g, y);
  S := rcf[1];
  K := rcf[2];
  s1, s2 := numer(S), denom(S);
  k1, k2 := numer(K), denom(K);
  T := normal(simplify(H/S));
  sp := FindSepcialPolyWildCase(k1, k2, y);
  

#### Shell reduction of S w.r.t. K ###
  sred := ShellReduction(S, K, y);
  r, a, b, i := sred[1], sred[2], sred[3], sred[4];

  
#### Polynomial reduction for a/(b*k2^i) ###

  ps := quo(a, normal(b*k2^i), y, 'as');
  if i = 0 then
     q1 := as;
     q2 := normal(ps*k2);
  else
     gex := ExtendedEuclidean(as, b, k2, y);
     q1 := gex[2];
     ## q2 := normal(gex[1]+ps*k2);
     q2 := MyRatSum(gex[1], ps*k2);
  end if;

  pred := PolynomialReduction2(q2, k1, k2, sp, y);
     v := normal(pred[2]);
   ##r := normal(r + pred[1]);
   r := MyRatSum(r, pred[1]);
  
  return([r, MyRatSum(q1/b, v/k2), T])
end proc;

#############################################################################
# Name: HermiteTelescoping
# Calling sequence:
#                           HermiteTelescoping(H, x, y, Dx, 'T')
# Input: H,      a bivariate hyperexponential function;
#        x, y,   two variable names;
#        Dx,     a operator name.
# Output: L,     a linear differential operator in Dx over F(x) and 
#                the fifth argument 'T' is a bivariate hyperexponential 
#                function such that
#
#                                  L(x, Dx)(H) = Dy(T).
#
#          T is of the form r*exp(int(udx + vdy)) 
#              with r a rational function and v differential-reduced w.r.t. y.
##############################################################################

HermiteTelescoping := proc(H, x, y, Dx, T)

 local f, g, vars, S, K, sp, s1, estord, s2, k1, k2, rcf, r, a, as, b, i, j, HL, Ts, Hs, opt,
       sred, kred, pred, sr, qL, vL, rL, tt, rcert, q1, q2, ps, gex, cL, testsol, L;

#### Initialize   ######
  rcert := 0;
  qL := table([]);
  vL := table([]);
  cL := table([]);
  rL := table([]);
  opt := 0;
 
#### Certificates ##### 
  HL := Hyperexp:-FromFunctionToList(H, [x, y]):
  f  := Hyperexp :-FindCertificate(HL, x)[1];
  g  := Hyperexp :-FindCertificate(HL, y)[1];
  vars := Hyperexp :-Variables(HL);

#### Trivial case: H is free of y ###
  if g = 0 then
     if nargs = 5 then
        T := [y, H];
        return(1);
     else
        return(1);
     end if;
  end if;

#### Compute the canonical DRNF of g ###
  rcf := HyperexpReduction:-RationalCanonicalForm(g, y);
  S := rcf[1];
  K := rcf[2];
  s1, s2 := numer(S), denom(S);
  k1, k2 := numer(K), denom(K);
  sp := FindSepcialPolyWildCase(k1, k2, y);
  Ts := normal(simplify(H/S));
  Hs := Hyperexp:-HyperexpQuotient(HL, Hyperexp:-FromFunctionToList(S, vars));

#### Shell reduction of S w.r.t. K ###
  sred := ShellReduction(S, K, y);
  r, a, b, i := sred[1], sred[2], sred[3], sred[4];

  if nargs = 5 then
     opt := 1;
     rcert := MyRatSum(rcert, r);
  end if;

  ### Estimate the order of minimal telescopers

     estord := OrderBound(b, k1, k2, y);

     printf("The estimated order of minimal telescopers is %d \n ", estord);

  
#### Polynomial reduction for a/(b*k2^i) ###

  ps := quo(a, b*k2^i, y, 'as');
  if i = 0 then
     qL[0] := as;
     q2 := normal(ps*k2);
  else
     gex := ExtendedEuclidean(as, b, k2, y);
     qL[0] := gex[2];
     q2 := MyRatSum(gex[1], ps*k2);
  end if;
  pred := PolynomialReduction2(q2, k1, k2, sp, y);
  vL[0] := normal(pred[2]);
  if nargs = 5 then
     rcert := MyRatSum(rcert, pred[1]);
     rL[0] := rcert;
  end if;
  

#### First test: whether q0= 0 and v0 = 0  ####
    if qL[0] = 0 and vL[0] = 0 then
       if nargs = 5 then
          T := [rcert, Ts];
          return(1);
       else
          return(1);
       end if;
    end if;
 
##### Find a minimal telescoper for H ####

      testsol := true;
      j := 0;
      ## print([qL[0], vL[0]]);

   while testsol do

         printf("Check order %d \n", j+1);

         ### Simple reduction ####
         sr := SimpleReduction(qL[j], vL[j], b, k1, k2, Hs, sp, x, y, opt);
         ### print(sr);

         qL[j+1] := sr[1];
         vL[j+1] := sr[2];
         if nargs = 5 then
            ## rL[j+1] := (diff(rL[j], x) + rL[j]*Hyperexp :-FindCertificate(Hs, x)[1] + sr[3]);
            ## rL[j+1] := normal(diff(rL[j], x) + rL[j]*Hyperexp :-FindCertificate(Hs, x)[1] + sr[3]);
               rL[j+1] := MyRatSum(MyRatSum(diff(rL[j], x),  normal(rL[j]*Hyperexp :-FindCertificate(Hs, x)[1])), sr[3]);
         end if;

         ### Find linear dependency among {pL, vL} ###

         ### tt := time():

         cL := LinearRelation:-FindLinearRelation([seq(qL[i], i=0..j+1)], [seq(vL[i], i=0..j+1)], y);
         ## tt := time()-tt;
         ###print(tt);

         ### print(cL);

         if nops(cL) = 0 then
            testsol := true;
            j := j+1;
         else
            L := add(cL[i]*Dx^(i-1), i=1..nops(cL));

            if nargs = 5 then
               tt := time():
               rcert := 0;
               for i from 1 to nops(cL) do
                   ## rcert := rcert + cL[i]*rL[i-1];
                   ##  rcert := normal(rcert + cL[i]*rL[i-1]);
                     rcert := MyRatSum(rcert, cL[i]*rL[i-1]);
               end do;
               tt := time()-tt; 
               print(time_for_certificates);
               print(tt);
               T := [rcert, Ts];
            end if;

            return(L);
         end if; 
   end do;

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 := FindSepcialPolyWildCase(a, b, y);
  PolynomialReduction2(p, a, b, sp, y)

end proc;

####################################################################
# Name: KernelReduction
# Calling sequence:
#                   KernelReduction(p, k1, k2, m, y)
# Input: p, k1, k2, three polynomials with K := k1/k2 differential-reduced 
#                   and gcd(p, k2)=1;
#        y,         a variable name.
# Output: [p1, p2], two polynomials such that
#                   p/k2^m = Dy(p1/k2^(m-1)) + p1/k2^(m-1)*K + p2/k2.
#
#####################################################################

KernelReduction := proc(p, k1, k2, m, y)
  local res, p1s, p2s, p1, p2, rds;
  if m = 1 then
     return([0, p]);
  end if;
  
  res := ExtendedEuclidean(p, MyRatSum(k1, -(m-1)*diff(k2, y)), k2, y);
  p1s := res[1];
  ## p2s := normal(res[2]-diff(res[1], y));
  p2s := MyRatSum(res[2], -diff(res[1], y));
  rds := KernelReduction(p2s, k1, k2, m-1, y);
  ###p1 := normal(rds[1]*k2+p1s);
  p1 := MyRatSum(rds[1]*k2, p1s);
  p2 := rds[2];
  
  [p1, p2];
end proc;


####################################################################
# Name: ShellReduction
# Calling sequence:
#                   ShellReduction(S, K, y)
# Input: S, K, two rational functions with K being differential-reduced
#                        and gcd(den(S), den(K)) = 1;
#        y,     a variable name.
# Output: f,   a rational function;
#         a, b, two polynomials;
#         i,    an integer in {0, 1}, such that
# 
#                    S= Dy(f) + f*K + a/(b*den(K)^i).
#####################################################################

ShellReduction := proc(S, K, y)
   local  s1, s2, k1, k2, f, f1, f2, a, b, r, qv, i, Rd, h, p, L1, L2;
   
   s1 := numer(S);
   s2 := denom(S);
   k1 := numer(K);
   k2 := denom(K);
   

   Rd := HyperexpReduction:-Original(S, K, y);
   (h, p, L1, L2) := Rd[1], Rd[2], Rd[3], Rd[4];
   f := List:-FromListToNormal(h);
   f1 := 0;
   f2 := 0;
   for i from 1 to nops(L1) do 
         f1 := MyRatSum(f1,  L1[i][2][1]*L1[i][2][2]/L1[i][1]);
    end do;
    for i from 1 to nops(L2) do 
          f2 := MyRatSum(f2,  L2[i][2][1]*L2[i][2][2]/(L2[i][1]));
     end do;
     r := MyRatSum(MyRatSum(p*k1, f1*k2), f2);
     (a, b) := numer(r), denom(r);
     if rem(a, k2, y, 'qv') = 0 then
         i := 0;
         a := qv;
      else
         i := 1;
      end if;
      [f, a, b, i];
   
end proc;



##################################
#                                #
#     Local Functions            #
#                                #
##################################





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


#####
##
#####
##################
## MyRatSum(f, g)
##
##################

MyRatSum :=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);
  normal(nf*cdg+ng*cdf)/(gd*cdf*cdg);
 
end proc;

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

OrderBound:=proc(b, k1, k2, y)
   if k1 = 0 then
      return(degree(b, y));
   end if;
   if degree(k1, y) < degree(k2, y)-1 then
      return(degree(b, y) + degree(k2, y)-1);
   else
      return(degree(b, y) + degree(k1, y))
   end if;
end proc;

###################################################################
## Name: SimpleReduction
## Calling sequence:
##                       SimpleReduction(q, v, b, k1, k2, Hs, sp, x, y, opt)
## Input: H = (q/b + v/k2)*Hs where Hs = [[a/k2, x], [k1/k2, y]]
##        opt = 0 no need compute the certificate, otherwise opt = 1.
## Output: [qs, vs, r],  such that
##         
##                  Dx(H) = Dy(r*Hs) +(qs/b + vs/k2)*Hs
##         
##
#####################################################################
SimpleReduction := proc(q, v, b, k1, k2, Hs, sp, x, y, opt)

    local f, a, nf, sr, kr, pr, ps, q1, q2, q3, q4, 
          i, u2, w1, w2, vs, r, as, qs, gex;
    
    ### Trivial case
      if q=0 and v=0 then
         return([0, 0, 0]);
      end if;
     
    f := Hyperexp :-FindCertificate(Hs, x)[1];
    nf := normal(f*k2);
    a := nf;
    ### Shell reduction
    sr  := ShellReduction(-q*diff(b, x)/b^2, k1/k2, y);
     r, i  := sr[1], sr[4];
     u2 := normal(sr[2]*b/sr[3]);
    
    ### Kernel reduction
    kr := KernelReduction(normal(nf-diff(k2, x))*v, k1, k2, 2, y);
    w1, w2 := kr[1], kr[2];

    ######################################
    qs := quo(normal(nf*q), b*k2, y, 'as');
    gex := ExtendedEuclidean(as, b, k2, y);
    q1 := gex[2]; 
    q2 := MyRatSum(qs*k2, gex[1]);

    
    qs := quo(u2, b*k2^i, y, 'as');   
    if i = 0 then
       q3 := as;
       q4 := normal(qs*k2);
    else
       gex := ExtendedEuclidean(as, b, k2, y);
       q3 := gex[2];
       q4 := MyRatSum(gex[1], qs*k2);
    end if;   

    ### Polynomial reduction
    ps := MyRatSum(MyRatSum(MyRatSum(diff(v, x), w2), q2),  q4);
    pr := PolynomialReduction2(ps, k1, k2, sp, y);
    if opt = 1 then
       ##r := r+w1/k2+pr[1];
       r := MyRatSum(MyRatSum(r, w1/k2), pr[1]);
    else
       r := r+w1/k2+pr[1];
    end if; 

    qs := MyRatSum(MyRatSum(diff(q, x) , q1),  q3);
    vs := pr[2];
    return([qs, vs, r])  

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 := MyRatSum(p1, lm);
             p2 := MyRatSum(MyRatSum(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 := MyRatSum(p1, lm);
             p2 := MyRatSum(MyRatSum(p2, -b*diff(lm, y)), -a*lm);
             dd := degree(p2, y) - db;
       end do;  
       lm := coeff(collect(p2, y), y, da)/la;
       p1 := MyRatSum(p1, lm);
       p2 := MyRatSum(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 := MyRatSum(p1 , lm);
             p2 := MyRatSum(MyRatSum(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 := MyRatSum(p1,  lm);
             p2 := MyRatSum(MyRatSum(p2, - b*diff(lm, y)), -a*lm);
             dd := degree(p2, y) - da;
       end do;
       [p1, p2]

end proc;

###########################################################
# FindSepcialPolyWildCase(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).  
###########################################################
FindSepcialPolyWildCase := 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: ExtendedEuclidean
# Calling sequence:
#                  ExtendedEuclidean(p, a, b, x)
# Input: p, a, b,  three polynomials with gcd(a, b)=1;
#        x,        a variable name.
# Output:  u, v,   two polynomial such that
#                   p = u*a + v*b. 
#
#####################################################################
ExtendedEuclidean := proc(p, a, b, x)
  local h, u, v, pp1, pp2, t, c1, c2, q;
  
  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 := divide(pp1, pp2, 'q');
  if t = false then 
      error "the division is not exact";
  end if;
  v := normal(c1/c2)*q;
  [u, v]
end proc;




end module:




############################
##
## Some auxiliary modules
##
##
############################

## Copyright (c) 2010, INRIA & KLMM. All rights reserved.
## Authors: Shaoshi Chen <schen@amss.ac.cn> and Ziming Li <zmli@mmrc.iss.ac.cn>.



Decomposition := module()
   description "this module is to compute squarefree partial fraction decomposition";
   export  SquareFree,
           SquareFreePartialFraction,
           VerifyPartialFraction;
   local 
         AdicExpansion,
         LaurentPart, 
         Position; 
   option package;

 ################################################################
 # Name: AdicExpansion                                          
 # Calling sequence:  
 #            AdicExpansion(F, x, t)              
 # Input: F, a nonzero fraction [a, d, n] representing a simple fraction
 #        a/d^n;                                       
 #        x, a name of variable.                                
 # Output: the list of the d-adic expansion of a/d^n,                           
 #         [[c1, a_1] ..., [cn, a_n]];                          
 # Optional: the third argument if presented is assigned the d-adic
 #           expansion in the usual form                       
 ################################################################

 AdicExpansion := proc(F, x, t)
    local a, c, d, n, i, r, q, T, l, g;
    (a, d, n) := F[1], F[2], F[3]; 
    T := table([]);
    if nargs = 3 then
       g := 0;
    end if;
    for i from n by -1 to 1 do 
       r := rem(a, d, x, 'a');
       r := primpart(r, x, 'c');
       if r <> 0 then
          T[i] := [c, r];
       else
          T[i] := [1, 0];
       end if;  
       if nargs = 3 then
          g := c*r/d^i + g;
       end if;
    end do;
    if nargs = 3 then
       t := g;
    end if;
    l := n;
    for i from n by -1 to 1 do
        if T[i][2] = 0 then
           l := l - 1;
        else
           break;
        end if;
    end do; 
    [seq(T[i], i=1..l)]; 
 end proc:

 ###########################################################################
 # Name: Position
 # Calling sequence:
 #       Position(L, n)
 # Input: L, a list of integers arranged increasingly;
 #        n, an integer
 # Output: l, a positive integer such that L[l-1] < n <= L[l];
 ###########################################################################

 Position := proc(L, n)
    local l, m, i; 
    l := 1; 
    m := nops(L);
    for i from 1 to m do
        if L[i] < n then
           l := l + 1;
        else
           break;
        end if; 
    end do;
    l;
 end proc;
      
 ###########################################################################
 # Name: SquareFree
 # Calling sequence:
 #         SquareFree(P, x)
 # Input: P, a nonzero polynomial
 #        x, a name
 # Output: a triple, c, L, M, where
 #         c is the content of P with respect to x;
 #         L is the list of squarefree factors
 #         M is the list of multiplicities corresponding to the factors in L
 #         (M is arranged increasingly)
 ###########################################################################

  SquareFree := proc(P, x)
     local T, L, M, Lp, Mp, Ls, Ms, c, n, i, f, m, k;  

     if degree(P, x) = 0 then
        return [P, []];
     end if;

     T := sqrfree(P, x);

     # rearrange the output

     c := T[1];
     n := nops(T[2]);
     (L, M) := [T[2][1][1]], [T[2][1][2]];
     for i from 2 to n do
         (f, m) := T[2][i][1], T[2][i][2]; 
         k := Position(M, m);
         Ls := [op(1..k-1, L)];
         Lp := [op(k..nops(L), L)];
         Ms := [op(1..k-1, M)];
         Mp := [op(k..nops(M), M)];
         while Mp <> [] and m = Mp[1] do
               f := f*Lp[1];
               Lp := [op(2..nops(Lp), Lp)];
               Mp := [op(2..nops(Mp), Mp)];
         end do;
         L := [op(Ls), f, op(Lp)];
         M := [op(Ms), m, op(Mp)];
     end do;
     [c, [seq([L[i], M[i]], i=1..nops(M))]];
  end proc;
    
###########################################################
# Name: LaurentPart
# Calling sequence:
#       LaurentPart(a, d, s, x, t)
# Input: a, a polynomial in x;
#        d, a list [u, v, n], representing a polynomial
#           u*v^n with gcd(u, v) = 1 and gcd(a, u*v)=1.
#        s, a polynomial such that s*u = 1 mod v
#        x, a name
# Output: c, L,  where c is the content, and L = [a1, .., an] such that
#         the Laurent serise of a/(u*v^n) at v truncated
#         at the term of order -1 is equal to
#
#         c*(a_1/v_1 + ... + a_n/v^n)
#
# Optional: the fifth argument is assigned a new numerator 
###########################################################

LaurentPart := proc(a, d, s, x, t)
   local nu, u, v, n, T, r, i, cp, w, ds, ns, l, rs, cs;
   nu := a;
   (u, v, n) := d[1], d[2], d[3];
   T := table([]); 
   cp := 1; ds := 1; 
   for i from n by -1 to 1 do
       cp := normal(cp*content(nu, x, 'nu')*(1/ds)); 
       r :=  rem(rem(nu, v, x)*s, v, x);
       if r = 0 then
          T[i] := [1, 0];
       else
          rs := primpart(r, x, 'cs');
          T[i] := [normal(cs*cp), rs];   #modiefied  on July 6
       end if; 
       (ds, ns) := denom(r), numer(r);
       nu := ds*nu - ns*u;
       nu := Division:-ExactDivision(nu, v, 1, x);
   end do;
   if nargs = 5 then
      t := normal(cp/ds)*nu;
   end if;
   l := n;
   for i from n by -1 to 1 do
       if T[i][2] = 0 then
           l := l - 1;
       else
           break;
       end if;
    end do; 
    [seq(T[i], i=1..l)]; 
end proc; 
 
##############################################################################
#Name: SquareFreePartialFraction
#Calling sequence
#         SquareFreePartialFraction(f, x, t)
#Input: f, a nonzero rational function of x;
#       x, a name
#Output: a pair c, G, where c is the content of f with respect to x,
#        and G is the list of partial fractions of f/c, p is the polynomial.
#Optional: the third argument, if presented, is assigned the usual expression
#          of squarefree partial fraction decompositions 
###############################################################################

SquareFreePartialFraction := proc(f, x, t)
   local g, a, b, d, c, p, i, S, L, M, n, Gt, w, r, 
         v, m, u, H, s, G, T, T1, T2, j, dp, cp, ap, bp, Sp, cs;

   # initialize

   g := normal(f);
   (a, d) := numer(g), denom(g);
   c := content(a, x, 'ap')/content(d, x, 'dp');

   # trivial case

   if degree(d, x) = 0 then
      if nargs = 3 then
         t := c*ap;
      end if;
      return(c, [x, ap, []]) # modified by Shaoshi 0703
   end if;

   # compute polynomial part

   p := quo(ap, dp, x, 'bp');
   ap := bp;

   # compute the squarefree factorization of the denominator
   
   Sp := SquareFree(dp, x);
   cs := Sp[1];
   dp := cs*dp;
   
   # compute the Laurent parts one by one

   S := Sp[2];
   n := nops(S);

   S := [seq(S[n-i], i=0..n-1)]; 

   Gt := table([]);

   
   for i from 1 to n-1 do

       # prepare for Laurent expansion

       v := S[i][1];
       m := S[i][2];
       u := Division:-ExactDivision(dp, v, m, x); 

       T := table([]);
       for j from i+1 to n do
           gcdex(S[j][1], v, x, 'w');
           T[j] := w$S[j][2]; 
       end do;
       H := [seq(T[j], j=i+1..n)];
       s := Division:-MergeRemainder(H, v, x);

       # compute Laurent expansion and update the numerator


       Gt[i] := LaurentPart(ap, [u, v, m], s, x, 'bp');
       ap := bp; 
     
       # update the denominator
              
       dp := u;

  end do;
  Gt[n] := AdicExpansion([ap, S[n][1], S[n][2]], x);

  #######################################
  # Prepare for return
  #######################################

  G := [];
  for i from 1 to n do
      if cs <> 1 then
         T1 := Gt[i];
         T2 := [];
         for j from 1 to nops(Gt[i]) do
             T2 := [ op(T2), [T1[j][1]*cs, T1[j][2]]];
         end do;
         Gt[i] := T2;
      end if;
      G := [op(G), [S[i][1], op(Gt[i])]];
  end do;
  G := [x, p, G];
  if nargs = 3 then
     t := c*List:-FromListToPartialFraction(G);
  end if; 
  c, G;
end proc; 


#######################################
# Verify(c, G, r)
# check the correctness of the output of
#  SquareFreePartialFraction
#######################################

VerifyPartialFraction := proc(c, G, r)
  local var, p, F, l, i, f, de, g, d, nu, h, j;
  (var, F) := G[1], G[3];
  l := nops(F);
  for i from 1 to l do 
      f := F[i];
      de := f[1];
      g := gcd(de, diff(de, var));
      if g <> 1 then 
         return false;
      end if; 
      d := degree(de, var);
      for j from 2 to nops(f) do
          nu := f[j][2];
          if degree(nu, var) >= d then
             return false;
          end if;
      end do;
  end do; 
  h := List:-FromListToPartialFraction(G);
  Testzero(c*h-r);
end proc; 
end module: 
                    
       
## Copyright (c) 2010, INRIA & KLMM. All rights reserved.
## Authors: Shaoshi Chen <schen@amss.ac.cn> and Ziming Li <zmli@mmrc.iss.ac.cn>.

Division := module()
  description "polynomial extended gcd computation"; 
  export 
         ExactDivision, 
         HalfEuclideanList,
         MergeInverse,
         MergeRemainder;
  option package;



 ####################################################################
 # Name: HalfEuclideanList                                          #
 # Calling sequence:                                                #
 #                    HalfEuclideanList(L, g, x)                    #
 # Input: L, a list of nonzero polynomials in x;                    #
 #        g, a nonzero polynomial in x such that f and g are        #
 #           co-prime for all f in L;                               #
 #        x, a name.                                                #
 # Output: C, a list of polynomials such that C[i]*L[i] =1 mod g;   #
 ####################################################################
 

 HalfEuclideanList := proc(L, g ,x)
     local C, i, h, s;
     C := [];
     for i from 1 to nops(L) do
       h := gcdex(L[i], g, x, 's');
       if degree(h, x) > 0 then
          error "a nontrivial gcd is found";
       end if;
       C := [op(C), s];
     end do; 
 end proc;


 ###############################################################
 # Name: ExactDivision                                         #
 # Calling sequence:                                           #
 #                    ExactDivision(P1, P2, m, x)              #
 # Input: P1, P2, two polynomials in x;                        #
 #        m, a positive integer                                #
 #        x, a name                                            #
 # Output: q,  q is a polynomial such that q is the            #  
 #             quotient of P1 by P2^m.                         #         
 ###############################################################


 ExactDivision :=proc(P1, P2, m, x) 
   local c1, c2, q, i, pp2, t;
  
       q := primpart(P1, x, 'c1');
       pp2 := primpart(P2, x, 'c2');
       for i from 1 to m do
           t := divide(q, pp2, 'q');
           if t = false then 
              error "the division is not exact";
           end if;
       end do;
       q*c1/c2^m;
 end proc;


#Division :=proc(P1, P2, m, x) 
#   local c1, c2, q, i, pp2, t;
# 
#       tt := time(); 
#       q, c1 := numer(P1), denom(P1); 
#       pp2, c2 := numer(P2), denom(P2);
#       print(time_taken_for_primpart);
#       tt := time()-tt;
#       print(tt); 
#       tt := time();
#       for i from 1 to m do
#           t := divide(q, pp2, 'q');
#           if t = false then 
#              error "the division is not exact";
#           end if;
#       end do;
#       tt := time()-tt;
#       print(time_taken_for_divide); print(tt);
#       q*c2^m/c1;
# end proc;
  

 ################## Merge Remainder ##########################
 # Calling sequence:                                         #
 #                   MergeRemainder(L, b, x)                 #
 # Input: L, a list of polynomials;                          #
 #        b, a nonzero polynomial;                           #
 #        x, a name of variable;                             #
 # Output:  the remainder of the product of all elements     #
 #          in L by b w.r.t x.                               #
 #############################################################

MergeRemainder := proc(L, b, x)
    local n, np, r, i, L1, L2;

    n := nops(L);

    if degree(b, x) = 0 then
       return 0;
    end if;

    if n = 1 then
       return rem(L[1], b, x);
    end if;
    np := floor(n/2);
    L1 := [seq(L[i], i=1..np)];
    L2 := [seq(L[i], i=np+1..n)];
    r  := rem(MergeRemainder(L1, b, x)*MergeRemainder(L2, b, x), b, x);
 end proc:


 ####################################################################
 # Name: MergeInverse                                               #
 # Calling sequence:                                                #
 #                   MergeInverse(L, g, C, x,'t')                   #
 #                                                                  #
 # Input: L, a list of polynomials co-prime with g;                 #
 #        g, a nonzero polynomial;                                  #
 #        x, a name;                                                #
 #        C, a list of polynomials such that                        #
 #           C[i]*L[i] =1 mod g;                                    #
 #        t, (optional) unevaluated name.                           #
 # Output: a polynomial s such that                                 #
 #                   s*(Prod_i L[i]) = 1 mod g.                     #
 # Option: if the fifth argument t is presented, then it is         #
 #         assigned the polynomial such that                        #
 #                 s*(Prod_i L[i]) + t*g = 1.                       #
 ####################################################################



 MergeInverse := proc(L, g, C, x, t)

   local s, h, f;
   s := MergeRemainder(C, g, x);
   if nargs = 5 then
      h := mul(L[i], i=1..nops(L));
      f := 1 - h*s;
      t := ExactDivision(f, g, x);
   end if;
   s; 
 end proc:

end module:


Hyperexp := module()
    description "Basic operations for hyperexponential functions";
    
    export 
          AreSimilar,
          AreCompatible,
          AreTheSame,
          Certificates, 
          FromFunctionToList,
          FindCertificate,
          HyperexpDiff,
          HyperexpProduct,
          HyperexpQuotient,
          HyperexpSum,
          IsConstant,
          IsPolynomial,
          IsRational,
          Variables;

    local 
          IsPolynomial_1,
          IsRational_1;
                     
    option package;



####################### Data  Description #####################################
# List representation:
# A hyperexponential function H(x1, ..., xn) is representated 
# by the following list
#
#              [[c1, x1], ..., [cn, xn]],  ----------------------  (*)
#
# where ck's are rational functions and xk's are variable names s.t.
#
#           ck = diff(H, xk)/H, for k = 1, ..., n, 
#
# and ck's satisfy compatible conditions
#
#           diff(ci, xj) = diff(cj, xi) for all i<>j in {1, ..., n}.
# Remark: if H = 0, we set [] its list representation.
#
###############################################################################

###############################################################################
# Any rational function f(x) can be decomposited into
#
#            f = diff(p, x)/p + q/r,
#
# where p in K(x) and 
#
#            gcd(r, q-t*diff(r, x)) = 1 for all t in Z.
#
# Here, we call [p, q, r] a Gosper form of f(x),
#
# Multiplicative decomposition:
# Let H(x1, ..., xn) be a hyperexponential function of x1, ..., xn, with certificate
#
#              [[c1, x1], ..., [cn, xn]],
#
# Let x in {x1, ..., xn} and f be the certificate of H with respect to x.
# Let [p, q, r] be a Goaper form of f. Then H can be written as
#
#                               H = p * H', -----------------------  (**)
#
# where H' is a hyperexponential function with certificate
#
#      [[c1-diff(p, x1)/p, x1], ..., [q/r, x], ...,  [cn-diff(p, xn)/p, xn]].
#
# Here, we call [p, H'] a Multiplicative decomposition of H, and 
# call p a rational factor of H.    
###############################################################################


################################################################################
# Export
# Name: FromFunctionToList
# Calling sequence:
#      FromFunctionToList(h, varlist)
# Input: h, a hyperexponential function
#       varlist, a list of variable
# Output: L, the list representation of h
# Remark: Hyperexponential functions are considered as the same one if they
#        differ from a nonzero multiplicative constant.
################################################################################

FromFunctionToList := proc(h, varlist)
    local T, l, i, t;
    if normal(h) = 0 then
       return([])
    end if;
    T := table([]);
    l := nops(varlist);
    for i from 1 to l do
        t := varlist[i];
        T[i] := [normal(diff(h, t)/h), t];
    end do;
    convert(T, list);
end proc; 


##############################################
# Export
# Name: Certificates
# Calling sequence:
#         Certificates(H)
# Input: H, a list [[c1, x1], ..., [cn, xn]].
# Output: Cer, a list [c1, ..., cn].
##############################################
Certificates := proc(H)
     [seq(H[i][1], i =1..nops(H))];
end proc;

##############################################
# Export
# Name: Variables
# Calling sequence:
#         Variables(H)
# Input: H, a list [[c1, x1], ..., [cn, xn]].
# Output: Var, a list [x1, ..., xn].
##############################################
Variables := proc(H)
     [seq(H[i][2], i =1..nops(H))];
end proc;


#####################################################################
# Export
# Name: AreCompatible
# Calling sequence:
#         AreCompatible(H)
# Input: H, a list [[c1, x1], ..., [cn, xn]].
# Output: true if ck's are compatible, otherwise, false.
#####################################################################
AreCompatible := proc(H)
   local n, i, j, tst;
   if H = [] then
      return(true)
   end if;
   n := nops(H);
   for i from 1 to n do
       for j from i+1 to n do
           tst := diff(H[i][1], H[j][2])-diff(H[j][1], H[i][2]);
           if Testzero(tst) = 'false' then
              return false
           end if;
       end do;
   end do;
   true;
end proc;

####################################################################
# Export 
# Name: IsConstant
# Calling sequence:
#            IsConstant(H)
# Input: H, a list [[c1, x1], ..., [cn, xn]].
# Output: true if H is a constant, otherwise, false.
#####################################################################
IsConstant := proc(H)
      local n, i;
      if H = [] then
         return(true)
      end if; 
      n := nops(H);
      for i from 1 to n do
          if normal(H[i][1]) <> 0 then
             return(false)
          end if;
      end do;
      true;
end proc;

####################################################################
# Local 
# Name: IsPolynomial_1
# Calling sequence:
#            IsPolynomial_1(f, x, p)
# Input: f, a rational function;
#        x, a name;
# Output: true if exp(int(f, x)) is a polynomial, otherwise, false.
# Option: the second argument, when present, is assigned the usual
#         polynomial form of exp(int(f, x)) if the output is true. 
#####################################################################
IsPolynomial_1 := proc(f, x, p)
  local Rcf, fp, df, S, K;
  fp := normal(f);
  df := denom(fp);
  if degree(gcd(df, diff(df, x)), x) <> 0 then
     return(false);
  else
     Rcf := HyperexpReduction:-RationalCanonicalForm(fp, x);
     if Rcf[2] <> 0 or degree(denom(Rcf[1]), x) <> 0 then
        return(false)
     else
        if nargs = 3 then
           p := Rcf[1];
        end if;
        true;
     end if;
  end if;

end proc;


####################################################################
# Export 
# Name: IsPolynomial
# Calling sequence:
#            IsPolynomial(H, p)
# Input: H, a list [[c1, x1], ..., [cn, xn]].
# Output: true if H is a polynomial, otherwise, false.
# Option: the second argument, when present, is assigned the usual
#         functional form of H if the output is true. 
#####################################################################
IsPolynomial := proc(H, p)
      local  n, t, p1, p2, p3, vars, p3L; 
      if H = [] then
         return(true);
         if nargs = 2 then
             p := 0;
         end if;
      end if; 
      n := nops(H);

      # recursive base
       if n = 1 then
          t := IsPolynomial_1(op(H[1]), p1);
          if t = 'false' then
               return 'false';
          end if;
          if nargs = 2 then
              p := p1;
          end if;
          return 'true';
       end if;

       # recursion
       t := IsPolynomial([op(1..n-1, H)], p2);
       if t = 'false' then
          return 'false';
       end if;

     # combine
       t := IsPolynomial_1(normal(H[n][1]-diff(p2, H[n][2])/p2),  H[n][2],  p3);
       if t = 'false' then
           return 'false';
       end if;
       vars := [seq(H[i][2], i=1..n-1)];
       p3L := FromFunctionToList(p3, vars); 
       if not IsConstant(p3L) then
          return 'false';
       end if;
       if nargs = 2 then
          p := p3*p2;
       end if;
       true;

end proc;

####################################################################
# Local 
# Name: IsRational_1
# Calling sequence:
#            IsRational_1(f, x, p)
# Input: f, a rational function;
#        x, a name;
# Output: true if exp(int(f, x)) is a rational, otherwise, false.
# Option: the second argument, when present, is assigned the usual
#         rational form of exp(int(f, x)) if the output is true. 
#####################################################################
IsRational_1 := proc(f, x, p)
  local Rcf, fp, df, S, K;
  fp := normal(f);
  df := denom(fp);
  if degree(gcd(df, diff(df, x)), x) <> 0 then
     return(false);
  else
     Rcf := HyperexpReduction:-RationalCanonicalForm(fp, x);
     if Rcf[2] <> 0 then
        return(false)
     else
        if nargs = 3 then
           p := Rcf[1];
        end if;
        true;
     end if;
  end if;

end proc;

#####################################################################
# Export
# Name: IsRational
# Calling sequence:
#         IsRational(H, p)
# Input: H, a list [[c1, x1], ..., [cn, xn]].
# Output: true if H represents a rational function, otherwise, false.
# Option: the second argument, when present, is assigned the usual
#         functional form of H if the output is true. 
#####################################################################
IsRational := proc(H, p)
      local  n, t, p1, p2, p3, vars, p3L; 
      if H = [] then
         return(true);
         if nargs = 2 then
             p := 0;
         end if;
      end if; 
      n := nops(H);

      # recursive base
       if n = 1 then
          t := IsRational_1(op(H[1]), p1);
          if t = 'false' then
               return 'false';
          end if;
          if nargs = 2 then
              p := p1;
          end if;
          return 'true';
       end if;

       # recursion
       t := IsRational([op(1..n-1, H)], p2);
       if t = 'false' then
          return 'false';
       end if;

     # combine
       t := IsRational_1(normal(H[n][1]-diff(p2, H[n][2])/p2),  H[n][2],  p3);
       if t = 'false' then
           return 'false';
       end if;
       vars := [seq(H[i][2], i=1..n-1)];
       p3L := FromFunctionToList(p3, vars); 
       if not IsConstant(p3L) then
          return 'false';
       end if;
       if nargs = 2 then
          p := p3*p2;
       end if;
       true;

end proc;


#####################################################################
# Export
# Name: AreSimilar
# Calling sequence:
#         AreSimilar(H1, H2, opt)
# Input: H1, H2, two hyperexponential functions in list representation;
# Output: true if H1, H2 are similar, otherwise, false.
# Optional: If the third optional argument is presented, it is assigned
#           to the ratio of H1 and H2.
#####################################################################
AreSimilar := proc(H1, H2, opt)
   local  r, Hq, t;
   if H1 = [] or H2 = [] then
      return(true);
      if nargs = 3 then
         opt := 0;
      end if;
   end if;  
   Hq := HyperexpQuotient(H1, H2);
   if nargs = 3 then
      t := IsRational(Hq, 'r');
      opt := r;  
      return(t);     
   else
      IsRational(Hq);
   end if

end proc;
    
 



#####################################################################
# Export prodedure
# Name: HyperexpProduct
# Calling sequence:
#                  HyperexpProduct(H1, H2)
# Input: H1, H2, two hyperexponential functions in list representation;
# Output: H, a hyperexponential function in list representation such that
#                 H = H1 * H2.
#
#####################################################################
HyperexpProduct := proc(H1, H2)
     local H, nL, b, n1, n2, i, j, k;
     if H1 = [] or H2 = [] then
        return([])
     end if;
     H := [];
     nL := {};
     n1 := nops(H1);
     n2 := nops(H2);
     for i from 1 to n1 do  
         b := true; 
         for j from 1 to n2 do
             if b and evalb(H1[i][2] = H2[j][2]) then
                H := [op(H), [normal(H1[i][1]+H2[j][1]), H1[i][2]]];
                nL := {op(nL), j};
                b := false;
             end if;
         end do;
         if b then
            H := [op(H), H1[i]];
         end if;
     end do;
     nL := convert({seq(k, k=1..n2)} minus nL, list);
     H := [op(H), seq(H2[nL[k]], k=1..nops(nL))];

end proc;

#####################################################################
# Export prodedure
# Name: HyperexpQuotient
# Calling sequence:
#                  HyperexpQuotient(H1, H2)
# Input: H1, H2, two hyperexponential functions in list representation;
# Output: H, a hyperexponential function in list representation such that
#                 H = H1 / H2.
#
#####################################################################
HyperexpQuotient := proc(H1, H2)
     local H, nL, b, n1, n2, i, j, k;
     if H2 = [] then
        error "numeric exception: division by zero";
     end if;
     if H1 = [] then
        return([])
     end if;

     H := [];
     nL := {};
     n1 := nops(H1);
     n2 := nops(H2);
     for i from 1 to n1 do  
         b := true; 
         for j from 1 to n2 do
             if b and evalb(H1[i][2] = H2[j][2]) then
                H := [op(H), [normal(H1[i][1] - H2[j][1]), H1[i][2]]];
                nL := {op(nL), j};
                b := false;
             end if;
         end do;
         if b then
            H := [op(H), H1[i]];
         end if;
     end do;
     nL := convert({seq(k, k=1..n2)} minus nL, list);
     H := [op(H), seq([-H2[nL[k]][1], H2[nL[k]][2]], k=1..nops(nL))];

end proc;





#####################################################################
# Export
# Name: HyperexpDiff
# Calling sequence:
#         HyperexpDiff(H, x)
# Input: H, a list [[c1, x1], ..., [cn, xn]];
#        x, a name;
# Output: Hp, a list which represents the derivative of H with respect
#             to x.
#####################################################################
HyperexpDiff := proc(H, x)
     local Hp, n, i, Hr, vars;
     if IsConstant(H) then
         return([])
     end if;
     n := nops(H);
     vars := {seq(H[i][2], i=1..n)};
     if not evalb(x in vars) then
        return([])
     end if;
     for i from 1 to n do
         if evalb(H[i][2] = x) then
            Hr := FromFunctionToList(H[i][1], convert(vars, list));
            Hp := HyperexpProduct(Hr, H);
            return(Hp);
         end if;
     end do;

end proc;



#####################################################################
# Export
# Name: HyperexpSum
# Calling sequence:
#         HyperexpSum(H1, H2)
# Input: H1, H2, two hyperexponential functions which are similar.
# Output: Hp, the sum of H1 and H2 in list representation.
#####################################################################
HyperexpSum := proc(H1, H2)
    local t, Hp, r, vars, Hr;
    if H1 = [] then
       return(H2);
    else
       if H2 = [] then
          return(H1)
       end if;
    end if;
    t := AreSimilar(H1, H2, 'r');
    if t = 'false' then
       error "The input hyperexponential functions are not similar";
    end if;
    vars := convert({op(Variables(H1))} union {op(Variables(H2))}, list);
    Hr := FromFunctionToList(r+1, vars);
    Hp := HyperexpProduct(Hr, H2);

end proc;

######################################################################
# Export
# Name: AreTheSame
# Calling sequence:
#          AreTheSame(H1, H2)
# Input: H1, H2, two hyperexponential functions in list representation;
# Output: true, if H1, H2 represents the same function, otherwise, false.
######################################################################
AreTheSame := proc(H1, H2)
   local n, i, t, idx, S;
   if H1 = [] and H2 = [] then
      return(true);
   else if (H1 = [] and H2 <> []) or (H1 <> [] and H2 = []) then
           return(false)
        end if;
   end if;
   n := nops(H1);
   idx := [];
   for i from 1 to n do 
       t := FindCertificate(H2, H1[i][2]);
       if t[2] = [] then
          if normal(H1[i][1]) <> 0 then
             return(false)           
          end if;
       else
          idx := [op(idx), t[2]];
          if not Testzero(t[1]-H1[i][1]) then
             return(false)
          end if;   
       end if;             
   end do;
   S := [op({seq(i, i=1..nops(H2))} minus {op(idx)})];
   for i from 1 to nops(S) do
       if not Testzero(H2[S[i]][1]) then
          return(false);
       end if;
   end do;
   true;  

end proc;

######################################################################
# Export
# Name: TakeVariables
# Calling sequence:
#         TakeVariables(H)
# Input: H, a hyperexponential function in list representation;
# Output: V, a list of variables of H.
######################################################################
#TakeVariables := proc(H)
#    local L, n;
#    n := nops(H);
#    L := [seq(H[i][2], i=1..n)];
#end proc:
######################################################################
# Export
# Name: FindCertificate
# Calling sequence:
#         FindCertificate(H, x)
# Input: H, a hyperexponential function in list representation;
#        x, a name;
# Output: c, the certificate of H w.r.t. x;
#         p, the position of c in H.
######################################################################
FindCertificate := proc(H, x)
    local n, i;
    n := nops(H);
    for i from 1 to n do
        if H[i][2] = x then
           return([H[i][1], i])
        end if;
    end do;
    [0, []];
end proc;


end module:


HyperexpReduction := module()
    description "Computing Hermite reduction of a hyperexponential function";
    
    export
          AreCompatible,
          HermiteLikeReduction,       
          Original,
          RatExpReduction,
          RationalCanonicalForm,
          SolvingGosperEquation,
          UndeterminedCoefficients,
          VerifyHLR,
          VerifyOriginal,
          VerifyRCF;
    local 
          BezoutRelation,
          FactorList,
          myExtendedGcd,
          MultOneSimpleFraction,
          ParFrac,
          ReduceOneComponent,
          SimpleReduce,
	  SolvingGosperEquation2,
          UpperBound,
          VerifyReduceOneComponent;
                     
    option package;



###############################################################################
#                                                                             #
#               Export procedures of HyperexpReduction module                 #
#                                                                             #
###############################################################################

###############################################################################
# Remark. 
# For testing this module, one should include the following modules:
#     Division.mm
#     Decomposition.mm
#     List.mm
###############################################################################
  
###############################################################################
# Export  
# Name: AreCompatible
# Calling sequence:
#         AreCompatible(H)
# Input: H, a list [[c1, x1], ..., [cn, xn]].
# Output: true if ck's are compatible, otherwise, false.
###############################################################################
AreCompatible := proc(H)
   local n, i, j, tst;
   if H = [] then
      return(true)
   end if;
   n := nops(H);
   for i from 1 to n do
       for j from i+1 to n do
           tst := diff(H[i][1], H[j][2])-diff(H[j][1], H[i][2]);
           if Testzero(tst) = 'false' then
              return false
           end if;
       end do;
   end do;
   true;
end proc;

###############################################################################
# Export
# Name: RatExpReduction
# Calling sequence:
#                  RatExpReduction(f, H, x)
# Input: f, a rational function of x1, ..., xn;
#        H, a function of the form exp(g), with g a rational function; 
#        x, a name in {x1, ..., xn};
# Output: [f1, H1], [f2, H2], two lists of functions such that
#                  f*H = diff(f1*H1, x) + f2*H2, 
#         where H1, H2 are functions of the form exp(g), with g a rational function,
#         and f2 has squarefree denominator. 
#         
################################################################################
RatExpReduction := proc(f, H, x)
       local S, K, tt, Rd, h, p, L1, L2, r1, f1, H1,
             k1, k2, g1, g2, i, r, v, u, j, r2, f2, H2, sol, qv;
       
       if Testzero(f) then
          return([[], []])
       end if;

       # Compute the rational canonical form of f*H
         S := f;
         K := normal(diff(H, x)/H);

       # Compute Hermite-like reduction of (S, K)
         #print(S);
         #print(K);
         tt := time():
         print(begin_HLR);
         Rd := Original(S, K, x);
         (h, p, L1, L2) := Rd[1], Rd[2], Rd[3], Rd[4];
         tt := time()-tt;
         print(time_for_HLR);
         print(tt);
          
        # Preapre for checking integrability
          tt := time():
          print(begin_prepare_for_check_integrability);
          r1 := List:-FromListToNormal(h);
          if r1 = 0 then
             f1 := 0;
             H1 := H;
          else
             f1 := r1;
             H1 := H;
          end if;
          (k1, k2) := numer(K), denom(K);
          g1 := 0;
          g2 := 0;
          for i from 1 to nops(L1) do 
              g1 := normal(g1 + L1[i][2][1]*L1[i][2][2]/L1[i][1]);
          end do;
          for i from 1 to nops(L2) do 
              g2 := normal(g2 + L2[i][2][1]*L2[i][2][2]/(L2[i][1]));
          end do;
          r := normal(p*k1+g1*k2+g2);
          (v, u) := numer(r), denom(r);
          if rem(v, k2, x, 'qv') = 0 then
             j := 0;
             v := qv;
          else
             j := 1;
          end if;
          r2 := [v, u, k2, j];
          tt := time()-tt;
          print(time_for_prepare_for_check_integrability);
          print(tt);

           # Check integrability
             tt := time():
             print(begin_check_integrability);
             if degree(u, x) > 0 then
                f2 := v/(u*k2^j);
                H2 := H;
                return([[f1, H1], [f2, H2]])
             end if; 
             
	     #-------------------------------------
             # Solve the Gosper equation 
             #       ry'+(r'+q)y = p  
             #-------------------------------------
              sol := SolvingGosperEquation(normal(v/u), normal(k1-j*diff(k2, x)), k2, x);
              print(time_for_checking_integrability);
              print(time()-tt);

              if evalb(sol = []) then
                 f2 := r2;
                 H2 := H;
                 return([[f1, H1], [f2, H2]])
              else
                 f1 := normal(f1 + normal(sol*k2/(k2^j)));
                 return([[f1, H], []]);
              end if;

end proc;

###############################################################################
# Export
# Name: HermiteLikeReduction
# Calling sequence:
#                  HermiteLikeReduction(H, x)
# Input: H, a hyperexponential function of x1, ..., xn in list representation;
#        x, a name in {x1, ..., xn};
# Output: H1, H2, two hyperexponential function such that
#                  H = diff(H1, x) + H2, 
#         where H1, H2 are represented as
#                 [r1, T1] and [r2, T2]
#         with r1, r2 in K(x1, ..., xn), and T1, T2 in list representation such that
#
#                 H1 = r1 * T1 and H2 = r2 * T2.
# Remark. r2 is a list [v, u, k2, i], where i is in {0, 1} and this list 
#         representing the rational function v/(u*k2)^i.         
################################################################################
HermiteLikeReduction := proc(H, x)
       local f, Ec, b, n, i, Rcf, Rd, h, p, L1, L2, S, K, k1, k2, 
             v, u, qv, j,  r1, r2, r, f1, f2, H1, H2, tt, solb, sol;
       
       if H = [] then
          return([[], []])
       end if;

       # Check the compatibility
        # if not AreCompatible(H) then
        #    error "The input certificates are not compatible";
        # end if;

       # Find the certificate of H w.r.t. x
         n := nops(H);
         b := true;
         for i from 1 to n do
             if b and evalb(H[i][2] = x) then 
                f := H[i][1];
                b := false;
             end if; 
         end do;
         if b then
            H1 := [x, H];
            H2 := [];
            return([H1, H2]);
         end if;
		 

        # Compute the rational canonical form of f
          tt := time():
          print(begin_RCF);
          Rcf := RationalCanonicalForm(f, x);
		  #print(Verifying_RCF);
		  #print(VerifyRCF(Rcf, f, x));
          tt := time()-tt:
          print(time_for_RCF);
          print(tt);
          (S, K) := Rcf[1], Rcf[2];
           
        # Compute the multiplicative decomposition of H
          tt := time():
          print(begin_MulDecomposition);
          Ec := [];
          for i from 1 to n do
              if evalb(H[i][2] = x) then
                 Ec := [op(Ec), [K, x]];
              else
                 Ec := [op(Ec), [normal(H[i][1]-diff(S, H[i][2])/S), H[i][2]]];
              end if;
          end do;
          tt := time()-tt:
          print(time_for_MulDecomposition);
          print(tt);
         
         # Compute Hermite-like reduction of (S, K)
           tt := time():
           #print(begin_HLR);
           Rd := Original(S, K, x);
           (h, p, L1, L2) := Rd[1], Rd[2], Rd[3], Rd[4];
	   #print(Verifying_Original);
	   #print(VerifyOriginal(Rd, S, K, x));
           tt := time()-tt;
           print(time_for_HLR);
           print(tt);
          
 
         # Preapre for checking integrability
           tt := time():
           print(begin_prepare_for_check_integrability);
           r1 := List:-FromListToNormal(h);
           if r1 = 0 then
              H1 := [];
           else
              H1 := [r1, Ec];
           end if;
           (k1, k2) := numer(K), denom(K);
           f1 := 0;
           f2 := 0;
           for i from 1 to nops(L1) do 
               f1 := normal(f1 + L1[i][2][1]*L1[i][2][2]/L1[i][1]);
           end do;
           for i from 1 to nops(L2) do 
               f2 := normal(f2 + L2[i][2][1]*L2[i][2][2]/(L2[i][1]));
           end do;
           r := normal(p*k1+f1*k2+f2);
           (v, u) := numer(r), denom(r);
           if rem(v, k2, x, 'qv') = 0 then
              j := 0;
              v := qv;
           else
              j := 1;
           end if;
           r2 := [v, u, k2, j];
           tt := time()-tt;
           print(time_for_prepare_for_check_integrability);
           print(tt);

           # Check integrability
             tt := time():
             print(begin_check_integrability);
             if degree(u, x) > 0 then
                H2 := [r2, Ec];
                return([H1, H2])
             end if; 
             
	     #-------------------------------------
             # Solve the Gosper equation 
             #       ry'+(r'+q)y = p  
             #-------------------------------------
              sol := SolvingGosperEquation(normal(v/u), normal(k1-j*diff(k2, x)), k2, x);
              print(time_for_checking_integrability);
              print(time()-tt);

              if evalb(sol = []) then
                 H2 := [r2, Ec];
                 return([H1, H2]);
              else
                 r1 := normal(r1 + normal(sol*k2/(k2^j)));
                 return([[r1, Ec], []]);
              end if;

           
end proc;

###############################################################################
# Export  
# Name: VerifyHLR
# Calling sequence:
#                 VerifyHLR(R, H, x)
# This function is for verifying the correctness of the previous function
#                  HermiteLikeReduction(R, H, x)
#
###############################################################################
VerifyHLR := proc(R, H, x)
      local Hp, r1, r2, r, H1, H2, H3, H4, Hr, vars, tt; 
      if R[1] = [] and R[2] = [] then
         if H = [] then
            return(true)
         else
            return(false)
         end if;
      else if R[1] = [] then
              vars := Hyperexp:-Variables(R[2][2]);
              r2 := normal(R[2][1][1]/(R[2][1][2]*R[2][1][3]^R[2][1][4]));
              H2 := Hyperexp:-FromFunctionToList(r2, vars);
              H4 := Hyperexp:-HyperexpProduct(H2, R[2][2]);
              return(Hyperexp:-IsConstant(Hyperexp:-HyperexpQuotient(H, H4)));
           end if;
      end if; 
  
      vars := Hyperexp:-Variables(R[1][2]);
      r1 := R[1][1];
      if R[2]= [] then
         H1 := Hyperexp:-FromFunctionToList(r1, vars);
         Hp := Hyperexp:-HyperexpDiff(Hyperexp:-HyperexpProduct(H1, R[1][2]), x);
         return(Hyperexp:-AreTheSame(H, Hp));
      else
        
         r2 := normal(R[2][1][1]/(R[2][1][2]*R[2][1][3]^R[2][1][4]));
         r := normal(diff(r1, x) + Hyperexp:-FindCertificate(R[1][2], x)[1]*r1+r2);
         Hr := Hyperexp:-FromFunctionToList(r, vars);
         Hp := Hyperexp:-HyperexpProduct(Hr, R[1][2]);
         return(Hyperexp:-AreTheSame(H, Hp));
      end if;  

end proc;


###############################################################################
# Note. Modified on 21/09/2009
# Export  
# Name: SolvingGosperEquation
# Calling sequence:
#         SolvingGosperEquation(p, q, r, x)
# Input: p, q, r, three polynomials with 
#                 gcd(r, q-t*diff(r, x)) = 1 for all t in Z;
#        x, a name;
# Output: S, a polynomial if there exists a polynomial solution of the equation
#                     r*diff(y, x) + (diff(r, x)+q)*y = p;
#         [], for otherwise.
###############################################################################
SolvingGosperEquation2 := proc(p, q, r, x)
       local sol, Data;
      sol := DETools[polysols]([normal(q+diff(r, x)), r], p, x);       
      if evalb(sol = []) then
         return([]);
      else
          sol[2];
      end if;
       
end proc;

SolvingGosperEquation := proc(p, q, r, x)
       local sol;
       sol := UndeterminedCoefficients(r, normal(q+diff(r, x)), p, x);
       
end proc;




###############################################################################
# Export
# Name: UndeterminedCoefficients
# Calling sequence:
#         UndeterminedCoefficients(a, b, c, x, opt)
# Input: a, b, c, three polynomials;
#        x, name;
# Output: A polynomial S in K[x] solving 
#            a*diff(y, x) + b*y = c,      
#         or [], which means there is no polynomial solution.
# Optional: If the fifth argument is presented, then assign it to the upper bound of polynomial
#           solutions.
###############################################################################
UndeterminedCoefficients := proc(a, b, c, x, opt)
       local g, ap, bp, cp, Bd, sol, i, s, S, eqns, vars, qq;
       if Testzero(a) then 
          if Testzero(b) then
             if Testzero(c) then
                return(1);
             else
                return([]);
             end if;
          else
             if Testzero(rem(c, b, x)) then
                return(normal(c/b));
	     else 
		return([]);
	     end if;
          end if;
       else 
          if Testzero(b) then
             if Testzero(rem(c, a, x)) then
                 return(int(c/a, x));
             else
                 return([]);
             end if;
          end if;    
       end if;
 
       g := gcd(a, b);
       if rem(c, g, x) <> 0 then
          return([])
       end if;
       ap := Division:-ExactDivision(a, g, 1, x);
       bp := Division:-ExactDivision(b, g, 1, x);
       cp := Division:-ExactDivision(c, g, 1, x);
       
       if nargs = 5 then
          Bd := args[5];
       else
          Bd := UpperBound(ap, bp, cp, x);
       end if;

       if Bd < 0 then
	      if degree(b, x) = degree(c, x) and rem(c, b, x, 'qq') = 0 then
		     return(qq);
	      else
             return([]);
		  end if;
       end if;
	   
       s := table([]);
       S := add(s[i]*x^i, i=0..Bd);
       eqns := PolynomialTools:-CoefficientList(collect(ap*diff(S, x)+bp*S-cp, x, distributed), x);
       vars := [seq(s[i], i=0..Bd)];
       sol := SolveTools:-Linear(eqns, vars);
       if sol = NULL then
          return([])
       else 
          S := eval(S, sol);
       end if;
         
end proc;





###############################################################################
# Export  
# Name: RationalCanonicalForm
# Calling sequence:
#         RationalCanonicalForm(f, x)
# Input: f, a rational function in x;
#        x, a name;
# Output: S, K, two rational function such that
#               f = K + diff(S, x)/S,
#               where K contains all simple fractions of class II.
###############################################################################
RationalCanonicalForm := proc(f, x)
     local R, a, d, L, g, ca, cd, P, n, i, S, K, m;
       
     if type(normal(f), freeof(x)) then
        return([1, f])
     end if; 
     R := MultOneSimpleFraction(f, x);
     (a, d, g) := R[1][1], R[1][2], R[2];
     if a = 0 then
        return([1, g])
     end if;
     L := FactorList(d, x);
     cd := L[1];
     P := ParFrac(L[2], primpart(a, x, 'ca'), x);
     n := nops(P);
     S := 1;
     K := g;
     for i from 1 to n do
         m := normal(ca*P[i][2]/cd/diff(P[i][1], x));
         if type(m, integer) then
            S := S*P[i][1]^m;
         else
            K := K + ca*P[i][2]/P[i][1]/cd;
         end if;   
     end do;
     [S, K];

end proc;


###############################################################################
# Export  
# Name: VerifyRCF
# Calling sequence:
#                 VerifyRCF(R, f, x)
# This function is for verifying the correctness of the previous function
#                   RationalCanonicalForm(f,x)
#
###############################################################################
VerifyRCF := proc(R, f, x)
      Testzero(f-diff(R[1], x)/R[1]-R[2]);
end proc;

###############################################################################
# Export  
# Name: original
# Calling sequence:
#                  Original(f, g, x)
# Input: f, g, two rational functions written as
#                      f = a/b, g = s/t,
#              such that gcd(b, t) = 1 and g is differential reduced, i.e.
#                gcd(t, s-i*diff(t, x)) = 1 for all i in Z;
#        x, a variable name.
# Output: h,  a rational function(in list representation);
#         P,  a polynomial; 
#         L1,  a list [[u1, [c1, v1]], ..., [um, [cm, vm]]] with ui sqrfree and ci free of x;
#         L2,  a list [[z1, [e1, w1]], ..., [zp, [ep, wp]]] with zi sqrfree and ei free of x;
#              such that
#               f = diff(h, x) + h*g + P*g + c1*v1/u1 + ... + cm*vm/um 
#                                          + e1*w1/(z1*t) + ... + em*wp/(zp*t).
#
###############################################################################

Original := proc(f, g, x)
   local pf, c, S, n, i, h, hp, P, L, L1, L2, v, u, j, R;
   
   # Computing squarefree partial fraction decomposition of f    
    c, S := Decomposition:-SquareFreePartialFraction(normal(f), x);
    n := nops(S[3]);

   # Case: f is polynomial
     if S[3] = [] then
        return([[x, 0, []], 0, [[1, [content(f, x, 'pf'), pf]]], [[1, [1, 0]]]])
     end if;

   # Prepare for reduction
     h := [x, int(S[2], x), []];
     P := -c*h[2];
     L1 := [];
     L2 := [];
     hp := [];

   # Process of reduction
     for i from 1 to n do
         L := S[3][i];
         R := ReduceOneComponent(L, g, x);
         hp := [op(hp), R[1]];
         L1 := [op(L1), [R[2][1], [R[2][2][1]*c, R[2][2][2]]]];
         L2 := [op(L2), [R[3][1], [R[3][2][1]*c, R[3][2][2]]]];
     end do;
     
   # Prepare for return
     h := List:-ScalarMultiply(c, [x, h[2], hp]);
     [h, P, L1, L2]; 
    
end proc;

###############################################################################
# Export  
# Name: VerifyOriginal
# Calling sequence:
#                 VerifyOriginal(R, f, g ,x)
# This function is for verifying the correctness of the previous function
#                   Original(f, g ,x)
#
###############################################################################
VerifyOriginal := proc(R, f, g ,x)
      local hL, hF, P, L1, L2, f1, f2, i, t;
      (hL, P, L1, L2) := seq(R[i], i=1..nops(R));
      t := denom(g);
      hF := List :-FromListToNormal(hL);
      f1 := 0;
      f2 := 0;
      for i from 1 to nops(L1) do 
          f1 := normal(f1 + L1[i][2][1]*L1[i][2][2]/L1[i][1]);
      end do;
      for i from 1 to nops(L2) do 
          f2 := normal(f2 + L2[i][2][1]*L2[i][2][2]/(L2[i][1]*t));
      end do;
      Testzero(diff(hF, x)+hF*g+P*g+f1+f2-f);

end proc;




###############################################################################
#                                                                             #
#                 Local procedures of HyperexpReduction module                #
#                                                                             #
###############################################################################

###############################################################################
# Local
# Name: UpperBound
# Calling sequence:
#                  UpperBound(a, b, c, x)
# Input: a, b, c, three polynomials;
#        x, a name;
# Output: k, a positive integer which is a upper bound of the degree of the 
#            polynomial solution of the following equation
#                       a*diff(y, x) + b*y = c.
#                   
###############################################################################
UpperBound := proc(a, b, c, x)
     local da, db, dc, k, k0;
     (da, db , dc) := degree(a, x), degree(b, x), degree(c, x);
     if db >= da then
        k := dc - db;
     else 
        if db < (da-1) then
           k := dc-da+1;
        else
           if db = da-1 then
              k0 := Normalizer(-lcoeff(b, x)/lcoeff(a, x));
              if type(k0, integer) then
                 k := max(k0, dc-da+1);
              end if;
           end if;
        end if
     end if;
     k;

end proc:




###############################################################################
# Local  
# Name: ParFrac
# Calling sequence:
#                  ParFrac(L, a, x)
# Input: L, a list of polynomials [p1, ..., pm] with gcd(pi, pj)=1 for i<>j;
#        a, a polynomial with gcd(pi, a)=1 for j=1..m;
#        x, a nanme;
# Output: Lp, a list [[p1, a1], ..., [pm, am]] such that
#               a/(p1*...*pm) = a1/p1 + ... + am/pm.
#
###############################################################################
ParFrac := proc(L, a, x)
   local n, i, j, Lp, ap, q, R;
   
   n := nops(L);
   if n = 1 then
      return([[L[1], a]])
   end if;
   Lp := [];
   ap := a;
   for i from 1 to n do
       q := mul(L[j], j=i+1..n);
       R := myExtendedGcd(L[i], q, ap, x);
       ap := R[1];
       Lp := [op(Lp), [L[i], R[2]]];
   end do;
   Lp;
end proc;


###############################################################################
# Local  
# Name: MultOneSimpleFraction
# Calling sequence: 
#                 MultOneSimpleFraction(f, x) 
# Input: f, a rational function;
#        x, a name;
# Output: [a, d], a list of polynomials where d is the product of factors of f with
#                 multiplicity 1;
#         g,      a rational function such that
#                          f = a/d + g.
# 
###############################################################################   
MultOneSimpleFraction := proc(f, x) 
     local nu, de, re, qu, cd, cr, ar, pd, co, R, s, t, a, d, h, d1, d2, g;
    
     (nu, de) := numer(f), denom(f);
     re := rem(nu, de, x, 'qu');
     if re = 0 then
        return([0, 1], qu);
     end if;
 
     cd := content(de, x, 'pd');
     h := gcd(pd, diff(pd, x), 'co');
     gcd(co, h, 'd1');
     divide(pd, d1, 'd2');
     cr := content(re, x, 'ar');
     R := myExtendedGcd(d1, d2, ar, x);
     (s, t) := R[1], R[2];
     a := cr*t/cd;
     g := qu+(cr*s/(cd*d2));
 
     [[a, d1], g];
 
end proc:

    
    

###############################################################################
# Local  
# Name: FactorList
# Calling sequence:
#                  FactorList(p, x)
# Input: p, a squarefree polynomial; 
#        x, a name;
# Output: [c, [p1, ..., pm]], a list with c free of x and pi is irreducible 
#                             factor of p and p = c*p1*...*pm.
###############################################################################
FactorList := proc(p, x)
      local c, fp, L1, L2, pp, i;
      c := content(p, x, 'pp');
      fp := factor(pp);
      if type(fp, `*`) then
         L1, L2 := selectremove(type, [op(fp)], freeof(x));
         c := c*mul(L1[i], i=1..nops(L1));
      else 
         L2 := [fp]; 
           #if type(fp, `+`) then
           #   L2 := [fp];
           #end if;
      end if;
      [c, L2];
end proc;



###############################################################################
# Local  
# Name: BezoutRelation
# Calling sequence:
#         BezoutRelation(P, Q, R, s, x)
# Input: P, Q, two coprime polynomails in x;
#        R, a list of polynomials;
#        s, a polynomial such that s*p = 1 mod Q;
#        x, a name;
# Output: B, C, two polynomials such that
#               B*P + C*Q = mul(R[i], i=1..nops(R)).   
#
###############################################################################
BezoutRelation := proc(P, Q, R, s, x)
      local B, C, RR, i;
      B := Division:-MergeRemainder([s, op(R)], Q, x);
      RR := mul(R[i], i=1..nops(R));
      C := Division:-ExactDivision(normal(RR-B*P), Q, 1, x);
      [B, C];
end proc; 

###############################################################################
# Local  
# Name: myExtendedGcd
# Calling sequence:
#          myExtendedGcd(p, q, a, x)
# Input: p, q, a, three polynomials with gcd(p, q) = 1 and deg(a, x) < deg(pq, x);
#        x, a name;
# Output: s, t, two polynomials such that
#               s*p + t*q = a.
#
################################################################################
myExtendedGcd := proc(p, q, a, x)
   local ind, S, c, H, C, L, Lp, n, i, s, t;
   
   ind := indets(p) union indets(q);
   if evalb(ind={x}) then
      gcdex(q, p, a, x, 't', 's');
      return([s, t]);
   end if;
   S := Decomposition:-SquareFree(q, x);
   c := S[1];
   n := nops(S[2]);
   L := [seq(S[2][i][1], i=1..n)];
   H := Division:-HalfEuclideanList(L, p, x);
   C := [a, seq(H[i]$S[2][i][2], i=1..n)];
   t := Division:-MergeRemainder(C, p, x); 
   s := Division:-ExactDivision(a-t*q/c, p, 1, x);
   [s, t/c]; 
 
end proc:
  
###############################################################################
# Local  
# Name: SimpleReduce
# Calling sequence: 
#                  SimpleReduce(L, i, K, C, x)
# Input: L, a list [c, a, d, m] representing a rational function
#           c*a/d^m with d squarefree;  
#        i, a integer in {0, 1};  
#        K, a list [k1, k2] representing a differential reduced 
#           rational function k1/k2 and gcd(d, k2) = 1;
#        C, a list [s0, s1] such that 
#                  s0 * d = 1 mod and s1 * d = 1 mod diff(d, x)*k2;
#        x, a name. 
# Output: L1, L2, L3, three lists of the form 
#                      [cj, aj, dj, mj], j = 1..3 such that
#                        c*a/(d^m*k2^i) = diff(S1, x) + S1*K + S2 + S3/k2,
#                      where Sj = cj*aj/(dj^mj) for j = 1..3.
###############################################################################  
SimpleReduce := proc(L, i, K, C, x)
      local c, a, d, j, dp, m, k1, k2, B, s0, s1, sp, tp, L1, L2, L3, pp1, pp2, pp3;

     # Initialize 
       (c, a, d, m) := L[1], L[2], L[3], L[4];
       (k1, k2) := K[1], K[2];
       if not evalb(i in {0, 1}) then
          error "the second argument should be in {0, 1}";
       end if;

     # Trivial Cases: m = 1 or a = 0
       if a = 0 then
          return([seq([1, 0, d, m-1], j=1..3)]);
       end if;
       if m = 1 then
          L1 := [1, 0, d, m];
          if i = 0 then
             L2 := [c, a, d, m];
             L3 := [1, 0, d, m];
          else
             L2 := [1, 0, d, m];
             L3 := [c, a, d, m];
          end if;
          return([L1, L2, L3]);
       end if;
 
     # Prepare for reduction
       dp := diff(d, x);
       (s0, s1) := C[1], C[2];

     # Case: i=0
       if i = 0 then
          B := BezoutRelation(d, dp, [a], s0, x);
          (sp, tp) := B[1], B[2];
          L1 := [c*content(tp/(1-m), x, 'pp1'), pp1, d, m-1];
          L2 := [c*content(normal(sp-diff(tp, x)/(1-m)), x, 'pp2'), pp2, d, m-1];
          L3 := [-L1[1]*content(k1, x, 'pp3'), L1[2]*pp3, d, m-1];
     # Case: i=1
       else
         B := BezoutRelation(d, dp*k2, [a], s1, x);
         (sp, tp) := B[1], B[2];
         L1 := [c*content(tp/(1-m), x, 'pp1'), pp1, d, m-1];
         L2 := [c*content(normal(-diff(tp, x)/(1-m)), x, 'pp2'), pp2, d, m-1];
         L3 := [c*content(sp-tp*k1/(1-m), x, 'pp3'), pp3, d, m-1];
       end if;
     
     return([L1, L2, L3]);
end proc;






###############################################################################
# Local prodecure
# Name: ReduceOneComponent
# Calling sequence:
#         ReduceOneComponent(L, K, x)
# Input: L, a list [d, [c1, a1], ..., [cn, an]];
#        K, a rational function k1/k2 such that 
#             gcd(k2, d) = 1;
#        x, a name;
# Output: S1, S2, S3, where S1 is a rational function in list representation
#                          [d, [e1, b1], ..., [em, bm]], with deg(ei, x)=0, 
#                     and S2, S3 are ratioanl functions in list representation
#                         [d, [cp2, ap2]] and [d, [cp3, ap3]] such that
#                     
#                     c1*a1/d+...+cn*an/d^n = diff(S1, x) + S1*K + S2 + S3/k2.     
#
################################################################################# 
ReduceOneComponent := proc(L, K, x)
      local d, dp, c, a, pp, k1, k2, KL, s0, s1, C, n, i, j, S1, S2, S3, L0, L1, SR0, SR1;
      
      n := nops(L)-1;
      d := L[1];

      # Case: n = 1
        if n = 1 then
           S1 := [d, [1, 0]];
           S2 := L;  
           S3 := [d, [1, 0]];
           return([S1, S2, S3]);   
        end if;

      # prepare for reduction
        dp := diff(d, x);
        k1 := numer(K);
        k2 := denom(K);
        KL := [k1, k2];
        gcdex(d, dp, x, 's0');
        gcdex(d, dp*k2, x, 's1');
        C := [s0, s1];

      # Reduction process
        S1 := [];
        S2 := [d, [1, 0]];
        S3 := [d, [1, 0]];
        for i from n by -1 to 2 do
            L0 := [content(S2[2][1]*S2[2][2]+L[i+1][1]*L[i+1][2], x, 'pp'), pp, d, i];
            L1 := [op(S3[2]), d, i];
            SR0 := SimpleReduce(L0, 0, KL, C, x);
            SR1 := SimpleReduce(L1, 1, KL, C, x);
            S1 := [[content(SR0[1][1]*SR0[1][2]+SR1[1][1]*SR1[1][2], x, 'pp'), pp], op(S1)];
            S2 := [d, [content(SR0[2][1]*SR0[2][2]+SR1[2][1]*SR1[2][2], x, 'pp'), pp]];
            S3 := [d, [content(SR0[3][1]*SR0[3][2]+SR1[3][1]*SR1[3][2], x, 'pp'), pp]];         
        end do;
          
      # prepare for return
        S1 := [d, op(S1)];
        S2 := [d, [content(S2[2][1]*S2[2][2]+L[2][1]*L[2][2], x, 'pp'), pp]];
        [S1, S2, S3];
 
 end proc;

###############################################################################
# Local  
# Name: VerifyReduceOneComponent
# Calling sequence: 
#                 VerifyReduceOneComponent(R, L, K, x)
# This function is for verifying the correctness of the previous function
#                 ReduceOneComponent(L, K, x)
###############################################################################

VerifyReduceOneComponent := proc(R, L, K, x)
      local nl, nr, i, le, re, S1, S2, S3;

      nl := nops(L)-1;
      le := normal(add(L[i+1][1]*L[i+1][2]/L[1]^i, i=1..nl));
      nr := nops(R[1])-1;
      S1 := normal(add(R[1][i+1][1]*R[1][i+1][2]/R[1][1]^i, i=1..nr));
      S2 := normal(R[2][2][1]*R[2][2][2]/R[2][1]);
      S3 := normal(R[3][2][1]*R[3][2][2]/R[3][1]);
      re := normal(diff(S1, x) + S1*K + S2 + S3/denom(K));
      Testzero(le-re);
     
end proc;
   
   
   
end module: 
   
   
   
   
   LinearRelation := module()
    description "Find linear dependency between polynomials";
    
    export
          FindLinearRelation,
          FindLinearRelation1,
          FindLinearRelation2;
    local       
          FindLinearRelationOneList, 
          IsZeroVector,
          IsZeroList,
          LowerBoundForRank;
                     
    option package;

###############################################
## Name: FindLinearRelation
##  
##           FindLinearRelation(pL, vL, y)
## Input: pL, vL, two same-size lists of polynomials;
##        y, a variable name.
## Output: cL, a list of constants (free of y) such that
##          
##         add(pL[i]*cL[i], i=1..n) =0 
##                   and
##         add(vL[i]*cL[i], i=1..n) =0
##
##              where n = nops(pL).
################################################

FindLinearRelation :=proc(pL, vL, y)
    
    local cL, tt, np, pp, clp, gm, MM, clv, pv, sys, vars, i, j, m, p, r, sol;

    cL := table([]);
    np := nops(pL);

    #### Generate a matrix of coefficients
    vars := [seq(cL[i], i=1..np)];
    pp   := collect(numer(add(cL[i]*pL[i], i=1..np)), y);
    pv   := collect(numer(add(cL[i]*vL[i], i=1..np)), y);


    clp  := convert(PolynomialTools:-CoefficientList(pp, y), set) minus {0};
    clv  := convert(PolynomialTools:-CoefficientList(pv, y), set) minus {0};  
    sys  := clp union clv;
    gm   := LinearAlgebra:-GenerateMatrix(sys, vars);
    if IsZeroVector(gm[2]) then
       MM := gm[1];
    else
       return([]);
    end if;
   
    #### Find linear dependency
    for j from 1 to 3 do
      m := rand(6..100):
      p := ithprime(m()):
      r := LowerBoundForRank(MM, p):
      if r = np then
         return([]):
      end if:
    end do: 
    ###tt := time():
    sol := LinearAlgebra:-NullSpace(MM);

    ####tt := time()-tt; print(tt);


    if nops(sol) =0 then
       return([]);
    else
       convert(sol[1], list);
    end if;

end proc;


FindLinearRelation1 :=proc(pL, vL, y)
    
    local cL, cpt, cvt, tt, np, pp, clp, gm, MM, clv, pv, sys, vars, i, j, m, p, r, sol;

    cL := table([]);
    np := nops(pL);

    #### Generate a matrix of coefficients
    vars := [seq(cL[i], i=1..np)];
    pp   := collect(numer(add(cL[i]*pL[i], i=1..np)), y);
    pv   := collect(numer(add(cL[i]*vL[i], i=1..np)), y);
    
    cpt := content(pp, [op(vars), y], 'pp');
    cvt := content(pv, [op(vars), y], 'pv');
    print([cpt, cvt]);
    

    clp  := convert(PolynomialTools:-CoefficientList(pp, y), set) minus {0};
    clv  := convert(PolynomialTools:-CoefficientList(pv, y), set) minus {0};  
    sys  := clp union clv;
    gm   := LinearAlgebra:-GenerateMatrix(sys, vars);
    if IsZeroVector(gm[2]) then
       MM := gm[1];
    else
       return([]);
    end if;
   
    #### Find linear dependency
    for j from 1 to 2 do
      m := rand(6..100):
      p := ithprime(m()):
      r := LowerBoundForRank(MM, p):
      if r = np then
         return([]):
      end if:
    end do: 
    tt := time():
    sol := LinearAlgebra:-NullSpace(MM);

    tt := time()-tt; print(tt);


    if nops(sol) =0 then
       return([]);
    else
       convert(sol[1], list);
    end if;

end proc;



FindLinearRelationOneList :=proc(pL, vars, y) 
 
    local cL, cpL, ppL, np, pp, poly, clp, gm, MM, sys,  i, j, m, p, r, sol;

    
    np := nops(pL);
    cpL := table([]);
    ppL := table([]);

    ### Find the contents of pL
    for i from 1 to np do
        cpL[i] := content(pL[i], y, 'pp');
        ppL[i] := pp; 
    end do;
 

    #### Generate a matrix of coefficients
    poly   := collect(normal(add(vars[i]*ppL[i], i=1..np)), y);
    sys  := convert(PolynomialTools:-CoefficientList(poly, y), set) minus {0};
    gm   := LinearAlgebra:-GenerateMatrix(sys, vars);
    if IsZeroVector(gm[2]) then
       MM := gm[1];
    else
       return([]);
    end if;
   
    #### Find linear dependency
    for j from 1 to 2 do
      m := rand(6..100):
      p := ithprime(m()):
      r := LowerBoundForRank(MM, p):
      if r = np then
         return([]):
      end if:
    end do: 

    ### sol := LinearAlgebra:-NullSpace(MM);
    sol := SolveTools:-Linear(sys, vars);
    if IsZeroList(eval(vars, sol)) then
       return([]);
    else
       ## convert(sol[1], list);
       eval([seq(vars[i]/cpL[i], i=1..np)], sol);      
    end if;
end proc;


FindLinearRelation2 := proc(pLL, vars, y)
    
    local sol, nvars, npLL, i, j, tt, k, pvars, idts, pLLs, ps, poly, su, sol2;


   ## pLL is a list of lists of polynomials

      npLL := nops(pLL);
      if npLL = 0 then
         return(vars);
      else
         tt := time():
         sol := FindLinearRelationOneList(pLL[1], vars, y);
         print(time_1);
         tt := time()-tt; print(tt);
        
         ### find the free variables
         nvars := [];
         pvars := [];
         idts := indets(sol); 
         for i from 1 to nops(vars) do
             if has(idts, vars[i]) then
                nvars := [op(nvars), vars[i]];
                pvars := [op(pvars), i];
             end if;
         end do;
    end if;


         ### update the other lists in pLL
         pLLs := [];
         for i from 2 to npLL do
             ps := [];
             for j from 1 to nops(pvars) do
                 poly := 0;
                 for k from 1 to nops(pLL[i]) do
                     poly := normal(poly + coeff(sol[k], nvars[j], 1)*pLL[i][k]);
                 end do;
                 ps := [op(ps), poly];
             end do;  
             pLLs := [op(pLLs), ps];
         end do;
         
tt := time():
        su := FindLinearRelation2(pLLs, nvars, y);
print(time_2);
tt := time()-tt;
print(tt);

        if nops(su)=0 then
           return([]);
        else
        sol2 := {seq(nvars[i]=su[i], i=1..nops(nvars))};
        eval(sol, sol2);
        end if;

end proc;         
 

###############################################################
  
  IsZeroVector :=proc(v)
    local i, vs;
    vs := convert(v, list);
    for i from 1 to nops(vs) do
        if not Testzero(vs[i]) then
           return(false)
        end if;
    end do;
    true;
  end proc;

IsZeroList :=proc(v)
    local i;
    for i from 1 to nops(v) do
        if not Testzero(v[i]) then
           return(false)
        end if;
    end do;
    true;
  end proc;

  ####################################################################
  # This code is from OreTools_ModularGCRDLCLM                          
  # Name: LowerBoundForRank                                              
  # Calling sequence:                                                   
  #                   LowerBoundForRank(M, p)                           
  #                                                                     
  # Input: M, A matrix with entries being multivariate polynomials      
  #           over p:                                                   
  #        p, a prime:                                                  
  # Output: b, a nonnegative integer less than or                       
  #            equal to the rank of M.                                         
  ####################################################################

   LowerBoundForRank := proc(M, p)
        local rn, cn, vars, g, e, N, i, j, en, s, b, de, dei, MM:

        #------------------------------
        # 1. Collect info about matrix
        #------------------------------

        rn := op(1, M)[1]:
        cn := op(1, M)[2]:
        vars := [op(indets(M))]:


        #---------------------------------------
        # Calculate rowwise common denominators
        #---------------------------------------

        de := []:
        for i from 1 to rn do
            dei := []:
            for j from 1 to cn do
                dei := [op(dei), denom(M[i,j])]:
            end do:
            de := [op(de), Lcm(op(dei)) mod p]:
        end do:

        MM := Matrix(rn, cn):
        for i from 1 to rn do
            for j from 1 to cn do
                MM[i,j] := Normal(de[i]*M[i,j]) mod p:
            end do:
        end do:


        #----------------------------------
        # 2. Choose two evaluation points
        #----------------------------------

        g := rand(1..p-1):
        e := [seq(g(), i=1..nops(vars))]:
        s := {seq(vars[i]=e[i], i=1..nops(vars))}:

        #---------------
        # 3. Evaluation
        #---------------

        N := Matrix(rn, cn):
        for i from 1 to rn do
                for j from 1 to cn do
                    N[i, j] := Expand(subs(s, MM[i, j])) mod p:
                end do:
        end do:

        #--------------------
        # 4. Compute rank(N)
        #--------------------

        LinearAlgebra:-Modular:-RowReduce(p, N, rn, cn, cn, 0, 
                      0, b, 0, 0, false):
        return(b):

   end proc:

end module:




   
   
   
   
   
  


## Copyright (c) 2010, INRIA & KLMM. All rights reserved.
## Authors: Shaoshi Chen <schen@amss.ac.cn> and Ziming Li <zmli@mmrc.iss.ac.cn>.

List := module()
   description "handling lists of partial fractions";
   export 
          AddSimpleFraction,
          AddTwoPartialFractions,
          DiffPartialFraction,
          FromListToNormal,
          FromListToPartialFraction,
          ScalarMultiply;
   local 
          AddBranch,
          AddBranch1,
          DiffBranch,
          DiffSimpleFraction,
          NormalizeBranch,
          SumOfSimpleFractions;

   option package;
 



#############################################################
# [y, p, fractionlist]
# fractionlist = [[d_1, [c_11, n_11], ..., [c_1n1, n_1k1]],
#                            .....
#                  [d_m, [c_m1, n_m1], ..., [c_m1, n_mkm]]]
##############################################################

##############################################################
# Name: FromListToPartialFraction
# Calling Sequence
#      
#       FromListToPartialFraction(L)
#
# Input: L, a list for a partial fraction decomposition
# Output: the partial fraction representation of the rational
#         function represented by L
##############################################################

FromListToPartialFraction := proc(L)
    local p, F, R, i, r, M, d, j;
    p := L[2];
    F := L[3];
    if F = [] then
       return p;
    end if;
    R := p;
    for i from 1 to nops(F) do
        r := 0;
        M := F[i];
        d := M[1];
        for j from 2 to nops(M) do
            r := r + M[j][1]*M[j][2]/d^(j-1);
        end do;
        R := R + r;
     end do;
     R;
end proc;

#############################################################
# Name: FromListToNormal
# Calling sequence: 
#                  FromListToNormal(L)
# Input: L, a list for a partial fraction decomposition
# Output: the usual normalized rational form of the function 
#         represented by L
##############################################################
FromListToNormal := proc(L)
    local p, F, R, i, r, M, d, j;
    p := L[2];
    F := L[3];
    if F = [] then
       return p;
    end if;
    R := p;
    for i from 1 to nops(F) do
        r := 0;
        M := F[i];
        d := M[1];
        for j from 2 to nops(M) do
            r := normal(r + normal(M[j][1]*M[j][2]/d^(j-1)));
        end do;
        R := R + r;
     end do;
     R;
end proc;



##############################################################
# Name: ScalarMultiply
# Calling Sequence
#      
#       ScalarMultiply(c, L)
#
# Input:  c, a rational function free of the main variable in L
#         L, a list for a partial fraction decomposition
# Output: c*L, in list representation
##############################################################

ScalarMultiply := proc(c, L)
    local var, p, F, cp, de, nu, G, i, M, T, j, cs, ns, k;
    (var, p, F) := L[1], L[2], L[3]; 
    cp := normal(c);
    (de, nu) := denom(c), numer(c);
    if degree(de, var) > 0 or degree(nu, var) > 0 then
       error "the first input is not a scalar with respect to the second";
    end if;
    if nu = 0 then
       return [var, 0, []];
    end if;
    G := []; 
    for i from 1 to nops(F) do
        M := F[i];
        T := table([]); 
        T[1] := M[1];
        for j from 2 to nops(M) do
            (cs, ns) := M[j][1], M[j][2];
            if ns  = 0 then
               T[j] := [1, 0];
            else
               T[j] := [normal(cp*cs), ns];
            end if;
        end do;
        G := [op(G), [seq(T[k], k=1..nops(M))]];
     end do;
     [var, cp*p, G];
end proc;



###############################################################
# Name: AddSimpleFraction
# Calling Sequence
#      
#       AddSimpleFraction(f, L)
#
# Input:  f, a simple fraction given by [c, de, nu, m]
#            representing c*nu/de^m;
#         where nu is nonzero, degree(nu, v) < degree(de, v)
#         (v the main variable of L), de is a dinominator in L
#         L, a list for a partial fraction decomposition
# Output: f + L, in list representation
###############################################################

AddSimpleFraction := proc(f, L)
    local c, de, nu, n, var, p, F, m, i, M, ds, r, dr, nr, l, G, np,
          N, newf, cp, Fp, k, R; 
    (c, de, nu, n) := f[1], f[2], f[3], f[4];
    (var, p, F) := L[1], L[2], L[3];
    if nu = 0 then
       return L;
    end if;
    Fp := 0;
    m := nops(F);  
    for i from 1 to m do
        M := F[i];
        ds := M[1];
        r := normal(ds/de);
        (dr, nr) := denom(r), numer(r);
        #print(degrees); print(degree(dr, var)); print(degree(nr, var));
        if degree(dr, var) = 0 and degree(nr, var) = 0 then
           #print(find_the_denom);
           c := normal(c/r);
           l := nops(M);
 
           # update a branch

           if l - 1 < n then
             #  G := [op(M), [1, 0]$(n-l-1), [c, nu]]; 
             G := [op(M), [1, 0]$(n-l), [c, nu]]; 
           else
              N := M[n+1];
              np := normal(N[1]*N[2] + c*nu);

              # update a simple fraction

              if np = 0 then
                 newf := [1, 0];
              else
                 cp := content(np, var, 'np');
                 newf := [cp, np];
              end if;

              # update a branch

              if np = 0 and l = n + 1 then
                 print(here);
                 G := [op(1..l-1, M)];
                 k := nops(G);
                 while k > 1 and G[k][2] = 0 do
                    G := [op(1..k-1, G)];
                    k := k - 1;
                 end do;
              else
                 G := [op(1..n, M), newf, op(n+2..l, M)];
              end if;
           end if;

           # update the list

          if nops(G) < 2 then
             Fp := [op(1..i-1, F), op(i+1..m, F)];
          else
             Fp := [op(1..i-1, F), G, op(i+1..m, F)];
          end if; 
          break;
       end if;
    end do;
    if Fp = 0 then
       R := [var, p, [op(F), [de, [1,0]$(n-1), [c, nu]]]];
    else
       R := [var, p, Fp];
    end if;
    R;
end proc;
                 
#########################################################
# Name: DiffSimpleFraction
# Calling sequence
#       DiffSimpleFraction(F, x, var)
# Input: F, a nonzero simple fraction in the list representation
#        x, a variable;
#        var, the main variable;
# Output: The result of diff(F, x)
#########################################################

DiffSimpleFraction := proc(F, x, var)
    local c, de, nu, m, np, cp, cs, ns;
    (c, de, nu, m) := F[1], F[2], F[3], F[4]; 
    if nu = 0 then
       return [[1,0], [1,0]];
    end if; 
    np := rem(-nu*diff(de, x), de, var, 'q');
    np := primpart(np, var, 'cp');
    cp := c*m*cp;
    ns := primpart(diff(c, x)*nu + c*diff(nu, x) + c*m*q, var, 'cs');
    if ns = 0 then
       cs := 1;
    end if;
    if np = 0 then
       cp := 1;
    end if;
    [[cs, ns], [cp, np]];
end proc;

#########################################################
# Name: SumOfSimpleFractions
# Calling sequence
#       SumOfSimpleFractions(F1, F2, r)
# Input: F1, F2, two simple fractions with the same denominators and 
#        multiplicities (in the list representation)
#        var, the main variable;
# Output: The result F1 + F2; 
#########################################################

SumOfSimpleFractions := proc(F1, F2, var)
    local c1, de, nu1, m, c2, nu2, nu, c;
    (c1, nu1) := F1[1], F1[2];
    (c2, nu2) := F2[1], F2[2]; 
    nu := primpart(c1*nu1+c2*nu2, var, 'c');
    if c = 0 then
       c := 1;
    end if;
    [c, nu];
end proc;
    
#########################################################
# Name: DiffBranch
# Calling sequence
#       DiffBranch(B, x, var)
# Input: B, a nonempty branch of a partial fraction 
#        x, a variable;
#        var, the main variable;
# Output: The result diff(B, x)
#########################################################

DiffBranch := proc(B, x, var)
    local de, M, m, T, i, l, F, u, S;
    m := nops(B);
    de := B[1];
    M := [op(2..m, B)];
    m := m - 1;
    T := table([]);
    for i from 1 to m do
        T[i] := DiffSimpleFraction([M[i][1], de, M[i][2], i], x, var);
    end do; 

    F := table([]);
    F[1] := T[1][1];
    u := T[1][2]; 
    for i from 2 to m do
        F[i] := SumOfSimpleFractions(u, T[i][1], var);
        u := T[i][2];  
    end do;
    F[m+1] := [T[m][2][1], T[m][2][2]];

    l := m + 1;
    for i from m+1 by -1 to 1 do
        if F[i][2] = 0 then
           l := l - 1;
        else
           break;
        end if;
    end do;
    [seq(F[i], i=1..l)];
   # [de, seq(F[i], i=1..l)];
end proc; 
             
######################################################################
# Name: DiffPartialFraction
# Calling sequence
#       DiffPartialFraction(F, x)
# Input: F, a partial fraction in list representation;
#        x, a name;
# Output: diff(F, x), in list representation
######################################################################

DiffPartialFraction := proc(F, x)
   local var, p, pp, G, l, H, i, de, r;
   (var, p) := F[1], F[2];
   pp := diff(p, x);
   G := F[3];
   if G = [] then
      return [var, pp, []];
   end if;
   l := nops(G);
   H := table([]);
   for i from 1 to l do
       de := G[i][1];
       r := DiffBranch(G[i], x, var);
       if r = [] then
          H[i] := op([]);
       else
          H[i] := [de, op(r)];
       end if;
   end do; 
   [var, pp, [seq(H[i], i=1..l)]];
end proc;

#test
NormalizeBranch := proc(B)
  local nu, l, de, i;
  nu := 0;
  l := nops(B);
  de := B[1];
  for i from 2 to l do
      nu := normal(nu + B[i][1]*B[i][2]*de^(l-i));
  end do;
  normal(nu/de^(l-1)); 
end proc;

########################################################################
# Name: AddBranch1
# Calling sequence:
#       AddBranch1(B, var)
# Input: B, a list of branches whose denominators are equal
#        var, the main variable
# Output: the sum of these branches
#########################################################################

AddBranch1 := proc(B, var)
   local b, de, nu1, C, nu2, l1, l2, l, T, nu, c, ll, d, i, m;
   b := nops(B);
   if b = 0 then
      error "the input is an empty list";
   end if;
   if b = 1 then
      return B[1];
   end if;
   m := nops(B[1]); 
   de := B[1][1];
   nu1 := [op(2..m, B[1])];
   C := AddBranch1([op(2..b, B)], y);
   if nu1 = [] then
      return C;
   end if;
   m := nops(B[2]);
   nu2 := [op(2..m, C)];
   if nu2 = [] then
      return B[1]; 
   end if;
   l1 := nops(nu1);
   l2 := nops(nu2); 
   l := min(l1, l2);
   ll := max(l1, l2);
   T := table([]);
   for i from 1 to l do
       nu := normal(nu1[i][1]*nu1[i][2]+nu2[i][1]*nu2[i][2]);
       if nu = 0 then
          T[i] := [1, 0];
       else
          nu := primpart(nu, var, 'c');
          T[i] := [c, nu];
       end if;
   end do;
   if l1 > l then
      for i from l+1 to ll do
          T[i] := nu1[i];
      end do;
   end if; 
   if l2 > l then
      for i from l+1 to ll do
          T[i] := nu2[i];
      end do;
   end if; 
   d := ll;
   for i from ll by -1 to 1 do
       if T[i][2] = 0 then
          d := d - 1;
       else
          break;
       end if;
   end do;
   [de, seq(T[i], i=1..d)];
end proc;  

###########################################################
#Name: AddBranch
#Calling sequence:
#      AddBranch(B, L)
#Input: B, a branch
#       L, a list of partial fractions
#Output: B + L in list representation
############################################################
  

AddBranch := proc(B, L)
   local var, p, de, F, l, i, H, dh, G;
   if nops(B) = 1 then
      return L;
   end if;
   (var, p) := L[1], L[2];
   F := L[3];
   if F = [] then
      return [var, p, [B]];
   end if;
   de := B[1];
   l := nops(F);
   for i from 1 to l do
       H := F[i];
       dh := H[1];
       if Testzero(de-dh) then
          G := AddBranch1([B, H], var);
          if nops(G) = 1 then
             return [var, p, [op(1..i-1, F), op(i+1..l, F)]];
          end if;
          return [var, p, [op(1..i-1, F), G, op(i+1..l, F)]];
        end if;
    end do;
    [var, p, [op(F), B]];
end proc; 

###########################################################
#Name: AddTwoPartialFractions
#Calling sequence:
#      AddTwoPartialFractions(B, L)
#Input: P, Q Two partial fractions
#       
#Output: P + Q in list representation
############################################################
  
AddTwoPartialFractions := proc(P, Q)
   local var, p, i, R;
   var := P[1];
   p := normal(P[2] + Q[2]);
   R := Q;
   for i from 1 to nops(P[3]) do
       R := AddBranch(P[3][i], R);
   end do;
   [var, p, R[3]];
end proc;
       
      
end module: