## Copyright (c) May 2024, FZU. All rights reserved.
## Author: Hui Huang <huanghui@fzu.edu.cn>
## Load the modules: QHypergeometricCT.mm, RationalTools.mm
################################################################################

QHypergeometricCT := module()
  description "Creative telescoping for q-hypergeometric terms via reduction":

  export
        QHypergeomTelescoping,
	QHypergeomReduction,
	QShellReduction,
	QPolynomialReduction,
	QOrderBounds,
	QShiftQuotient,
	QRNF,
	QSplitFactor,
	QDispersionSet,
	QShiftFactorization,
	VerifyQHypergeomTelescoping,
	VerifyQHypergeomReduction,
	VerifyQShellReduction;
  local
	QSpecialReduction,
	QNormalReduction,
	LocalQNormalReduction,
	QSimpleReduction,
	LinearizeQRemainders,
	FindLinearRelation,
	LowerBoundForRank,
	QUpperBound,
	QLowerBound,
	MinPosIntSet,
	SubsPower,
	LocalSubsPower,
	QNormalRNF,
	IsIrredQIntegerLinear,
	IsQShiftEquivalent,
	IsQPower;
  option package:
  
##########################  Task  Description ##################################
# Creative telescoping problem for q-hypergeometric terms:
# Given: a q-hypergeometric term T over QQ(q)(q^n,q^k) with
# 	 q an indeterminate over QQ;
# Find:  a nonzero linear q-difference operator L in QQ(q,q^n)[Qn]
#        such that
#                   L(T) = (Qk-1)(G),
#        where G is a q-hypergeometric term over QQ(q)(q^n,q^k).
#
# Notation: F = QQ(q,q^n), Qn(q^n) = q^(n+1), Qk(q^k) = q^(k+1),
# 	    Sx(x) = q*x, Sy(y) = q*y.
################################################################################

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

################################################################################
# Name: QHypergeomTelescoping
# Calling sequence:
#                  QHypergeomTelescoping(T,n,k,q,Qn,'cert')
# Input: T,     a q-hypergeometric term over QQ(q)(q^n,q^k);
#        n,k,   two variable names;
#        q,     an indeterminate over QQ used for the parameter q;
#        Qn,    an operator name used for the q-shift operator w.r.t. n;
#        cert,  (optional) an unevaluated name used for storing the certificate.
# Output: L,  a nonzero linear q-difference operator in QQ(q,q^n)[Qn]
#             such that
#                      L(T) = (Qk-1)(G),
#             where G is a q-hypergeometric term over QQ(q)(q^n,q^k), 
#	      if such an L exists.
#	      The optional 6th argument cert is assigned [g,H] with G = g*H,
#	      where g is a rational function in QQ(q)(q^n,q^k) and H is a 
#	      q-hypergeometric term over QQ(q)(q^n,q^k) whose Qk-quotient
#	      is Qk-standard.
#         "No telescoper exists!", 
#             if there is no such L.
################################################################################

QHypergeomTelescoping := proc(T,n,k,q,Qn)
  local fk,u,v,s,t,S,H,g,h,p,m,fn,d,tp,i,N,gL,hL,pL,gg,hh,pp,cL,z;

  #------------------------------------
  # compute a standard RNF of Qk(T)/T
  #------------------------------------
  fk := QShiftQuotient(T,k,q,y);
  fk := normal(SubsPower(q^n=x,fk));
  u,v,s,t := op(QRNF(fk,q,y,'standard'));

  #-------------------------------------------------
  # compute the multiplicative decomposition of T
  #-------------------------------------------------
  S := subs(x=q^n,y=q^k,s/t);
  H := normal(T/S);

  #------------------------------------------------------
  # apply the shell reduction to s/t w.r.t. u/v
  #------------------------------------------------------
  g,h,p := op(QShellReduction(s,t,u,v,q,y,'special-free'));

  #------------------------------------------
  # trivial case: T is q-summable w.r.t. k
  #------------------------------------------
  m := nops(h);
  if m = 0 and p = 0 then
    if nargs = 6 then
      assign(args[6]=[subs(x=q^n,y=q^k,RationalTools:-NormalizeRatSum(g)),H]);
    fi;
    return(1);
  fi;

  #-----------------------------------------------------------
  # trivial case: T is free of n and non-q-summable w.r.t. k
  #-----------------------------------------------------------
  fn := QShiftQuotient(T,n,q,x);
  if fn = 1 then
    if nargs = 6 then assign(args[6]=[0,H]); fi;
    return(-1+Qn);
  fi;
  
  #----------------------------------------------------
  # check if the denominator of h is q-integer linear
  #----------------------------------------------------
  d := map(z->z[2],h);
  tp := table([]);
  for i from 1 to m do
    if not IsIrredQIntegerLinear(d[i],x,y,'tp[i]') then
      return("No telescoper exists!");
    fi;
  od;
  tp := convert(tp,list);

  #--------------------------------------------
  # compute the q-shift quotient of H w.r.t. n
  #--------------------------------------------
  fn := normal(SubsPower(q^k=y,fn));
  N := normal(fn*(s/t)/subs(x=q*x,s/t));

  #----------------------------------------------------------
  # iteratively compute a minimal telescoper for T w.r.t. n
  #----------------------------------------------------------
  if nargs = 6 then gL := [RationalTools:-NormalizeRatSum(g)]; fi;
  hL := [h];
  pL := [p];
  i := 1;
  while i > 0 do
    gg,hh,pp := op(QSimpleReduction(hL[i],pL[i],u,v,N,q,x,y));
    g,hh,p := op(LinearizeQRemainders(hh,d,tp,u,v,q,y));
    if nargs = 6 then
      gL := [op(gL),RationalTools:-NormalizeRatSum(subs(x=q*x,gL[i])*N+gg+g)];
    fi;
    hL := [op(hL),hh];
    pL := [op(pL),RationalTools:-NormalizeRatSum(p+pp)];
    cL := FindLinearRelation(hL,pL,y);
    m := nops(cL);
    if m = 0 then
      i := i+1;
    else
      if nargs = 6 then
	 g := RationalTools:-NormalizeRatSum(add(cL[j]*gL[j],j=1..m));
         assign(args[6]=[subs(x=q^n,y=q^k,g),H]);
      fi;
      return(add(subs(x=q^n,cL[j])*Qn^(j-1),j=1..m));
    fi;
  od;
end:

################################################################################
# Name: QHypergeomReduction
# Calling sequence:
#                  QHypergeomReduction(T,k,q)
# Input: T,  a q-hypergeometric term over F(q^k);
#        k,  a variable name;
#	 q,  an indeterminate over QQ used for the parameter q.
# Output: [[g,r],H],  where
# 	    g,r - two rational functions in F(q^k),
#	    H   - a q-hypergeometric term over F(q^k) with Qk(H)/H = K
#	     	  being Qk-standard,
#           such that
#           (1) T = (Qk-1)(g*H)+r*H;
#           (2) r is a Qk-remainder w.r.t. K.
################################################################################

QHypergeomReduction := proc(T,k,q)
  local f,u,v,s,t,S,H,g,h,p,r,y,z;
  
  #------------------------------------
  # compute a standard RNF of Qk(T)/T
  #------------------------------------
  f := QShiftQuotient(T,k,q,y);
  u,v,s,t := op(QRNF(f,q,y,'standard'));

  #-------------------------------------
  # compute the q-hypergeometric term H
  #-------------------------------------
  S := subs(y=q^k,s/t);
  H := normal(T/S);

  #------------------------------------------------------
  # apply the shell reduction to s/t w.r.t. u/v
  #------------------------------------------------------
  g,h,p := op(QShellReduction(s,t,u,v,q,y,'special-free'));
  r := RationalTools:-NormalizeRatSum(add(map(z->z[1]/z[2]^z[3],h))+p/v);
  
  [map2(subs,y=q^k,[RationalTools:-NormalizeRatSum(g),r]),H];   
end:

################################################################################
# Name: QShellReduction
# Calling sequence:
#                  QShellReduction(s,t,u,v,q,y,'special-free')
# Input: s,t,u,v,         four polynomials in F[y] with
#	              	  (1) u/v being Sy-standard;
#                     	  (2) gcd(s,t) = 1;
#		      	  (3) t(0) <> 0 if the optional 7th argument
#		      	      'special-free' is given;
#	 q,           	  an indeterminate over QQ used for the parameter q;
#	 y,          	  a variable name;
#	 'special-free',  (optional) an argument name used for indicating
#	 	      	  the input.
# Output: [g,h,p],  where
# 	    g - a rational function in F(y),
#	    h - a list of the form [[a_1,d_1,k_1],...,[a_m,d_m,k_m]] for
#	      	nonzero polynomials a_i,d_i in F[y] and positive integers k_i,
#	    p - a polynomial in F[y],
#	    such that
#           (1) s/t = ((u/v)*Sy-1)(g) + sum(a_i/d_i^(k_i),i=1..m) + p/v;
#           (2) degree(a_i,y) < degree(d_i^(k_i),y) and gcd(a_i,d_i) = 1;
#	    (3) d_i is q-primitive, irreducible and strongly coprime with u/v;
#	    (4) the d_i are pairwise q-shift inequivalent;
#           (5) p is in the standard complement of the image of the map
#	    	for polynomial reduction w.r.t. u/v.
################################################################################

QShellReduction := proc(s,t,u,v,q,y)
  local a,p,an,dn,gs,bs,c,ds,as,g,h,b,f;

  #------------------------------------------------
  # compute the canonical representation of s/t
  #------------------------------------------------
  a := rem(s,t,y,'p');
  if nargs = 7 and type(args[7],'identical(special-free)') then
    an,dn := a,t;
    gs,bs := 0,0;
  else
    c,ds,dn := op(QSplitFactor(t,q,y));
    as,an := op(RationalTools:-ParFracDecomp(a/c,[ds,dn],y));
    #------------------------------------------------
    # apply the special reduction to the special 
    # part of s/t w.r.t. u/v
    #------------------------------------------------
    gs,bs := op(QSpecialReduction(as,ds,u,v,q,y));
  fi;
  
  #------------------------------------------------
  # apply the normal reduction to the normal 
  # part of s/t w.r.t. u/v
  #------------------------------------------------
  g,h,b := op(QNormalReduction(an,dn,u,v,q,y));

  #--------------------------------------------------------
  # apply the polynomial reduction to p*v+bs+b w.r.t. u/v
  #--------------------------------------------------------
  f,p := op(QPolynomialReduction(RationalTools:-NormalizeRatSum(p*v+bs+b),u,v,q,y));
  g := gs+g+f;
  
  [g,h,p];
end:

################################################################################
# Name: QPolynomialReduction
# Calling sequence:
#                  QPolynomialReduction(b,u,v,q,y)
# Input: b,u,v,  three polynomials in F[y] with u/v being Sy-standard;
#	 q,      an indeterminate over QQ used for the parameter q;
#	 y,      a variable name.
# Output: [a,p],  where
# 	    a,p - two polynomials in F[y]
#	    such that
#           (1) b = u*Sy(a) - v*a + p;
#           (2) p is in the standard complement of the image of the map
#	    	for polynomial reduction w.r.t. u/v.
################################################################################

QPolynomialReduction := proc(b,u,v,q,y)
  local du,dv,lu,lv,m,d,a,p,i,g,pg,lp,bm,pm,r,gr,dr,pr,lr,z;
  
  #-----------------------
  # trivial case: b = 0
  #-----------------------
  if b = 0 then return([0,0]); fi;
  
  #--------------------------------------------------------
  # construct an echelon basis for the image of the map
  # for polynomial reduction w.r.t. u/v
  #--------------------------------------------------------
  du := degree(u,y);
  dv := degree(v,y);
  lu := lcoeff(u,y);
  lv := lcoeff(v,y);
    
  if du <> dv or not IsQPower(lv/lu,q,'m') or not type(m,nonnegative) then
    #-------------------------------------------------
    # case I: either degree(u,y) <> degree(v,y) 
    # 	    or lcoeff(u,y)*q^i - lcoeff(v,y) <> 0
    #	    for all nonnegative integers i
    #-------------------------------------------------
    d := max(du,dv);
    lu := coeff(u,y,d);
    lv := coeff(v,y,d);
    a := 0;
    p := b;
    i := degree(p,y)-d;    
    while i >= 0 do
      g := y^i/(lu*q^i-lv);
      pg := RationalTools:-NormalizeRatSum(u*q^i*g-v*g);
      lp := lcoeff(p,y);
      a := a+lp*g;
      p := RationalTools:-NormalizeRatSum(p-lp*pg);
      i := degree(p,y)-d;
    od;
  elif du = 0 then
    #-------------------------------------------------------
    # case II: u/v = q^(-m) for some nonnegative integer m
    #-------------------------------------------------------
    p := coeff(b,y,m)*y^m;
    a := collect(b-p,y);
    a := add(map(z->normal(z/((q^degree(z,y)-q^m)*u)),`if`(type(a,`+`),[ob(a)],[a])));
  else
    #-------------------------------------------------
    # case III: degree(u,y) = degree(v,y) > 0 and
    # 	    lcoeff(u,y)*q^m - lcoeff(v,y) = 0
    #	    for some nonnegative integer m
    #-------------------------------------------------
    a := 0;
    bm := coeff(b,y,du+m)*y^(du+m);
    p := RationalTools:-NormalizeRatSum(b-bm);
    i := degree(p,y)-du;
    while i >= 0 do
      g := y^i/(lu*q^i-lv);
      pg := RationalTools:-NormalizeRatSum(u*q^i*g-v*g);
      lp := lcoeff(p,y);
      a := a+lp*g;
      p := RationalTools:-NormalizeRatSum(p-lp*pg);
      pm := coeff(p,y,du+m)*y^(du+m);
      bm := bm+pm;
      p := RationalTools:-NormalizeRatSum(p-pm);
      i := degree(p,y)-du;
    od;
    
    r := RationalTools:-NormalizeRatSum(u*q^m*y^m-v*y^m);
    gr := y^m;
    i := degree(r,y)-du;
    while i >= 0 do
      g := lcoeff(r,y)*y^i/(lu*q^i-lv);
      gr := gr-g;
      r := RationalTools:-NormalizeRatSum(r-u*q^i*g+v*g);
      i := degree(r,y)-du;
    od;

    dr := i+du;
    lr := lcoeff(r,y);

    pr := coeff(p,y,dr);
    a := a+pr*gr/lr;
    p := RationalTools:-NormalizeRatSum(p+bm-pr*r/lr);
  fi;
  
  [a,p];  
end:

################################################################################
# Name: QOrderBounds
# Calling sequence:
#                  QOrderBounds(T,n,k,q,A)
# Input: T,    a q-hypergeometric term over QQ(q)(q^n,q^k);
#        n,k,  two variable names;
#        q,    an indeterminate over QQ used for the parameter q;
#        A,    (optional) an argument name used for specifying the output.
# Output: [up,lo],  if T has a telescoper and the optional 5th argument is
# 	  	    missing, where
# 	   up,lo - two nonnegative integers,
#	   such that up and lo are respectively upper and lower bounds
# 	   for the order of minimal telescopers for T.
#	  up,  	    if T has a telescoper and the optional 5th argument
#	  	    'upper' is given.
#         lo,	    if T has a telescoper and the optional 5th argument
#	  	    'lower' is given.
#         "No telescoper exists!",
#	            if T has no telescopers.
################################################################################

QOrderBounds := proc(T,n,k,q)
  local fk,y,x,u,v,s,t,h,p,m,tab,i,d,kk,tp,po,l,PL,up;

  #------------------------------------
  # compute a standard RNF of Qk(T)/T
  #------------------------------------
  fk := QShiftQuotient(T,k,q,y);
  fk := normal(SubsPower(q^n=x,fk));
  u,v,s,t := op(QRNF(fk,q,y,'standard'));
  
  #------------------------------------------------------
  # apply the shell reduction to s/t w.r.t. u/v
  #------------------------------------------------------
  h,p := op(QShellReduction(s,t,u,v,q,y,'special-free')[2..-1]);

  #------------------------------------------
  # trivial case: h is an empty list
  #------------------------------------------
  m := nops(h);
  if m = 0 then
    if nargs = 5 and type(args[5],'identical(upper)') then
      return `if`(p=0,0,max(degree(u,y),degree(v,y),1));
    elif nargs = 5 and type(args[5],'identical(lower)') then
      return `if`(p=0,0,1);
    else
      return `if`(p=0,[0,0],[max(degree(u,y),degree(v,y),1),1]);
    fi;
  fi;

  #----------------------------------------------------
  # check if the denominator of h is q-integer linear
  #----------------------------------------------------
  tab = table([]);
  for i from 1 to m do
    d,kk := op(2..-1,h[i]);
    if not IsIrredQIntegerLinear(d,x,y,'tp','po') then
      return("No telescoper exists!");
    fi;
    if assigned(tab[tp],'PL') then
      i := 1;
      while i <= nops(PL) and not IsQShiftEquivalent(po,PL[i,1],q,y,'l') do
        i := i+1;
      od;
      if i <= nops(PL) then
        PL[i,2] := [op(PL[i,2]),[kk,l]];
      else
        tab[tp] := [op(PL),[po,[[kk,0]]]];
      fi;
    else
      tab[tp] := [[po,[[kk,0]]]];
    fi;
  od;
  
  #----------------------------------
  # compute the order bounds
  #----------------------------------
  if nargs = 5 and type(args[5],'identical(upper)') then
    up := max(degree(u,y),degree(v,y));
    if up = 0 and IsQPower(u/v,q) then up := 1; fi;
    return(up+QUpperBound(tab,y));
  elif nargs = 5 and type(args[5],'identical(lower)') then
    return(QLowerBound(tab));
  fi;

  up := max(degree(u,y),degree(v,y));
  if up = 0 and IsQPower(u/v,q) then up := 1; fi;

  [up+QUpperBound(tab,y),QLowerBound(tab)];
end:

################################################################################
# Name: QShiftQuotient
# Calling sequence:
#                  QShiftQuotient(T,k,q,y)
# Input: T,  a term over F(q^k);
#        k,  a variable name;
#	 q,  an indeterminate over QQ used for the parameter q;
#	 y,  a variable name used for replacing q^k.
# Output: f,  a rational function in F(y) such that f(q^k) = Qk(T)/T,
#	      if T is a q-hypergeometric term over F(q^k);
#	  "Input is not a q-hypergeometric term in q^k",
#	      otherwise.
################################################################################

QShiftQuotient := proc(T,k,q,y)
  local f;

  if not QDifferenceEquations:-IsQHypergeometricTerm(T,k,q^k=y,'f') then
   error "Input is not a q-hypergeometric term in %1",q^k;
  fi:
  
  f;
end:

################################################################################
# Name: QRNF
# Calling sequence:
#                  QRNF(f,q,y,'standard')
# Input: f,           a nonzero rational function in F(y);
# 	 q,           an indeterminate over QQ used for the parameter q;
#        y,           a variable name;
#	 'standard',  (optional) an argument name used for specifying
#	 	      the output.
# Output: [u,v,s,t],  where
#           u,v,s,t - four polynomials in F[y],
#	    such that
#           (1) f = K*Sy(S)/S with K = u/v and S = s/t;
#	    (2) K is Sy-reduced;
#           (3) gcd(s,t) = gcd(u,s*Sy(t)) = gcd(v,Sy(s)*t) = 1;
#	    (4) v,s,t are primitive w.r.t. y and t(0) <> 0;
#	    (5) K is Sy-standard if the optional 4th argument
#	    	'standard' is given.
################################################################################

QRNF := proc(f,q,y)
  local c,fs,fn,u,v,s,t,k,uv0,m;

  #--------------------------------------------
  # compute the splitting factorization of f
  #--------------------------------------------
  c,fs,fn := op(QSplitFactor(f,q,y));
  
  #---------------------------------------
  # compute an RNF of the normal part fn
  #---------------------------------------
  u,v,s,t := op(QNormalRNF(fn,q,y));

  #--------------------------
  # standardize the RNF
  #--------------------------
  k := degree(fs,y);
  uv0 := normal(coeff(u,y,0)/coeff(v,y,0));
  if k = 0 and nargs = 4 and type(args[4],'identical(standard)') and\
     IsQPower(c*uv0,q,'m') and m > 0 then
    u := normal(u/uv0);
    s := y^m*s;
  elif k >= 0 then
    u := y^k*normal(c*u);
  else
    u := normal(c*u);
    v := y^(-k)*v;
  fi;
  
  [u,v,s,t];
end:

################################################################################
# Name: QSplitFactor
# Calling sequence:
#                  QSplitFactor(f,q,y)
# Input: f,  a rational function in F(y);
# 	 q,  an indeterminate over QQ used for the parameter q;
#        y,  a variable name.
# Output: [c,fs,fn],  where
#           c  - an element in F,
#	    fs - a power of y,
#	    fn - a q-primitive rational function in F(y),
#	    such that
#	             f = c*fs*fn.
################################################################################

QSplitFactor := proc(f,q,y)
  local ff,c,nf,df,nk,dk,fs,fn;

  #-------------------------------
  # trivial case: f is free of y
  #-------------------------------
  ff := normal(f);
  if type(ff,freeof(y)) then
    return([ff,1,1]);
  fi;

  #-----------------------------
  # compute the special part
  #------------------------------
  c := normal(content(numer(ff),y,'nf')/content(denom(ff),y,'df'));
  nk := ldegree(nf,y);
  dk := ldegree(df,y);
  if nk > 0 then
    fs := y^nk;
    nf := normal(nf*y^(-nk));
  elif dk > 0 then
    fs := y^(-dk);
    df := normal(df*y^(-dk));
  else
    fs := 1;
  fi;

  #-----------------------------
  # compute the normal part
  #-----------------------------
  fn := normal(nf/df);

  [c,fs,fn];
end:

################################################################################
# Name: QDispersionSet
# Calling sequence:
#                  QDispersionSet(a,b,q,y)
# Input: a,b,  two q-primitive polynomials in F[y];
# 	 q,    an indeterminate over QQ used for the parameter q;
#        y,    a variable name.
# Output: DS,  a finite set of all integers i such that gcd(a,Sy^i(b)) <> 1.
################################################################################

QDispersionSet := proc(a,b,q,y)
  local afac,bfac,DS,i,j,m,z,zz;

  #---------------------------------------
  # trivial case: a or b is free of y
  #---------------------------------------
  if type(a,freeof(y)) or type(b,freeof(y)) then
    return({});
  fi;

  #---------------------------
  # factor the sqrfree parts
  #---------------------------
  afac := map(z->op(map(zz->zz[1],factors(z[1])[2])),sqrfree(a,y)[2]);
  bfac := map(z->op(map(zz->zz[1],factors(z[1])[2])),sqrfree(b,y)[2]);

  #----------------------------------------
  # force all factors to be q-primitive
  #----------------------------------------
  afac := map(z->numer(z),afac);
  bfac := map(z->numer(z),bfac);
  
  #---------------------------
  # find the dispersion set
  #---------------------------
  DS := {};
  for i from 1 to nops(afac) do
    for j from 1 to nops(bfac) do
      if IsQShiftEquivalent(afac[i],bfac[j],q,y,'m') then
        DS := DS union {m};
      fi;
    od;
  od;

  DS;
end:

################################################################################
# Name: QShiftFactorization
# Calling sequence:
#                  QShiftFactorization(d,q,y,[u,v]);
# Input: d,      a polynomial in F[y] with degree(d,y) > 0 and d(0) <> 0;
#	 q,      an indeterminate over QQ used for the parameter q;
#	 y,      a variable name;
#	 [u,v],  (optional) a list of polynomials in F[y] with
#	 	 u/v being Sy-reduced, used for specifying the output.
# Output: [c,[[d_1,[[k_{11},l_{11}],...,[k_{1n_1},l_{1n_1}]]],
#	      ...
#	      [d_m,[[k_{m1},l_{m1}],...,[k_{mn_m},l_{mn_m}]]]]],  where
#	    c      - a nonzero element in F,
#	    d_i    - a q-primitive irreducible polynomial in F[y],
#	    k_{ij} - a positive integer,
#	    l_{ij} - an integer,
#	    such that
#	    (1) d = c*prod(Sy^(l_{ij})(d_i^(k_{ij})),j=1..n_i,i=1..m);
#	    (2) the d_i are pairwise q-shift inequivalent;
#	    (3) l_{i1},...,l_{in_i} are distinct for each i;
#	    (4) d_i is strongly coprime with u/v if the optional 4th 
#	    	argument [u,v] is given.
################################################################################

QShiftFactorization := proc(d,q,y)
  local pd,cd,c,dfac,S,n,p,k,kl,temp,i,p1,k1,l,co,u,v,un,vn,udset,vdset,m,z;

  #--------------------------------------
  # compute the full factorization of d
  #--------------------------------------
  pd := primpart(d,y,'cd');
  c,dfac := op(factors(pd));

  #--------------------------------------------------
  # force all irreducible factors to be q-primitive
  #--------------------------------------------------
  c := normal(cd*c/mul(map(z->denom(z[1])^(z[2]),dfac)));
  dfac := map(z->[numer(z[1]),z[2]],dfac);

  #-----------------------------------------------
  # group together all q-shift equivalent factors
  #-----------------------------------------------
  S := [];
  n := nops(dfac);
  while n <> 0 do
    p,k := op(dfac[1]);
    kl := [[k,0]];
    temp := [];
    for i from 2 to n do
      p1,k1 := op(dfac[i]);
      if IsQShiftEquivalent(p1,p,q,y,'l','co') then
        kl := [op(kl),[k1,l]];
	c := normal(c*co^(k1));
      else
        temp := [op(temp),[p1,k1]];
      fi;
    od;
    S := [op(S),[p,kl]];
    dfac := temp;
    n := nops(dfac);
  od;
  
  #--------------------------------------------------------
  # make the representative polynomial in each q-shift
  # equivalence class strongly coprime with u/v
  # if the optional 4th argument [u,v] is given
  #--------------------------------------------------------
  if nargs = 4 then
    u,v := op(args[4]);
    un := QSplitFactor(u,q,y)[-1];
    vn := QSplitFactor(v,q,y)[-1];
    for i from 1 to nops(S) do
      p,kl := op(S[i]);
      l := map(z->z[2],kl);
      udset := QDispersionSet(un,p,q,y);
      if udset = {} then
        vdset := QDispersionSet(vn,p,q,y);
	m := min(min(l),min(vdset)-1);
      else
        m := max(max(l),max(udset)+1);
      fi;
      p := subs(y=q^m*y,p);
      if m < 0 then
	c := normal(c/denom(p)^(add(z[1],z=kl)));
	p := numer(p);
      fi;
      S[i] := [p,map(z->[z[1],z[2]-m],kl)];
    od;
  fi;
  
  [c,S];
end:

######## Verification

################################################################################
# Name: VerifyQHypergeomTelescoping
# Calling sequence:
#                  VerifyQHypergeomTelescoping([L,[g,H]],T,n,k,q,Qn)
# Input: [L,[g,H]],  a list with L a nonzero linear q-difference operator
# 	 	     in QQ(q,q^n)[Qn], g in QQ(q)(q^n,q^k) and H a
#		     q-hypergeometric term over QQ(q)(q^n,q^k);
#	 T,          a q-hypergeometric term over QQ(q)(q^n,q^k);
#        n,k,        two variable names;
#        q,          an indeterminate over QQ used for the parameter q;
#        Qn,         an operator name used for the q-shift operator w.r.t. n.
# Output: true,   if the following conditions are satisfied:
# 	  	  (1) L(T) = (Qk-1)(g*H);
#		  (2) the Qk-quotient of H is Qk-standard;
#         false,  otherwise.
################################################################################

VerifyQHypergeomTelescoping := proc(res,T,n,k,q,Qn);
  local L,cert,g,H,S,K,x,y,N,r,i,j,u,v,us,un,vs,vn,u0,v0,m,z;
  
  #----------------------------------------------
  # compute the kernel and the shell of Qk(T)/T
  #----------------------------------------------
  L,cert := op(res);
  g,H := op(cert);
  L := SubsPower(q^n=x,L);
  g := SubsPower({q^n=x,q^k=y},g);
  S := SubsPower({q^n=x,q^k=y},normal(T/H));
  K := SubsPower(q^n=x,QShiftQuotient(H,k,q,y));

  #----------------------------------------------
  # compute the shell of Qk(L(T))/L(T) w.r.t. K
  #----------------------------------------------
  N := SubsPower(q^k=y,QShiftQuotient(H,n,q,x));
  r := add(coeff(L,Qn,i)*subs(x=q^i*x,S)*mul(subs(x=q^j*x,N),j=0..i-1),i=0..degree(L,Qn));

  #----------------------------------------------------
  # check the first condition, which is equivalent to
  # r = (K*Sy-1)(g)
  #----------------------------------------------------
  r := RationalTools:-NormalizeRatSum(r-K*subs(y=q*y,g)+g);
  if r <> 0 then
    print(The_1st_condition_is_violated);
    return(false);
  fi;

  #-------------------------------
  # check the second condition
  #-------------------------------
  u,v := numer(K),denom(K);
  us,un := op(QSplitFactor(u,q,y)[2..-1]);
  vs,vn := op(QSplitFactor(v,q,y)[2..-1]);
  if degree(us,y)*degree(vs,y) <> 0 or nops(QDispersionSet(un,vn,q,y)) <> 0 then
    print(The_2nd_condition_1_is_violated);
    return(false);
  fi;
  
  u0,v0 := seq(subs(y=0,z),z=[u,v]);
  if u0*v0 <> 0 and IsQPower(u0/v0,q,'m') and type(m,positive) then
    print(The_2nd_condition_2_is_violated);
    return(false);
  fi;
  
  true;
end:

################################################################################
# Name: VerifyQHypergeomReduction
# Calling sequence:
#                  VerifyQHypergeomReduction([[g,r],H],T,k,q)
# Input: [[g,r],H],  a list with g,r in F(q^k) and H a q-hypergeometric 
# 	 	     term over F(q^k);
#        T,  	     a q-hypergeometric term over F(q^k);
#        k,  	     a variable name;
#	 q,  	     an indeterminate over QQ used for the parameter q.
# Output: true,   if the following conditions are satisfied:
# 	    	  (1) T = (Qk-1)(g*H)+r*H;
#           	  (2) the Qk-quotient K of H is Qk-standard;
#		  (3) r is a Qk-remainder w.r.t. K;
#         false,  otherwise.
################################################################################

VerifyQHypergeomReduction := proc(res,T,k,q)
  local gr,H,S,K,x,y,g,r,u,v,us,un,vs,vn,u0,v0,m,nr,dr,b,cv,a,p,du,dv,lu,lv,d,\
  	i,z;
  
  #----------------------------------------------
  # compute the kernel and the shell of Qk(T)/T
  #----------------------------------------------
  gr,H := op(res);
  S := SubsPower({q^n=x,q^k=y},normal(T/H));
  K := SubsPower(q^n=x,QShiftQuotient(H,k,q,y));
  g,r := seq(SubsPower({q^n=x,q^k=y},z),z=gr);

  #----------------------------------------------------
  # check the first condition, which is equivalent to
  # S = (K*Sy-1)(g)+r
  #----------------------------------------------------
  if RationalTools:-NormalizeRatSum(S-K*subs(y=q*y,g)+g-r) <> 0 then
    print(The_1st_condition_is_violated);
    return(false);
  fi;

  #-------------------------------
  # check the second condition
  #-------------------------------
  u,v := numer(K),denom(K);
  us,un := op(QSplitFactor(u,q,y)[2..-1]);
  vs,vn := op(QSplitFactor(v,q,y)[2..-1]);
  if degree(us,y)*degree(vs,y) <> 0 or nops(QDispersionSet(un,vn,q,y)) <> 0 then
    print(The_2nd_condition_1_is_violated);
    return(false);
  fi;

  u0,v0 := seq(subs(y=0,z),z=[u,v]);
  if u0*v0 <> 0 and IsQPower(u0/v0,q,'m') and type(m,positive) then
    print(The_2nd_condition_2_is_violated);
    return(false);
  fi;
  
  #-------------------------------
  # check the third condition
  #-------------------------------
  nr,dr := numer(r),denom(r);
  g := gcd(dr,v,'b','cv');
  a,p := op(RationalTools:-ParFracDecomp(nr,[b,g],y));
  p := normal(p*cv);

  if degree(a,y) >= degree(b,y) then 
    print(The_3rd_condition_1_is_violated);
    return(false);
  fi;

  if subs(y=0,b) = 0 or QDispersionSet(b,b,q,y) <> {0} then 
    print(The_3rd_condition_2_is_violated);
    return(false);
  fi;

  if max(QDispersionSet(un,b,q,y)) >= 0 then
    print(The_3rd_condition_3_is_violated);
    return(false);
  fi;

  if min(QDispersionSet(vn,b,q,y)) <= 0 then
    print(The_3rd_condition_4_is_violated);
    return(false);
  fi;

  du := degree(u,y);
  dv := degree(v,y);
  lu := lcoeff(u,y);
  lv := lcoeff(v,y);
  if du <> dv or not IsQPower(lv/lu,q,'m') or not type(m,nonnegative) then
    d := max(du,dv);
    if degree(p,y) >= d then
      print(The_3rd_condition_5_is_violated);
      return(false);
    fi;
  elif du = 0 then
    if degree(p,y) <> m or not type(p,monomial) then
      print(The_3rd_condition_5_is_violated);
      return(false);
    fi;
  else
    r := RationalTools:-NormalizeRatSum(u*q^m*y^m-v*y^m);
    i := degree(r,y)-du;
    while i >= 0 do
      g := lcoeff(r,y)*y^i/(lu*q^i-lv);
      r := RationalTools:-NormalizeRatSum(r-u*q^i*g+v*g);
      i := degree(r,y)-du;
    od;
    dr := i+du;
    p := RationalTools:-NormalizeRatSum(p-coeff(p,y,du+m)*y^(du+m));
    if degree(p,y) >= du or coeff(p,y,dr) <> 0 then
      print(The_3rd_condition_5_is_violated);
      return(false);
    fi;
  fi;
  
  true;
end:

################################################################################
# Name: VerifyQShellReduction
# Calling sequence:
#                  VerifyQShellReduction([g,h,p],s,t,u,v,q,y)
# Input: [g,h,p],  a list with g in F(y), h a list of the form
# 	 	   [[a_1,d_1,k_1],...,[a_m,d_m,k_m]] for nonzero polynomials
#		   a_i,d_i in F[y] and positive integers k_i, and p in F[y];    
#        s,t,u,v,  four polynomials in F[y] with
#	           (1) u/v being Sy-standard;
#                  (2) gcd(s,t) = 1;
#	 q,        an indeterminate over QQ used for the parameter q;
#	 y,        a variable name.
# Output: true,  if the following conditions are satisfied:
#           	 (1) s/t = ((u/v)*Sy-1)(g) + sum(a_i/d_i^(k_i),i=1..m) + p/v;
#           	 (2) degree(a_i,y) < degree(d_i^(k_i),y) and gcd(a_i,d_i) = 1;
#	    	 (3) d_i is q-primitive, irreducible and strongly coprime with u/v;
#	    	 (4) the d_i are pairwise q-shift inequivalent;
#           	 (5) p is in the standard complement of the image of the map
#	    	     for polynomial reduction w.r.t. u/v.
#         false,  otherwise.
################################################################################

VerifyQShellReduction := proc(res,s,t,u,v,q,y)
  local g,h,p,r,i,a,d,k,un,vn,b,du,dv,lu,lv,m,dr,z;
  
  #-------------------------------
  # check the first condition
  #-------------------------------
  g,h,p := op(res);
  r := RationalTools:-NormalizeRatSum(u/v*subs(y=q*y,g)-g+add(map(z->z[1]/z[2]^(z[3]),h))+p/v-s/t);
  if r <> 0 then
    print(The_1st_condition_is_violated);
    return(false);
  fi;

  #---------------------------------------
  # check the second and third conditions
  #---------------------------------------
  for i from 1 to nops(h) do
    a,d,k := op(h[i]);
    
    if degree(a,y) >= k*degree(d,y) then 
      print(The_2nd_condition_1_is_violated);
      return(false);
    fi;
    
    if degree(RationalTools:-UnivariateGCD(a,d,y),y) <> 0 then
      print(The_2nd_condition_2_is_violated);
      return(false);
    fi;
    
    if coeff(d,y,0) = 0 then
      print(The_3rd_condition_1_is_violated);
      return(false);
    fi;
    
    if not irreduc(d) then
      print(The_3rd_condition_2_is_violated);
      return(false);
    fi;

    un := QSplitFactor(u,q,y)[-1];
    if max(QDispersionSet(un,d,q,y)) >= 0 then
      print(The_3rd_condition_3_is_violated);
      return(false);
    fi;
    
    vn := QSplitFactor(v,q,y)[-1];
    if min(QDispersionSet(vn,d,q,y)) <= 0 then
      print(The_3rd_condition_4_is_violated);
      return(false);
    fi;
  od;

  #-------------------------------
  # check the fourth condition
  #-------------------------------
  b := mul(map(z->z[2],h));
  if nops(QDispersionSet(b,b,q,y)) > 1 then 
    print(The_4th_condition_is_violated);
    return(false);
  fi;

  #-------------------------------
  # check the fifth condition
  #-------------------------------
  du := degree(u,y);
  dv := degree(v,y);
  lu := lcoeff(u,y);
  lv := lcoeff(v,y);
  if du <> dv or not IsQPower(lv/lu,q,'m') or not type(m,nonnegative) then
    d := max(du,dv);
    if degree(p,y) >= d then
      print(The_5th_condition_is_violated);
      return(false);
    fi;
  elif du = 0 then
    if degree(p,y) <> m or not type(p,monomial) then
      print(The_5th_condition_is_violated);
      return(false);
    fi;
  else
    r := RationalTools:-NormalizeRatSum(u*q^m*y^m-v*y^m);
    i := degree(r,y)-du;
    while i >= 0 do
      g := lcoeff(r,y)*y^i/(lu*q^i-lv);
      r := RationalTools:-NormalizeRatSum(r-u*q^i*g+v*g);
      i := degree(r,y)-du;
    od;
    dr := i+du;
    p := RationalTools:-NormalizeRatSum(p-coeff(p,y,du+m)*y^(du+m));
    if degree(p,y) >= du or coeff(p,y,dr) <> 0 then
      print(The_5th_condition_is_violated);
      return(false);
    fi;
  fi;
  
  true;
end:

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

################################################################################
# Name: QSpecialReduction
# Calling sequence:
#                  QSpecialReduction(s,t,u,v,q,y)
# Input: s,t,u,v,  four polynomials in F[y] with
#	           (1) u/v being Sy-standard;
#                  (2) degree(s,y) < degree(t,y) and gcd(s,t) = 1;
#		   (3) t is a power of y;
#	 q,        an indeterminate over QQ used for the parameter q;
#	 y,        a variable name.
# Output: [g,b],  where
# 	    g  - a Laurent polynomial in F[y^(-1)],
#	    b  - a polynomial in F[y],
#	    such that
#                    s/t = ((u/v)*Sy-1)(g) + b/v.
################################################################################

QSpecialReduction := proc(s,t,u,v,q,y)
  local g,b,m,u0,v0,g0;

  #-------------------------------------------
  # trivial case: s = 0
  # including the case when degree(t,y) = 0
  #-------------------------------------------
  if s = 0 then
    return([0,0]);
  fi;
  
  g := 0;
  b := normal(v*s/t);
  m := ldegree(b,y);
  u0 := coeff(u,y,0);
  v0 := coeff(v,y,0);
  
  while m < 0 do
    g0 := coeff(b,y,m)/(u0*q^m-v0)*y^m;
    g := g+g0;
    b := RationalTools:-NormalizeRatSum(b-u*subs(y=q*y,g0)+v*g0);
    m := ldegree(b,y);
  od;

  [g,b];
end:

################################################################################
# Name: QNormalReduction
# Calling sequence:
#                  QNormalReduction(s,t,u,v,q,y)
# Input: s,t,u,v,  four polynomials in F[y] with
#	           (1) u/v being Sy-reduced;
#                  (2) degree(s,y) < degree(t,y), gcd(s,t) = 1, and t(0)<>0;
#	 q,        an indeterminate over QQ used for the parameter q;
#	 y,        a variable name.
# Output: [g,h,b],  where
# 	    g - a rational function in F(y),
#	    h - a list of the form [[a_1,d_1,k_1],...,[a_m,d_m,k_m]] for
#	      	nonzero polynomials a_i,d_i in F[y] and positive integers k_i,
#	    b - a polynomial in F[y],
#	    such that
#           (1) s/t = ((u/v)*Sy-1)(g) + sum(a_i/d_i^(k_i),i=1..m) + b/v;
#           (2) degree(a_i,y) < degree(d_i^(k_i),y) and gcd(a_i,d_i) = 1;
#	    (3) d_i is q-primitive, irreducible and strongly coprime with u/v;
#	    (4) the d_i are pairwise q-shift inequivalent;.
################################################################################

QNormalReduction := proc(s,t,u,v,q,y)
  local c,tfac,m,dL,i,d,kl,aL,g,h,b,n,a,k,j,kk,l,gg,aa,bb,z;

  #-------------------------------------------
  # trivial case: s = 0
  # including the case when degree(t,y) = 0
  #-------------------------------------------
  if s = 0 then
    return([0,[],0]);
  fi;

  #---------------------------------------------------
  # compute a strong Sy-factorization of t w.r.t. u/v
  #---------------------------------------------------
  c,tfac := op(QShiftFactorization(t,q,y,[u,v]));

  #---------------------------------------------------
  # compute the partial fraction decomposition of s/t 
  # w.r.t the above strong Sy-factorization of t
  #---------------------------------------------------
  m := nops(tfac);
  dL := [];
  for i from 1 to m do
    d,kl := op(tfac[i]);
    dL := [op(dL),op(map(z->subs(y=q^(z[2])*y,d^(z[1])),kl))];
  od;
  aL := RationalTools:-ParFracDecomp(s/c,dL,y);

  #---------------------------------------------------
  # perform the normal reduction to each fraction
  #---------------------------------------------------
  g := 0; h := []; b := 0;
  for i from 1 to m do
    d,kl := op(tfac[i]);
    n := nops(kl);
    a := 0;
    k := 0;
    for j from 1 to n do
      kk,l := op(kl[j]);
      gg,aa,kk,bb := op(LocalQNormalReduction(aL[j],d,kk,l,u,v,q,y));
      if aa <> 0 then
        if kk >= k then
	  a := RationalTools:-NormalizeRatSum(a*d^(kk-k)+aa);
	  k := kk;
	else
	  a := RationalTools:-NormalizeRatSum(a+aa*d^(k-kk));
	fi;
      fi;
      g := g+gg;
      b := RationalTools:-NormalizeRatSum(b+bb);
    od;
    if a <> 0 then
      h := [op(h),[a,d,k]];
    fi;
    aL := aL[n+1..-1];
  od;

  [g,h,b];
end:

################################################################################
# Name: LocalQNormalReduction
# Calling sequence:
#                  LocalQNormalReduction(c,d,k,l,u,v,q,y)
# Input: c,d,  two polynomials in F[y] with
# 	       (1) 0 <= degree(c,y) < degree(d^k,y);
#	       (2) d being q-primitive, irreducible and
#	       	   strongly coprime with u/v;
# 	 k,    a positive integer;
#	 l,    an integer;
#	 u,v,  two polynomials in F[y] with v <> 0;
#	 q,    an indeterminate over QQ used for the parameter q;
#	 y,    a variable name.
# Output: [g,a,kk,b],  where
# 	    g   - a rational function in F(y),
#	    a,b - two polynomials in F[y] with degree(a,y) < degree(d^kk,y)
#	    	  and gcd(a,d^kk) = 1,
#	    kk  - a nonnegative integer no more than k,
#	    such that
#                    c/Sy^l(d^k) = ((u/v)*Sy-1)(g) + a/d^kk + b/v.
################################################################################

LocalQNormalReduction := proc(c,d,k,l,u,v,q,y)
  local g,kk,b,dd,a,di,pv,vi,i,ki,s,t,pu,cu,ui;

  #------------------------
  # trivial case: l = 0
  #------------------------
  if l = 0 then return([0,c,k,0]); fi;

  #----------------------------------------
  # if l > 0, then gcd(u,Sy^l(d^k)) = 1,
  # else gcd(v,Sy^(l+1)(d^k)) = 1
  #----------------------------------------
  g := 0;
  kk := k;
  b := 0;

  if l > 0 then
    dd := d;
    a := c;
    di := subs(y=q^l*y,dd);
    pv := v;
    vi := 1;
    for i from l by -1 to 1 do
      ki := kk;
      while ki > 0 and RationalTools:-IsDivided(pv,di,y,'pv') do
        ki := ki-1;
      od;
      a := normal(pv*vi*a);
      if ki = 0 then
        return([g,0,0,RationalTools:-NormalizeRatSum(b+a)]);
      else
        s,t := RationalTools:-ExtendedEuclidean(a,u,di^ki,y);
	a := subs(y=q^(-1)*y,s);
	vi := normal(vi*di^(kk-ki));
	kk := ki;
	di := subs(y=q^(i-1)*y,dd);
	g := g+a/di^kk;
	b := RationalTools:-NormalizeRatSum(b+t);
      fi;
    od;
  else
    dd := subs(y=q^l*y,d);
    a := normal(c*denom(dd)^k);
    di := numer(dd);
    pu := primpart(u,y,'cu');
    ui := 1;
    for i from l to -1 do
      ki := kk;
      g := g-a/di^kk;
      di := subs(y=q*y,di);
      while ki > 0 and RationalTools:-IsDivided(pu,di,y,'pu') do
       ki := ki-1;
      od;
      a := normal(cu*pu*ui*subs(y=q*y,a));
      if ki = 0 then
        return([g,0,0,RationalTools:-NormalizeRatSum(b+a)]);
      else
        s,t := RationalTools:-ExtendedEuclidean(a,v,di^ki,y);
	a := s;
	b := RationalTools:-NormalizeRatSum(b+t);
	ui := normal(ui*di^(kk-ki));
	kk := ki;
      fi;
    od;
    a := normal(a/denom(dd)^kk);
  fi;

  [g,a,kk,b];
end:

################################################################################
# Name: QSimpleReduction
# Calling sequence:
#                  QSimpleReduction(h,p,u,v,N,q,x,y)
# Input: h,      a list of the form [[a_1,d_1,k_1],...,[a_m,d_m,k_m]] with
# 	     	 (1) a_i,d_i in F[y]\{0} and k_i positive integers;
#	     	 (2) degree(a_i,y) < degree(d_i^(k_i),y) and gcd(a_i,d_i) = 1;
#	     	 (3) d_i being q-primitive, irreducible and strongly coprime
#		     with u/v;
#        p,u,v,  three polynomials in F[y] with
#	 	 (1) u/v being Sy-standard;
#		 (2) p in the standard complement of the image of the map
#	    	     for polynomial reduction w.r.t. u/v.
#	 N,      a nonzero rational function in F(y) with
#	 	    Sy(N)/N = Sx(u/v)/(u/v);
#	 q,      an indeterminate over QQ used for the parameter q;
#	 x,y,    two variable names.
# Output: [g,hh,pp],  where
# 	    g  - a rational function in F(y),
#	    hh - a list of the form [[aa_1,dd_1,k_1],...,[aa_m,dd_m,k_m]] for
#	      	 nonzero polynomials aa_i,dd_i in F[y],
#	    pp - a polynomial in F[y],
#	    such that
#           (1) Sx(sum(a_i/d_i^(k_i),i=1..m)+p/v)*N
#	    	= ((u/v)*Sy-1)(g) + sum(aa_i/dd_i^(k_i),i=1..m) + pp/v;
#           (2) degree(aa_i,y) < degree(dd_i^(k_i),y) and gcd(aa_i,dd_i) = 1;
#	    (3) dd_i is q-primitive, irreducible and strongly coprime with u/v;
#	    (4) the dd_i are pairwise q-shift inequivalent;
#           (5) pp is in the standard complement of the image of the map
#	    	for polynomial reduction w.r.t. u/v.
################################################################################

QSimpleReduction := proc(h,p,u,v,N,q,x,y)
  local NN,r,g,h0,pp,hh,i,a,d,k,gi,hi,pi;

  #----------------------------------------------------
  # apply the shell reduction to Sx(p/v)*N w.r.t. u/v
  #----------------------------------------------------
  NN := numer(N)/factor(denom(N));
  r := normal(subs(x=q*x,p/v)*NN);
  g,h0,pp := op(QShellReduction(numer(r),denom(r),u,v,q,y));
  
  #---------------------------------------------------------
  # apply the shell reduction to each Sx(a_i/d_i^(k_i))*N
  # w.r.t. u/v and then combine all
  #---------------------------------------------------------
  hh := table([]);
  for i from 1 to nops(h) do
    a,d,k := op(h[i]);
    r := normal(subs(x=q*x,a/d^k)*NN);
    gi,hi,pi := op(QShellReduction(numer(r),denom(r),u,v,q,y));
    g := g+gi;
    hh[i] := op(hi);
    pp := RationalTools:-NormalizeRatSum(pp+pi);
  od;
  hh := convert(hh,list);

  [g,hh,pp];
end:

################################################################################
# Name: LinearizeQRemainders
# Calling sequence:
#                  LinearizeQRemainders(h,d,tp,u,v,q,y)
# Input: h,    a list of the form [[a_1,d_1,k_1],...,[a_m,d_m,k_m]] with
# 	       (1) a_i,d_i in F[y]\{0} and k_i positive integers;
#	       (2) degree(a_i,y) < degree(d_i^(k_i),y) and gcd(a_i,d_i) = 1;
#	       (3) d_i being q-primitive, irreducible and strongly coprime
#	       	   with u/v;
#        d,    a list of the form [d0_1,...,d0_m] with d0_i being q-primitive,
#	       irreducible and strongly coprime with u/v;
#	 tp,   a list of the form [[lambda_1,mu_1],...,[lambda_m,mu_m]] with
#	       (1) lambda_i,mu_i integers;
#	       (2) mu_i > 0 and gcd(lambda_i,mu_i) = 1,
#	       such that d0_i is q-integer linear of type (lambda_i,mu_i);
#        u,v,  two polynomials in F[y] with u/v being Sy-standard;
#	 q,    an indeterminate over QQ used for the parameter q;
#	 y,    a variable name.
# Output: [g,hh,p],  where
# 	   g  - a rational function in QQ(q)(x,y),
#	   hh - a list of the form [[aa_1,dd_1,k_1],...,[aa_m,dd_m,k_m]] for
#	      	 nonzero polynomials aa_i,dd_i in F[y],
#	   p  - a polynomial in F[y],
#	   such that
#          (1) sum(a_i/d_i^(k_i),i=1..m)
#	       = ((u/v)*Sy-1)(g) + sum(aa_i/dd_i^(k_i),i=1..m) + p/v;
#          (2) degree(aa_i,y) < degree(dd_i^(k_i),y) and gcd(aa_i,dd_i) = 1;
#	   (3) dd_i is q-primitive, irreducible, strongly coprime with u/v;
#	   (4) the dd_i are pairwise q-shift inequivalent;
#	   (5) dd_1*...*dd_m is Sy-coprime with d0_1*...*d0_m;
#          (6) p is in the standard complement of the image of the map
#	       for polynomial reduction w.r.t. u/v.
################################################################################

LinearizeQRemainders := proc(h,d,tp,u,v,q,y)
  local g,p,hh,n,id,m,i,a,dd,k,lm,j,l,co,d0,gg,b,f;
  
  #--------------------------------------------------
  # transform h to a Sy-remainder whose significant
  # denominator is Sy-coprime with d0_1*...*d0_m
  #--------------------------------------------------
  g := 0;
  p := 0;
  hh := table([]);
  n := nops(d);
  id := [seq(i,i=1..n)];
  m := n;
  for i from 1 to n do
    a,dd,k := op(h[i]);
    lm := tp[i];
    j := 1;
    while j <= m and\
          (lm <> tp[id[j]] or not IsQShiftEquivalent(dd,d[id[j]],q,y,'l','co'))\
    do
      j := j+1;
    od;
    if j <= m then
      d0 := d[id[j]];
      gg,a,k,b := op(LocalQNormalReduction(normal(a/co),d0,k,l,u,v,q,y));
      g := g+gg;
      hh[i] := [a,d0,k];
      p := RationalTools:-NormalizeRatSum(p+b);
      id := [op(id[1..j-1]),op(id[j+1..-1])];
      m := m-1;
    else
      hh[i] := [a,dd,k];
    fi;
  od;
  hh := convert(hh,list);
  
  #-------------------------------------------------
  # apply the polynomial reduction to p w.r.t. u/v
  #-------------------------------------------------
  f,p := op(QPolynomialReduction(p,u,v,q,y));
  g := g+f;

  [g,hh,p];
end:

################################################################################
# Name: FindLinearRelation
# Calling sequence:
#                  FindLinearRelation([h_1,...,h_m],[p_1,...,p_m],y)
# Input: [h_1,...,h_m],  a list of lists, where h_i is of the form
# 	       	       	 [[a_{i1},d_{i1},k_{i1}],...,[a_{in},d_{in},k_{in}]]
#		       	 with
# 	      	       	 (1) a_{ij},d_{ij} in F[y]\{0} and k_{ij} positive integers;
#	      	       	 (2) degree(a_{ij},y) < degree(d_{ij}^(k_{ij}),y)
#	      	       	     and gcd(a_{ij},d_{ij}) = 1;
#	      	       	 (3) d_{ij} being q-primitive and irreducible;
#			 (4) d_{i1},...,d_{in} are pairwise q-shift inequivalent;
# 	 [p_1,...,p_m],  a list of polynomials, where p_i is in F[y];
#        y,   	    	 a variable name.
# Output: [c_1,...,c_m],  a list of elements in F, not all zero, such that
#              		  sum(c_i*sum(a_{ij}/d_{ij}^(k_{ij}),j=1..n),i=1..m)=0,
#              		  and sum(c_i*p_i,i=1..m) = 0.
################################################################################

FindLinearRelation := proc(hL,pL,y)
  local cL,m,vars,rL,eqns,sys,mat,i,gp,p,r,sol,c,z,zz;
  
  #----------------------------------------
  # generate a matrix of coefficients
  #----------------------------------------
  cL := table([]);
  m := nops(hL);
  vars := [seq(cL[i],i=1..m)];
  rL := map(z->RationalTools:-NormalizeRatSum(add(map(zz->zz[1]/zz[2]^(zz[3]),z))),hL);
  eqns := [add(cL[i]*rL[i],i=1..m),add(cL[i]*pL[i],i=1..m)];
  sys := `union`(op(map(z->
	 convert(PolynomialTools:-CoefficientList(collect(numer(z),y),y),set)\
	 minus {0},eqns)));
  mat := LinearAlgebra:-GenerateMatrix(sys,vars)[1];

  #------------------------------
  # find the linear dependency
  #------------------------------
  try
    for i from 1 to 3 do
      gp := rand(6..100);
      p := ithprime(gp());
      r := LowerBoundForRank(mat,p);
      if r = m then
        return([]);
      fi;
    od;
  catch:
  end try;
  sol := SolveTools:-Linear(sys,vars);

  c := rhs(op(select(z->lhs(z)=cL[m],sol)));
  if c = 0 then
    return([]);
  else
    assign(map(z->lhs(z)=normal(rhs(z)/c),sol));
    convert(cL,list);
  fi;
end:

################################################################################
# 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: r,  a nonnegative integer less than or equal to the rank of M.         
################################################################################

LowerBoundForRank := proc(M,p)
  local rn,cn,vars,d,i,di,j,MM,g,n,e,s,N,r;
  
  #--------------------------------
  # collect info about the matrix
  #--------------------------------
  rn := op(1,M)[1]:
  cn := op(1,M)[2]:
  vars := [op(indets(M))];

  #---------------------------------------------
  # compute the common denominator for each row
  #---------------------------------------------
  d := []:
  for i from 1 to rn do
    di := []:
    for j from 1 to cn do
      di := [op(di),denom(M[i,j])]:
    od:
    d := [op(d),Lcm(op(di)) mod p]:
  od:

  MM := Matrix(rn,cn):
  for i from 1 to rn do
    for j from 1 to cn do
      MM[i,j] := Normal(d[i]*M[i,j]) mod p:
    od;
  od;
  
  #-------------------------------
  # choose two evaluation points
  #-------------------------------
  g := rand(1..p-1);
  n := nops(vars);
  e := [seq(g(),i=1..n)];
  s := {seq(vars[i]=e[i],i=1..n)};

  #-------------------------
  # perform the 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:
    od;
  od;

  #------------------------
  # compute the rank
  #------------------------
  LinearAlgebra:-Modular:-RowReduce(p,N,rn,cn,cn,0,0,r,0,0,false);

  r;
end:

################################################################################
# Name: QUpperBound
# Calling sequence:
#                  QUpperBound(tab,y)
# Input: tab,  a table whose indices are distinct and with entries of the form
# 	               tab[[lambda,mu]]
#		     = [
#		        [P_1,[[k_{11},l_{11}],...,[k_{1n_1},l_{1n_1}]]],
#			       ...
#	      		[P_m,[[k_{m1},l_{m1}],...,[k_{mn_m},l_{mn_m}]]]
#		       ],
#              where lambda,mu are integers with mu > 0 and gcd(lambda,mu) = 1,
#	       each P_i is in QQ(q)[y] with degree(P_i,y) > 0, and k_{ij},l_{ij}
#	       are integers with k_{ij} > 0;
#	 y,    a variable name.
# Output: ord,  the sum of all
# 	  	    sum(max(k_{i1},...,k_{in_i})*mu^2*degree(P_i,y),i=1..m)
#		for mu running through the indices of tab.
################################################################################

QUpperBound := proc(tab,y)
  local ord,id,i,lambda,mu,z,zz;
  
  ord := 0;
  id := [indices(tab,'nolist')];
  for i from 1 to nops(id) do
    lambda,mu := op(id[i]);
    ord := normal(ord+\
           add(map(z->mu^2*max(map(zz->zz[1],z[2]))*degree(z[1],y),tab[[lambda,mu]])));
  od;

  ord;
end:

################################################################################
# Name: QLowerBound
# Calling sequence:
#                  QLowerBound(tab)
# Input: tab,  a table whose indices are distinct and with entries of the form
# 	               tab[[lambda,mu]]
#		     = [
#		        [P_1,[[k_{11},l_{11}],...,[k_{1n_1},l_{1n_1}]]],
#			       ...
#	      		[P_m,[[k_{m1},l_{m1}],...,[k_{mn_m},l_{mn_m}]]]
#		       ],
#              where lambda,mu are integers with mu > 0 and gcd(lambda,mu) = 1,
#	       each P_i is in QQ(q)[y] with degree(P_i,y) > 0, and k_{ij},l_{ij}
#	       are integers with k_{ij} > 0.
# Output: ord,  the nonnegative integer which is equal to
# 	  	max_{i=1..m,j=1..n_i}(
#		 min(rho in ZZ^+ | k_{ij}<=k_{ij'} and
#		  l_{ij}=l_{ij'}+lambda*rho mod mu for some j' in {1,...,n_i}))
#               for mu running through the indices of tab.
################################################################################

QLowerBound := proc(tab)
  local S,id,i,lambda,mu,rho,z;

  S := {};
  id := [indices(tab,'nolist')];
  for i from 1 to nops(id) do
    lambda,mu := op(id[i]);
    igcdex(lambda,mu,'rho');
    S := S union\
         `union`(op(map(z->MinPosIntSet(z[2],mu,rho),tab[[lambda,mu]])));
  od;

  max(S);
end:

################################################################################
# Name: MinPosIntSet
# Calling sequence:
#                  MinPosIntSet([[k_1,l_1],...,[k_m,l_m]],mu,rho)
# Input: [[k_1,l_1],...,[k_m,l_m]],  a list of integer pairs, where k_i,l_i
# 	 			     are integers and k_i > 0;
#	 mu,   			     a positive integer;
#	 rho,  			     a nonzero integer with gcd(rho,mu) = 1.
# Output: {n_1,...,n_m},  where
# 	   n_i - the minimal positive integer such that
#	       	        mu divides rho*(l_i-l_j)-n_i
#		 for some j in {1,...,m}.
################################################################################

MinPosIntSet := proc(kl,mu,rho)
  local m,nL,i,k,l,j,kk,ll,c,n;

  m := nops(kl);
  nL := [seq(mu,i=1..m)];
  for i from 1 to m do
    k,l := op(kl[i]);
    for j from i to m do
      kk,ll := op(kl[j]);
      c := l-ll;
      if k <= kk then
        n := modp(c*rho,mu);
	if n <> 0 then nL[i] := min(nL[i],n); fi;
      fi;
      if k >= kk then
        n := modp(-c*rho,mu);
	if n <> 0 then nL[j] := min(nL[j],n); fi;
      fi;
    od;
  od;

  convert(nL,set);
end:

################################################################################
# This code is based on SubsPowers in Maple_QDifferenceEquations.
# Name: SubsPower
# Calling sequence:
#                  SubsPower(sub,expr)
# Input: sub,   a single substitution or a list/set of substitutions;  
#        expr,  any expression.
# Output: nexpr,  the expression obtained by substituting sub into expr.
################################################################################

SubsPower := proc(sub,expr)
  local s,nexpr,names,sL,z;

  if type(sub,equation) then
    s := {sub};
  elif type(sub,{list(equation),set(equation)}) then
    s := {op(sub)};
  fi;

  nexpr := subs(s,expr);

  #----------------------------------------
  # check if all powers were substituted
  #----------------------------------------
  s := select(type,s,name^name=anything);
  names := map(z->op(2,lhs(z)),s) intersect indets(nexpr,name);
  if names <> {} then
    s := select(unapply('has'('op'(1,z),names),z), s);
    sL := [seq([op(op(1,z)),op(2,z)],z=s)];
    names := indets(nexpr,anything^anything) minus\
             indets(nexpr,{anything^(-1),name^name,name^integer,\
	                   (name^name)^integer,identical(I)});
    nexpr := eval(nexpr,[seq(z=LocalSubsPower(sL,Normalizer(op(1,z)),\
					      expand(normal(op(2,z)))),z=names)]);
  fi;

  nexpr;
end:

################################################################################
# Name: LocalSubsPower
# Calling sequence:
#                  LocalSubsPower(sL,b,e)
# Input: sL,   a list of the form [[b_1,e_1,x_1],[b_2,e_2,x_2],...]
# 	        representing the substitutions
#			     sub = {b_1^e_1=x_1,b_2^e_2=x_2,...};
#        b,e,  any expressions.
# Output: expr,  the expression obtained by substituting sub into b^e.
################################################################################

LocalSubsPower := proc(sL,b,e)
  local bL,z;

  bL := map(z->z[1],sL);
  if type(e,`+`) then
    return `*`(seq(procname(sL,b,z),z=e));
  elif type(b,`*`) then
    return map2(procname,sL,b,e);
  elif type(b,`^`) then
    return procname(sL,op(1,b),expand(normal(op(2,b)*e)));
  elif type(b,name) and member(b,bL,'z') then
    z := sL[z];
    if denom(numer(e)/z[2]) = 1 then
      return expand(z[3]^(e/z[2]));
    else
      return expand(b^e);
    fi;
  elif hastype(b,`^`) and has(b,bL) then
    return subs([seq(z=procname(sL,op(z)),z=indets(b,`^`) minus\
    	    	                            indets(b,anything^(-1)))],b)^e;
  else
    return expand(b^e);
  fi;
  
end:

################################################################################
# Name: QNormalRNF
# Calling sequence:
#                  QNormalRNF(f,q,y)
# Input: f,  a q-primitive rational function in F(y);
# 	 q,  an indeterminate over QQ used for the parameter q;
#        y,  a variable name.
# Output: [u,v,s,t],  where
#           u,v,s,t - four q-primitive polynomials in F[y],
#	    such that
#           (1) f = K*Sy(S)/S with K = u/v and S = s/t;
#	    (2) K is Sy-reduced;
#           (3) gcd(s,t) = gcd(u,s*Sy(t)) = gcd(v,Sy(s)*t) = 1.
################################################################################

QNormalRNF := proc(f,q,y)
  local ff,u,v,s,t,DS,pDS,nDS,i,k,p,pk,j,st;

  #------------------------------------------
  # trivial case: either f or 1/f is in F[y]
  #------------------------------------------
  ff := normal(f);
  u := numer(ff);
  v := denom(ff);
  s := 1;
  t := 1;
  if type(u,freeof(y)) or type(v,freeof(y)) then
    return([u,v,s,t]);
  fi;
  
  #-----------------------------------------
  # compute the dispersion set of u and v
  #-----------------------------------------
  DS := QDispersionSet(u,v,q,y);
  pDS,nDS := selectremove(type,DS,positive);

  #--------------------------------------
  # deal with positive shifts
  #--------------------------------------
  pDS := sort(convert(pDS,list));
  for i from 1 to nops(pDS) do
    k := pDS[i];
    p := RationalTools:-UnivariateGCD(u,subs(y=q^k*y,v),y,'u');
    pk := subs(y=q^(-k)*y,p);
    v := normal(v/pk);
    s := s*mul(subs(y=q^(-j)*y,p),j=1..k-1)*pk;
  od;

  #--------------------------------------
  # deal with negative shifts
  #--------------------------------------
  nDS := sort(convert(nDS,list));
  for i from nops(nDS) by -1 to 1 do
    k := nDS[i];
    p := RationalTools:-UnivariateGCD(u,subs(y=q^k*y,v),y,'u');
    pk := subs(y=q^(-k)*y,p);
    v := normal(v/pk);
    t := t*mul(subs(y=q^j*y,p),j=0..-k-1);
  od;
  st := normal(s/t);
  
  [u,v,numer(st),denom(st)];
end:

################################################################################
# Name: IsIrredQIntegerLinear
# Calling sequence:
#                  IsIrredQIntegerLinear(p,x,y,'tp','po')
# Input: p,      an irreducible polynomial in QQ(q)[x,y] with degree(p,y) > 0;
# 	 x,y,    two variable names;
#	 tp,po,  (optional) two unevaluated names.
# Output: true,   if p is q-integer linear, and then the optional 4th and 5th 
# 	  	  arguments are assigned the list [lambda,mu] and a polynomial P,
#		  respectively, where lambda,mu are integers with mu > 0 and
#		  gcd(lambda,mu) = 1, and P is in QQ(q)[y] such that
#		                   p = numer(P(x^lambda*y^mu));
# 	  false,  otherwise.
################################################################################

IsIrredQIntegerLinear := proc(p,x,y)
  local supp,xsupp,ysupp,n,s,lm,lambda,mu,i,a,b,P,z;

  #------------------------------------
  # trivial case: p is free of x
  #------------------------------------
  if type(p,freeof(x)) then
    if nargs = 4 then assign(args[4]=[0,1]); fi;
    if nargs = 5 then assign(args[4]=[0,1],args[5]=p); fi;
    return true;
  fi;
  
  #------------------------------------------------------
  # compute the support of p w.r.t. x,y and sort it in
  # lexicographic order with y > x
  #------------------------------------------------------
  supp := Groebner:-Support(p,[y,x]);
  supp := map(z->[degree(z,x),degree(z,y)],supp);

  #------------------------------------------------------
  # check if the x-exponents, as well as y-exponents,
  # of all nonzero terms of p are distinct
  #------------------------------------------------------
  xsupp := convert(map(z->z[1],supp),set);
  ysupp := convert(map(z->z[2],supp),set);
  n := nops(supp);
  if nops(xsupp) < n or nops(ysupp) < n then
    return false;
  fi;

  #---------------------------------------------
  # find the possible q-integer linear type
  # Note that s[2] = 0 since p is irreducible
  #---------------------------------------------
  s := supp[-1];
  supp := map(z->z-s,supp[1..-2]);
  lm := `/`(op(supp[-1]));
  lambda := numer(lm);
  mu := denom(lm);

  #-----------------------------------------------
  # check if all pairs in supp are of the form
  # (lambda*i,mu*i) for some integer i
  #-----------------------------------------------
  for i from n-2 by -1 to 1 do
    a,b := op(supp[i]);
    if RationalTools:-NormalizeRatSum(mu*a-lambda*b) <> 0 then
      return false;
    fi;
  od;

  if nargs = 4 then assign(args[4]=[lambda,mu]); fi;
  if nargs = 5 then
    P := normal(subs(x=x^mu,y=y*x^(-lambda),p/x^s[1]));
    assign(args[4]=[lambda,mu],args[5]=normal(subs(y=y^(1/mu),P)));
  fi;
  
  true;
end:

################################################################################
# Name: IsQShiftEquivalent
# Calling sequence:
#                  IsQShiftEquivalent(a,b,q,y,'m','co');
# Input: a,b,  two q-primitive polynomials in F[y] of positive degrees in y;
# 	 q,    an indeterminate over QQ used for the parameter q;
#        y,    a variable name;
#	 m,co,  (optional) unevaluated names.
# Output: true,   if a = c*Sy^k(b) in F for some c in F and some integer k, 
# 	  	  and then the optional 5th and 6th arguments m and co are
#		  assigned the integer k and the element c, respectively.
#         false,  otherwise.
################################################################################

IsQShiftEquivalent := proc(a,b,q,y)
  local d,c,m;

  #---------------------
  # check the degree
  #---------------------
  d := degree(a,y);
  if d <> degree(b,y) then return false; fi;

  #---------------------
  # find the candidate
  #---------------------
  c := normal(coeff(a,y,0)/coeff(b,y,0));
  if IsQPower(lcoeff(a,y)/lcoeff(b,y)/c,q,'m') then
    m := normal(m/d); 
    if type(m,integer) and RationalTools:-NormalizeRatSum(a-c*subs(y=q^m*y,b)) = 0 then
      if nargs = 5 then assign(args[5]=m); fi;
      if nargs = 6 then assign(args[5]=m,args[6]=c); fi;
      return true;
    fi;
  fi;

  false;
end:

################################################################################
# Name: IsQPower
# Calling sequence:
#                  IsQPower(c,q,'m')
# Input: c,  an element in F;
#        q,  an indeterminate over QQ used for the parameter q;
#	 m,  (optional) an unevaluated name.
# Output: true,   if c = q^k for some integer k, and then the optional 
# 	  	  3rd argument m is assgined the integer k.
# 	  false,  otherwise.
################################################################################

IsQPower := proc(c,q)
  local cc;

  cc := normal(c);
  if cc = 1 then
    if nargs = 3 then assign(args[3]=0); fi;
    return true;
  elif cc = q then
    if nargs = 3 then assign(args[3]=1); fi;
    return true;
  elif type(cc,`^`) and op(1,cc) = q and type(op(2,cc),integer) then
    if nargs = 3 then assign(args[3]=op(2,cc)); fi;
    return true;
  fi;
  
  false;
end:

end module:

######## RationalTools

RationalTools := module()
  description "Tools for dealing with rational functions":

  export
        UnivariateGCD,
	ExtendedEuclidean,
	IsDivided,
	NormalizeRatSum,
	ParFracDecomp;
  local
	ExactDivision,
	ModularInverse,
	ListModularInverse,
	MultiPowerInverse,
	PowerInverse,
	QuadLifting;
  option package:

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

################################################################################
# Name: UnivariateGCD
# Calling sequence:
#                  UnivariateGCD(a,b,y,'cofa','cofb')
# Input: a,b,        two polynomials in F[y];
#        y,          a variable name;
# 	 cofa,cofb,  (optional) two unevaluated names.
# Output: g,  the primitive gcd of a and b in F[y], and the optional
# 	      4th and 5th arguments cofa and cofb are assigned the
#	      cofactors a/g and b/g, respectively.
################################################################################

UnivariateGCD := proc(a,b,y)
  local pa,ca,pb,cb,g,cofa,cofb;

  pa := primpart(a,y,'ca');
  pb := primpart(b,y,'cb');
  if nargs = 5 then
    g := gcd(pa,pb,'cofa','cofb');
    assign(args[4]=normal(ca*cofa),args[5]=normal(cb*cofb));
  elif nargs = 4 then
    g := gcd(pa,pb,'cofa');
    assign(args[4]=normal(ca*cofa));
  else
    g := gcd(pa,pb);
  fi;
  
  g;
end:

################################################################################
# Name: ExtendedEuclidean
# Calling sequence:
#                  ExtendedEuclidean(c,a,b,y,'half')
# Input: c,a,b,  three polynomials in F[y] with gcd(a,b)=1;
#        y,      a variable name;
#	 'half'  (optional) an argument name used for specifying the output.
# Output: s,t,  if the optional 5th argument 'half' is missing, where
# 	   s,t - polynomials in F[y] with degree(s,y) < degree(b,y)
#	   such that
#          	    c = s*a + t*b.
#         s,    otherwise.
################################################################################

ExtendedEuclidean := proc(c,a,b,y)
  local s,t;

  #-------------------------------------------------------
  # compute s in F[y] such that degree(s,y) < degree(b,y)
  # and s*a = c mod b
  #-------------------------------------------------------  
  s := ModularInverse(a,b,y);
  s := rem(normal(s*c),b,y);
  if nargs = 5 and type(args[5],'identical(half)') then
    return s;
  fi;

  #--------------------------------------------
  # compute t in F[y] such that c = s*a + t*b
  #--------------------------------------------
  t := ExactDivision(NormalizeRatSum(c-s*a),b,y);

  s,t;
end:

################################################################################
# Name: IsDivided
# Calling sequence:
#                  IsDivided(a,b,y,'c')
# Input: a,b,  two polynomials in F[y];
# 	 y,    a variable name;
#	 c,    (optional) an unevaluated name.
# Output: true,   if a is divided by b, and then the optional 4th argument
# 	  	  is assigned the quotent of a by b;
#         false,  otherwise.
################################################################################

IsDivided := proc(a,b,y)
  local c,r;

  c := quo(a,b,y,'r');
  if r <> 0 then
    return false;
  fi;

  if nargs = 4 then assign(args[4]=c); fi;

  true;
end:

################################################################################
# Name: NormalizeRatSum
# Calling sequence:
#                  NormalizeRatSum(r)        
# Input: r,  a rational function in F(y);
# Output: rr,  the reduced form of r.
################################################################################

NormalizeRatSum := proc(r)
  local rL,rr,i,nr,dr,f,nf,df,g,cdr,cdf;

  rL := `if`(type(r,`+`),[op(r)],[r]);
  rr := normal(rL[1]);
  for i from 2 to nops(rL) do
    nr,dr := numer(rr),denom(rr);
    f := normal(rL[i]);
    nf,df := numer(f),denom(f);
    g := gcd(dr,df,'cdr','cdf');
    rr := normal(nr*cdf+nf*cdr)/(g*cdr*cdf);
  od;
  
  rr;
end:

################################################################################
# Name: ParFracDecomp
# Calling sequence:
#                  ParFracDecomp(a,[d_1,...,d_m],y)
# Input: a,              a polynomial in F[y];
# 	 [d_1,...,d_m],  a list of polynomials in F[y] with
#	 		 (1) degree(a,y) < degree(d_1*...*d_m,y)
#			     and gcd(a,d_1*...*d_m) = 1;
#			 (2) degree(d_i) > 0 and gcd(d_i,d_j) = 1 whenever i<>j;
#	 y,              a variable name.
# Output: [a_1,...,a_m],  where
# 	    a_i - a polynomial in F[y] with degree(a_i,y) < degree(d_i,y)
#	    	  and gcd(a_i,d_i) = 1,
#	    such that
#	          a/(d_1*...*d_m) = a_1/d_1+...+a_m/d_m.
################################################################################

ParFracDecomp := proc(a,dL,y)
  local aL,m,i,d1,j,d2,s,t,z;

  #----------------------------------------------------------
  # construct [a_1,...,a_m] by
  # a_i = a*prod(d_j,j=1..m and j<>i)^(-1) mod d_i
  #----------------------------------------------------------
  aL := map(z->rem(a,z,y),dL);
  m := nops(dL);
  for i from 1 to m-1 do
    d1 := dL[i];
    for j from i+1 to m do
      d2 := dL[j];
      s := ModularInverse(d2,d1,y);
      aL[i] := rem(aL[i]*s,d1,y);
      t := ExactDivision(NormalizeRatSum(1-s*d2),d1,y);
      aL[j] := rem(aL[j]*t,d2,y);
    od;
  od;

  aL;
end:

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

################################################################################
# Name: ExactDivision
# Calling sequence:
#                  ExactDivision(a,b,y)
# Input: a,b,  two polynomials in F[y] with b dividing a;
# 	 y,    a variable name.
# Output: c,  a polynomial in F[y] such that a = b*c.
################################################################################

ExactDivision := proc(a,b,y)
  local pa,ca,pb,cb,c;

  pa := primpart(a,y,'ca');
  pb := primpart(b,y,'cb');
  c := quo(pa,pb,y);

  normal(c*ca/cb);
end:

################################################################################
# Name: ModularInverse
# Calling sequence:
#                  ModularInverse(f,g,y)
# Input: f,g,  two coprime polynomials in F[y];
#        y,    a variable name.
# Output: r,  a polynomial in F[y] of degree less than degree(g,y) such that
# 	      g divides (1-r*f).
################################################################################

ModularInverse := proc(f,g,y)
  local cf,pf,ffac,pg,gfac,r,z;

  #----------------------------------------------------------
  # trivial case: degree(f,y) = 0 including the case g = 0 
  #----------------------------------------------------------
  if degree(f,y) = 0 then return(1/f); fi;

  #----------------------------------------------------------
  # trivial case: degree(g,y) = 0 including the case f = 0
  #----------------------------------------------------------
  if degree(g,y) = 0 then return(0); fi;

  #--------------------------------------------------------
  # explore the possible factorization forms of
  # the primitive parts of f and g
  #--------------------------------------------------------
  cf := content(f,y,'pf');
  ffac := `if`(type(pf,`*`),[op(pf)],[pf]);
  ffac := map(z->`if`(type(z,`^`),[op(z)],[z,1]),ffac);

  pg := primpart(g,y);
  gfac := `if`(type(pg,`*`),[op(pg)],[pg]);
  gfac := map(z->`if`(type(z,`^`),[op(z)],[z,1]),gfac);

  #---------------------------------------------------------------
  # compute the inverse of pf modulo pg using modular arithmetics
  #---------------------------------------------------------------
  r := ListModularInverse(ffac,gfac,y);

  #-------------------------------------------
  # incorporate the content of f and output
  #-------------------------------------------
  normal(r/cf);
end:

################################################################################
# Name: ListModularInverse
# Calling sequence:
#   ListModularInverse([[f_1,d_1],...,[f_m,d_m]],[[g_1,e_1],...,[g_n,e_n]],y)
# Input: [[f_1,d_1],...,[f_m,d_m]],[[g_1,e_1],...,[g_n,e_n]],
#	     two nonempty lists with
#	     (1) f_i,g_i in F[y]\F and gcd(f_i,g_j) = 1 for all i,j;
#	     (2) d_i,e_i positive integers;
#        y,  a variable name.
# Output: r,  a polynomial in F[y] with
# 	             degree(r,y) < degree(g_1^(e_1)*...*g_n^(e_n),y)
#	      such that
# 	    r*(f_1^(d_1)*...*f_m^(d_m)) = 1 mod (g_1^(e_1)*...*g_n^(e_n)).
################################################################################

ListModularInverse := proc(ffac,gfac,y)
  local s,g,r,i;

  #-----------------------------------------------------------------------
  # compute the inverse of each f_i^(d_i) modulo g_1^(e_1)*...*g_n^(e_n)
  #-----------------------------------------------------------------------
  s := map(MultiPowerInverse,ffac,gfac,y,'g');

  #--------------------------------------------------------------------------
  # accumulate the modular inverses iteratively using
  # (r*s[i])*(f_1^(d_1)*...*f_i^(d_i)) = 1 mod g_1^(e_1)*...*g_n^(e_n))
  # with r*(f_1^(d_1)*...*f_{i-1}^(d_{i-1})) = 1 mod g_1^(e_1)*...*g_n^(e_n))
  # for i = 2,...,m
  #--------------------------------------------------------------------------
  r := s[1];

  for i from 2 to nops(s) do
    r := rem(r*s[i],g,y);
  od;

  r;
end:

################################################################################
# Name: MultiPowerInverse
# Calling sequence:
#              MultiPowerInverse([f,d],[[g_1,e_1],...,[g_n,e_n]],y,'g')
# Input: [f,d],  a list with f in F[y]\F and d a positive integer;
# 	 [[g_1,e_1],...,[g_n,e_n]],
#		 a nonempty list with g_i in F[y]\F and e_i positive
#		 integers such that gcd(f,g_i) = 1 for all i;
#        y,   	 a variable name;
#	 g,      (optional) an unevaluated name.
# Output: r,  a polynomial in F[y] with
# 	             degree(r,y) < degree(g_1^(e_1)*...*g_n^(e_n),y)
#	      such that
# 	           r*f^d = 1 mod (g_1^(e_1)*...*g_n^(e_n)).
#	      The optional 4th argument g is assigned the product
# 	                  g_1^(e_1)*...*g_n^(e_n).
################################################################################

MultiPowerInverse := proc(ffac,gfac,y)
  local s,r,f,g,i,t;

  #----------------------------------------------------
  # compute the inverse of f^d modulo each g_i^(e_i)
  #----------------------------------------------------
  s := map2(PowerInverse,ffac,gfac,y);

  #--------------------------------------------------------
  # accumulate the modular inverses iteratively using
  # (r+s[i]-r*s[i]*f^d)*f^d = 1 mod g_1^(e_1)*...*g_i^(e_i)
  # with r*f^d = 1 mod g_1^(e_1)*...*g_{i-1}^(e_{i-1})
  # for i = 2,...,n
  #--------------------------------------------------------
  r := s[1];
  f := `^`(op(ffac));
  g := `^`(op(gfac[1]));

  for i from 2 to nops(s) do
    t := s[i];
    g := normal(g*`^`(op(gfac[i])));
    r := rem(NormalizeRatSum(r+t-r*t*f),g,y);
  od;

  if nargs = 4 then assign(args[4]=g); fi;

  r;
end:

################################################################################
# Name: PowerInverse
# Calling sequence:
#                  PowerInverse([f,d],[g,e],y)
# Input: [f,d],[g,e],  two lists with
# 	 	       (1) f,g in F[y]\F and gcd(f,g)=1;
# 	 	       (2) d,e positive integers;
#        y,   	       a variable name.
# Output: r,  a polynomial in F[y] with degree(r,y) < degree(g^e,y)
# 	      such that
# 	                r*f^d = 1 mod g^e.
################################################################################

PowerInverse := proc(ffac,gfac,y)
  local f,d,g,e,r,b,br,s,gg,i;

  #-----------------------------------
  # compute the inverse of f modulo g
  #-----------------------------------
  f,d := op(ffac);
  g,e := op(gfac);
  gcdex(f,g,y,'r');

  #------------------------------------------------------
  # quadratically lift the inverse of f modulo g
  # to the inverse of f modulo g^e
  #------------------------------------------------------
  r := QuadLifting(r,f,g,e,y);

  #--------------------------
  # trivial case: d = 1
  #--------------------------
  if d = 1 then return(r); fi;

  #-------------------------------------------------------
  # compute the binary representation of d, namely find 
  # [d_0,...,d_{m-1}] with d_i in {0,1} such that
  #     d = 2^m+d_{m-1}*2^(m-1)+...+d_1*2+d_0
  #-------------------------------------------------------
  b := convert(d,binary);
  br := [];
  while b <> 1 do
    br := [op(br),irem(b,10,'b')];
  od;

  #----------------------------------------------------------------------
  # iteratively and quadratically lift the inverse of f modulo g^e 
  # to the inverse of f^d modulo g^e using
  # (1) r^2*f^(d_m*2^(m-i)+...+d_(i+1)*2) = 1 mod g^e
  # (2) (r*s)*f^(d_m*2^(m-i)+...+d_{i+1}*2+d_i) = 1 mod g^e
  # with s*f = 1 mod g^e and r*f^(d_m*2^(m-i-1)+...+d_{i+1}) = 1 mod g^e
  # for i = m-1,...,0 (note that d_m = 1)
  #-----------------------------------------------------------------------
  s := r;
  gg := g^e;
  for i from nops(br) by -1 to 1 do
    r := rem(r^2,gg,y);
    if br[i] = 1 then
      r := rem(r*s,gg,y);
    fi;
  od;

  r;
end:

################################################################################
# Name: QuadLifting
# Calling sequence:
#                  QuadLifting(s,f,g,e,y)
# Input: s,    a polynomial in F[y] with
# 	       (1) degree(s,y) < degree(g,y);
#	       (2) s*f = 1 mod g;
# 	 f,g,  two coprime polynomials in F[y]\F;
#	 e,    a positive integer;
#        y,    a variable name.
# Output: r,  a polynomial in F[y] with degree(r,y) < degree(g^e,y)
# 	      such that
# 	                r*f = 1 mod g^e.
################################################################################

QuadLifting := proc(s,f,g,e,y)
  local b,br,gg,r,c,i;

  #-----------------------
  # trivial case: e = 1
  #-----------------------
  if e = 1 then return(s); fi;

  #-------------------------------------------------------
  # compute the binary representation of e, namely find 
  # [e_0,...,e_{m-1}] with e_i in {0,1} such that
  #     e = 2^m+e_{m-1}*2^(m-1)+...+e_1*2+e_0
  #-------------------------------------------------------
  b := convert(e,binary);
  br := [];
  while b <> 1 do
    br := [op(br),irem(b,10,'b')];
  od;

  #------------------------------------------------------------------
  # iteratively and quadratically lift the inverse of f modulo g
  # to the inverse of f modulo g^e using
  # (1) (r*(2-r*f))*f = 1 mod g^(e_n*2^(n-i)+...+e_{i+1}*2)
  # (2) (s+r*(1-s*f))*f = 1 mod g^(e_n*2^(n-i)+...+e_{i+1}*2+e_i)
  # with s*f = 1 mod g and r*f = 1 mod g^(e_n*2^(n-i-1)+...+e_{i+1})
  # for i = n-1,...,0 (note that e_n = 1)
  #------------------------------------------------------------------
  gg := g;
  r := s;
  c := NormalizeRatSum(1-s*f);
  for i from nops(br) by -1 to 1 do
    gg := normal(gg^2);
    r := rem(NormalizeRatSum(r*(2-r*f)),gg,y);
    if br[i] = 1 then
      gg := normal(g*gg);
      r := rem(NormalizeRatSum(s+r*c),gg,y);
    fi;
  od;

  r;
end:

end module:
