## Copyright (c) May 2024, FZU. All rights reserved.
## Author: Hui Huang <huanghui@fzu.edu.cn>
## Load the modules: HypergeometricCT.mm, RationalTools.mm
################################################################################
				   
HypergeometricCT := module()
  description "Creative telescoping for hypergeometric terms via reduction":

  export
        HypergeomTelescoping,
	HypergeomReduction,
	ShellReduction,
	PolynomialReduction,
	OrderBounds,
	ShiftQuotient,
	RNF,
	DispersionSet,
	ShiftFactorization,
	VerifyHypergeomTelescoping,
	VerifyHypergeomReduction,
	VerifyShellReduction;
  local
	NormalReduction,
	LocalNormalReduction,
	SimpleReduction,
	LinearizeRemainders,
	FindLinearRelation,
	LowerBoundForRank,
	UpperBound,
	LowerBound,
	MinPosIntSet,
	PrimitiveRNF,
	IsIrredIntegerLinear,
	IsShiftEquivalent;
  option package:

##########################  Task  Description ##################################
# Creative telescoping problem for hypergeometric terms:
# Given: a hypergeometric term T over QQ(x,y);
# Find:  a nonzero linear difference operator L in QQ(x)[Sx]
#        such that
#                   L(T) = (Sy-1)(G),
#        where G is a hypergeometric term over QQ(x,y).
#
# Notation: F = QQ(x), Sx(x) = x+1, Sy(y) = y+1.
################################################################################

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

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

HypergeomTelescoping := proc(T,x,y,Sx)
  local fy,u,v,s,t,H,g,h,p,m,fx,d,i,N,gL,hL,pL,gg,hh,pp,cL,z,tp;

  #----------------------------------
  # trivial case: T is free of y
  #----------------------------------
  fy := ShiftQuotient(T,y);
  if fy = 1 then
    if nargs = 5 then assign(args[5]=[y,T]); fi;
    return(1);
  fi;
  
  #------------------------------------
  # compute an RNF of Sy(T)/T
  #------------------------------------
  u,v,s,t := op(RNF(fy,y));

  #-------------------------------------------------
  # compute the multiplicative decomposition of T
  #-------------------------------------------------
  H := normal(T/(s/t));

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

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

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

  #--------------------------------------------
  # compute the shift quotient of H w.r.t. x
  #--------------------------------------------
  N := normal(fx*(s/t)/subs(x=x+1,s/t));

  #----------------------------------------------------------
  # iteratively compute a minimal telescoper for T w.r.t. x
  #----------------------------------------------------------
  if nargs = 5 then gL := [RationalTools:-NormalizeRatSum(g)]; fi;
  hL := [h];
  pL := [p];
  i := 1;
  while i > 0 do
    gg,hh,pp := op(SimpleReduction(hL[i],pL[i],u,v,N,x,y));
    g,hh,p := op(LinearizeRemainders(hh,d,tp,u,v,y));
    if nargs = 5 then
      gL := [op(gL),RationalTools:-NormalizeRatSum(subs(x=x+1,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 = 5 then
	 g := RationalTools:-NormalizeRatSum(add(cL[j]*gL[j],j=1..m));
         assign(args[5]=[g,H]);
      fi;
      return(add(cL[j]*Sx^(j-1),j=1..m));
    fi;
  od;
end:

################################################################################
# Name: HypergeomReduction
# Calling sequence:
#                  HypergeomReduction(T,y)
# Input: T,  a hypergeometric term over F(y);
#        y,  a variable name.
# Output: [[g,r],H],  where
# 	    g,r - two rational functions in F(y),
#	    H   - a hypergeometric term over F(y) with Sy(H)/H = K
#	     	  being Sy-reduced,
#           such that
#           (1) T = (Sy-1)(g*H)+r*H;
#           (2) r is a Sy-remainder w.r.t. K.
################################################################################

HypergeomReduction := proc(T,y)
  local f,u,v,s,t,H,g,h,p,r,z;
  
  #------------------------------------
  # compute an RNF of Sy(T)/T
  #------------------------------------
  f := ShiftQuotient(T,y);
  u,v,s,t := op(RNF(f,y));

  #-------------------------------------
  # compute the hypergeometric term H
  #-------------------------------------
  H := normal(T/(s/t));

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

################################################################################
# Name: ShellReduction
# Calling sequence:
#                  ShellReduction(s,t,u,v,y)
# Input: s,t,u,v,         four polynomials in F[y] with
#	              	  (1) u/v being Sy-reduced;
#                     	  (2) gcd(s,t) = 1;
#	 y,          	  a variable name.
# 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 primitive, irreducible and strongly coprime with u/v;
#	    (4) the d_i are pairwise shift inequivalent;
#           (5) p is in the standard complement of the image of the map
#	    	for polynomial reduction w.r.t. u/v.
################################################################################

ShellReduction := proc(s,t,u,v,y)
  local a,p,g,h,b,f;

  #------------------------------------------------
  # compute the canonical representation of s/t
  #------------------------------------------------
  a := rem(s,t,y,'p');
  
  #------------------------------------------------
  # apply the normal reduction to the normal 
  # part of s/t w.r.t. u/v
  #------------------------------------------------
  g,h,b := op(NormalReduction(a,t,u,v,y));

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

################################################################################
# Name: PolynomialReduction
# Calling sequence:
#                  PolynomialReduction(b,u,v,y)
# Input: b,u,v,  three polynomials in F[y] with u/v being Sy-reduced;
#	 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.
################################################################################

PolynomialReduction := proc(b,u,v,y)
  local du,dv,lu,lv,a,p,d,i,g,pg,lp,cu,cv,m,bm,pm,r,gr,dr,lr,pr;
  
  #-----------------------
  # 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);
  a := 0;
  p := b;
    
  if du <> dv or lu <> lv then
    #-------------------------------------------------
    # case I: either degree(u,y) <> degree(v,y) 
    # 	    or lcoeff(u,y) <> lcoeff(v,y)
    #-------------------------------------------------
    d := max(du,dv);
    lu := coeff(u,y,d);
    lv := coeff(v,y,d);
    i := degree(p,y)-d;    
    while i >= 0 do
      g := y^i/(lu-lv);
      pg := RationalTools:-NormalizeRatSum(u*subs(y=y+1,g)-v*g);
      lp := lcoeff(p,y);
      a := a+lp*g;
      p := RationalTools:-NormalizeRatSum(p-lp*pg);
      i := degree(p,y)-d;
    od;
    return([a,p]);
  elif du = 0 then
    #---------------------------
    # case II: u/v = 1
    #---------------------------
    i := degree(p,y)+1;
    while i > 0 do
      g := y^i/(i*u);
      pg := normal(RationalTools:-NormalizeRatSum((y+1)^i-y^i)/i);
      lp := lcoeff(p,y);
      a := a+lp*g;
      p := RationalTools:-NormalizeRatSum(p-lp*pg);
      i := degree(p,y)+1;
    od;
    return([a,p]);
  fi;

  cu := coeff(u,y,du-1);
  cv := coeff(v,y,du-1);
  m := normal((cv-cu)/lu);
  if type(m,integer) and m >= 0 then
    #---------------------------------------------------------------------
    # case IV: degree(u,y) = degree(v,y) > 0, lcoeff(u,y) = lcoeff(v,y)
    #       and m*lu+coeff(u,y,du-1)-coeff(v,y,du-1) = 0
    #	    for some nonnegative integer m
    #---------------------------------------------------------------------
    bm := coeff(b,y,du+m-1)*y^(du+m-1);
    p := RationalTools:-NormalizeRatSum(p-bm);
    i := degree(p,y)-du+1;
    while i >= 0 do
      g := y^i/(i*lu+cu-cv);
      pg := RationalTools:-NormalizeRatSum(u*subs(y=y+1,g)-v*g);
      lp := lcoeff(p,y);
      a := a+lp*g;
      p := RationalTools:-NormalizeRatSum(p-lp*pg);
      pm := coeff(p,y,du+m-1)*y^(du+m-1);
      bm := bm+pm;
      p := RationalTools:-NormalizeRatSum(p-pm);
      i := degree(p,y)-du+1;
    od;
    
    r := RationalTools:-NormalizeRatSum(u*(y+1)^m-v*y^m);
    gr := y^m;
    i := degree(r,y)-du+1;
    while i >= 0 do
      g := lcoeff(r,y)*y^i/(i*lu+cu-cv);
      gr := gr-g;
      r := RationalTools:-NormalizeRatSum(r-u*subs(y=y+1,g)+v*g);
      i := degree(r,y)-du+1;
    od;

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

    pr := coeff(p,y,dr);
    a := a+pr*gr/lr;
    p := RationalTools:-NormalizeRatSum(p+bm-pr*r/lr);
  else
    #---------------------------------------------------------------------
    # case III: degree(u,y) = degree(v,y) > 0, lcoeff(u,y) = lcoeff(v,y)
    #       and i*lu+coeff(u,y,du-1)-coeff(v,y,du-1) <> 0
    #	    for all nonnegative integers i
    #---------------------------------------------------------------------
    i := degree(p,y)-du+1;
    while i >= 0 do
      g := y^i/(i*lu+cu-cv);
      pg := RationalTools:-NormalizeRatSum(u*subs(y=y+1,g)-v*g);
      lp := lcoeff(p,y);
      a := a+lp*g;
      p := RationalTools:-NormalizeRatSum(p-lp*pg);
      i := degree(p,y)-du+1;
    od;
  fi;
  
  [a,p];
end:

################################################################################
# Name: OrderBounds
# Calling sequence:
#                  OrderBounds(T,x,y,A)
# Input: T,    a hypergeometric term over QQ(x,y);
#        x,y,  two variable names;
#        A,    (optional) an argument name used for specifying the output.
# Output: [up,lo],  if T has a telescoper and the optional 4th 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 4th argument 
#		    'upper' is given.
#         lo,       if T has a telescoper and the optional 4th argument
#		    'lower' is given.
#         "No telescoper exists!",
#	            if T has no telescopers.
################################################################################

OrderBounds := proc(T,x,y)
  local fy,u,v,s,t,du,dv,lu,lv,tap,h,p,m,tab,i,d,k,tp,po,PL,l,up;

  #------------------------------------
  # compute an RNF of Sy(T)/T
  #------------------------------------
  fy := ShiftQuotient(T,y);
  u,v,s,t := op(RNF(fy,y));
  du := degree(u,y);
  dv := degree(v,y);
  lu := lcoeff(u,y);
  lv := lcoeff(v,y);
  tap := `if`(du>0 and du=dv and lu=lv,1,0);
  
  #------------------------------------------------------
  # apply the shell reduction to s/t w.r.t. u/v
  #------------------------------------------------------
  h,p := op(ShellReduction(s,t,u,v,y)[2..-1]);

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

  #----------------------------------------------------
  # check if the denominator of h is integer linear
  #----------------------------------------------------
  tab = table([]);
  for i from 1 to m do
    d,k := op(2..-1,h[i]);
    if not IsIrredIntegerLinear(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 IsShiftEquivalent(po,PL[i,1],y,'l') do
        i := i+1;
      od;
      if i <= nops(PL) then
        PL[i,2] := [op(PL[i,2]),[k,l]];
      else
        tab[tp] := [op(PL),[po,[[k,0]]]];
      fi;
    else
      tab[tp] := [[po,[[k,0]]]];
    fi;
  od;
  
  #----------------------------------
  # compute the order bounds
  #----------------------------------
  if nargs = 4 and type(args[4],'identical(upper)') then
    up := max(du,dv)-tap;
    return(up+UpperBound(tab,y));
  elif nargs = 4 and type(args[4],'identical(lower)') then
    return(LowerBound(tab));
  fi;

  up := max(du,dv)-tap;

  [up+UpperBound(tab,y),LowerBound(tab)];
end:

################################################################################
# Name: ShiftQuotient
# Calling sequence:
#                  ShiftQuotient(T,y)
# Input: T,  a term over F(y);
#        y,  a variable name.
# Output: f,  a rational function in F(y) such that f = Sy(T)/T,
#	      if T is a hypergeometric term over F(y);
#	  "Input is not a hypergeometric term in y",
#	      otherwise.
################################################################################

ShiftQuotient := proc(T,y)
  local f;

  if not SumTools:-Hypergeometric:-IsHypergeometricTerm(T,y,'f') then
   error "Input is not a hypergeometric term in %1",y;
  fi:
  
  f;
end:

################################################################################
# Name: RNF
# Calling sequence:
#                  RNF(f,y)
# Input: f,           a nonzero rational function in F(y);
#        y,           a variable name.
# 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.
################################################################################

RNF := proc(f,y)
  local ff,c,nf,df,u,v,s,t;

  #-------------------------------
  # trivial case: f is free of y
  #-------------------------------
  ff := normal(f);
  if type(ff,freeof(y)) then
    return([ff,1,1,1]);
  fi;
  
  #--------------------------------------------
  # compute an RNF of the primitive part
  #--------------------------------------------
  c := normal(content(numer(ff),y,'nf')/content(denom(ff),y,'df'));
  u,v,s,t := op(PrimitiveRNF(nf,df,y));

  #--------------------------
  # standardize the RNF
  #--------------------------
  u := normal(c*u);
  
  [u,v,s,t];
end:

################################################################################
# Name: DispersionSet
# Calling sequence:
#                  DispersionSet(a,b,y)
# Input: a,b,  two polynomials in F[y];
#        y,    a variable name.
# Output: DS,  a finite set of all integers i such that gcd(a,Sy^i(b)) <> 1.
################################################################################

DispersionSet := proc(a,b,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 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 IsShiftEquivalent(afac[i],bfac[j],y,'m') then
        DS := DS union {m};
      fi;
    od;
  od;

  DS;
end:

################################################################################
# Name: ShiftFactorization
# Calling sequence:
#                  ShiftFactorization(d,y,[u,v]);
# Input: d,      a polynomial in F[y] with degree(d,y) > 0;
#	 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 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 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 3th 
#	    	argument [u,v] is given.
################################################################################

ShiftFactorization := proc(d,y)
  local pd,cd,c,dfac,S,n,p,k,kl,temp,i,p1,k1,l,u,v,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 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 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 IsShiftEquivalent(p1,p,y,'l') then
        kl := [op(kl),[k1,l]];
      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 shift
  # equivalence class strongly coprime with u/v
  # if the optional 3th argument [u,v] is given
  #--------------------------------------------------------
  if nargs = 3 then
    u,v := op(args[3]);
    for i from 1 to nops(S) do
      p,kl := op(S[i]);
      l := map(z->z[2],kl);
      udset := DispersionSet(u,p,y);
      if udset = {} then
        vdset := DispersionSet(v,p,y);
	m := min(min(l),min(vdset)-1);
      else
        m := max(max(l),max(udset)+1);
      fi;
      p := subs(y=y+m,p);
      S[i] := [p,map(z->[z[1],z[2]-m],kl)];
    od;
  fi;
  
  [c,S];
end:

######## Verification

################################################################################
# Name: VerifyHypergeomTelescoping
# Calling sequence:
#                  VerifyHypergeomTelescoping([L,[g,H]],T,x,y,Sx)
# Input: [L,[g,H]],  a list with L a nonzero linear difference operator
# 	 	     in QQ(x)[Sx], g in QQ(x,y) and H a hypergeometric
#		     term over QQ(x,y);
#	 T,          a hypergeometric term over QQ(x,y);
#        x,y,        two variable names;
#        Sx,         an operator name used for the shift operator w.r.t. x.
# Output: true,   if the following conditions are satisfied:
# 	  	  (1) L(T) = (Sy-1)(g*H);
#		  (2) the Sy-quotient of H is Sy-reduced;
#         false,  otherwise.
################################################################################

VerifyHypergeomTelescoping := proc(res,T,x,y,Sx);
  local L,cert,g,H,S,K,N,r,i,j;
  
  #----------------------------------------------
  # compute the kernel and the shell of Sy(T)/T
  #----------------------------------------------
  L,cert := op(res);
  g,H := op(cert);
  S := normal(T/H);
  K := ShiftQuotient(H,y);

  #----------------------------------------------
  # compute the shell of Sy(L(T))/L(T) w.r.t. K
  #----------------------------------------------
  N := ShiftQuotient(H,x);
  r := add(coeff(L,Sx,i)*subs(x=x+i,S)*mul(subs(x=x+j,N),j=0..i-1),i=0..degree(L,Sx));

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

  #-------------------------------
  # check the second condition
  #-------------------------------
  if nops(DispersionSet(numer(K),denom(K),y)) <> 0 then
    print(The_2nd_condition_is_violated);
    return(false);
  fi;
  
  true;
end:

################################################################################
# Name: VerifyHypergeomReduction
# Calling sequence:
#                  VerifyHypergeomReduction([[g,r],H],T,y)
# Input: [[g,r],H],  a list with g,r in F(y) and H a hypergeometric 
# 	 	     term over F(y);
#        T,  	     a hypergeometric term over F(y);
#        y,  	     a variable name.
# Output: true,   if the following conditions are satisfied:
# 	    	  (1) T = (Sy-1)(g*H)+r*H;
#           	  (2) the Sy-quotient K of H is Sy-reduced;
#		  (3) r is a Sy-remainder w.r.t. K;
#         false,  otherwise.
################################################################################

VerifyHypergeomReduction := proc(res,T,y)
  local gr,H,S,K,g,r,u,v,nr,dr,b,cv,a,p,du,dv,lu,lv,d,cu,m,i;
  
  #----------------------------------------------
  # compute the kernel and the shell of Sy(T)/T
  #----------------------------------------------
  gr,H := op(res);
  S := normal(T/H);
  K := ShiftQuotient(H,y);
  g,r := op(gr);

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

  #-------------------------------
  # check the second condition
  #-------------------------------
  u,v := numer(K),denom(K);
  if nops(DispersionSet(u,v,y)) <> 0 then
    print(The_2nd_condition_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 DispersionSet(b,b,y) <> {0} then 
    print(The_3rd_condition_2_is_violated);
    return(false);
  fi;

  if max(DispersionSet(u,b,y)) >= 0 then
    print(The_3rd_condition_3_is_violated);
    return(false);
  fi;

  if min(DispersionSet(v,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 lu <> lv 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 p <> 0 then
      print(The_3rd_condition_5_is_violated);
      return(false);
    fi;
  else
    cu := coeff(u,y,du-1);
    cv := coeff(v,y,du-1);
    m := normal((cv-cu)/lu);
    if type(m,integer) and m >= 0 then
      r := RationalTools:-NormalizeRatSum(u*(y+1)^m-v*y^m);
      i := degree(r,y)-du+1;
      while i >= 0 do
        g := lcoeff(r,y)*y^i/(i*lu+cu-cv);
        r := RationalTools:-NormalizeRatSum(r-u*subs(y=y+1,g)+v*g);
        i := degree(r,y)-du+1;
      od;
      dr := i+du-1;
      p := RationalTools:-NormalizeRatSum(p-coeff(p,y,du+m-1)*y^(du+m-1));
      if degree(p,y) >= du-1 or coeff(p,y,dr) <> 0 then
        print(The_5th_condition_is_violated);
        return(false);
      fi;
    else
      if degree(p,y) >= du-1 then
        print(The_5th_condition_is_violated);
        return(false);
      fi;
    fi;
  fi;
  
  true;
end:

################################################################################
# Name: VerifyShellReduction
# Calling sequence:
#                  VerifyShellReduction([g,h,p],s,t,u,v,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-reduced;
#                  (2) gcd(s,t) = 1;
#	 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 primitive, irreducible and strongly coprime with u/v;
#	    	 (4) the d_i are pairwise 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.
################################################################################

VerifyShellReduction := proc(res,s,t,u,v,y)
  local g,h,p,r,i,a,d,k,b,du,dv,lu,lv,cu,cv,m,dr,z;
  
  #-------------------------------
  # check the first condition
  #-------------------------------
  g,h,p := op(res);
  r := RationalTools:-NormalizeRatSum(u/v*subs(y=y+1,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 the 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 not irreduc(d) then
      print(The_3rd_condition_1_is_violated);
      return(false);
    fi;

    if max(DispersionSet(u,d,y)) >= 0 then
      print(The_3rd_condition_2_is_violated);
      return(false);
    fi;
    
    if min(DispersionSet(v,d,y)) <= 0 then
      print(The_3rd_condition_3_is_violated);
      return(false);
    fi;
  od;

  #-------------------------------
  # check the fourth condition
  #-------------------------------
  b := mul(map(z->z[2],h));
  if nops(DispersionSet(b,b,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 lu <> lv 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 p <> 0 then
      print(The_3rd_condition_5_is_violated);
      return(false);
    fi;
  else
    cu := coeff(u,y,du-1);
    cv := coeff(v,y,du-1);
    m := normal((cv-cu)/lu);
    if type(m,integer) and m >= 0 then
      r := RationalTools:-NormalizeRatSum(u*(y+1)^m-v*y^m);
      i := degree(r,y)-du+1;
      while i >= 0 do
        g := lcoeff(r,y)*y^i/(i*lu+cu-cv);
        r := RationalTools:-NormalizeRatSum(r-u*subs(y=y+1,g)+v*g);
        i := degree(r,y)-du+1;
      od;
      dr := i+du-1;
      p := RationalTools:-NormalizeRatSum(p-coeff(p,y,du+m-1)*y^(du+m-1));
      if degree(p,y) >= du-1 or coeff(p,y,dr) <> 0 then
        print(The_5th_condition_is_violated);
        return(false);
      fi;
    else
      if degree(p,y) >= du-1 then
        print(The_5th_condition_is_violated);
        return(false);
      fi;
    fi;
  fi;
  
  true;
end:

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

################################################################################
# Name: NormalReduction
# Calling sequence:
#                  NormalReduction(s,t,u,v,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) and gcd(s,t) = 1;
#	 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 primitive, irreducible and strongly coprime with u/v;
#	    (4) the d_i are pairwise shift inequivalent;.
################################################################################

NormalReduction := proc(s,t,u,v,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(ShiftFactorization(t,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=y+z[2],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(LocalNormalReduction(aL[j],d,kk,l,u,v,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: LocalNormalReduction
# Calling sequence:
#                  LocalNormalReduction(c,d,k,l,u,v,y)
# Input: c,d,  two polynomials in F[y] with
# 	       (1) 0 <= degree(c,y) < degree(d^k,y);
#	       (2) d being 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;
#	 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.
################################################################################

LocalNormalReduction := proc(c,d,k,l,u,v,y)
  local g,kk,a,b,dd,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;
  a := c;
  b := 0;

  if l > 0 then
    dd := d;
    di := subs(y=y+l,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=y-1,s);
	vi := normal(vi*di^(kk-ki));
	kk := ki;
	di := subs(y=y+i-1,dd);
	g := g+a/di^kk;
	b := RationalTools:-NormalizeRatSum(b+t);
      fi;
    od;
  else
    dd := subs(y=y+l,d);
    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=y+1,di);
      while ki > 0 and RationalTools:-IsDivided(pu,di,y,'pu') do
       ki := ki-1;
      od;
      a := normal(cu*pu*ui*subs(y=y+1,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;
  fi;

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

################################################################################
# Name: SimpleReduction
# Calling sequence:
#                  SimpleReduction(h,p,u,v,N,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 primitive, irreducible and strongly coprime
#		     with u/v;
#        p,u,v,  three polynomials in F[y] with
#	 	 (1) u/v being Sy-reduced;
#		 (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);
#	 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 primitive, irreducible and strongly coprime with u/v;
#	    (4) the dd_i are pairwise shift inequivalent;
#           (5) pp is in the standard complement of the image of the map
#	    	for polynomial reduction w.r.t. u/v.
################################################################################

SimpleReduction := proc(h,p,u,v,N,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=x+1,p/v)*NN);
  g,h0,pp := op(ShellReduction(numer(r),denom(r),u,v,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=x+1,a/d^k)*NN);
    gi,hi,pi := op(ShellReduction(numer(r),denom(r),u,v,y));
    g := g+gi;
    hh[i] := op(hi);
    pp := RationalTools:-NormalizeRatSum(pp+pi);
  od;
  hh := convert(hh,list);

  [g,hh,pp];
end:

################################################################################
# Name: LinearizeRemainders
# Calling sequence:
#                  LinearizeRemainders(h,d,tp,u,v,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 primitive, irreducible and strongly coprime
#	       	   with u/v;
#        d,    a list of the form [d0_1,...,d0_m] with d0_i being 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 integer linear of type (lambda_i,mu_i);
#        u,v,  two polynomials in F[y] with u/v being Sy-reduced;
#	 y,    a variable name.
# Output: [g,hh,p],  where
# 	   g  - a rational function in QQ(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 primitive, irreducible, strongly coprime with u/v;
#	   (4) the dd_i are pairwise 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.
################################################################################

LinearizeRemainders := proc(h,d,tp,u,v,y)
  local g,p,hh,n,id,m,i,a,dd,k,lm,j,l,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 IsShiftEquivalent(dd,d[id[j]],y,'l'))\
    do
      j := j+1;
    od;
    if j <= m then
      d0 := d[id[j]];
      gg,a,k,b := op(LocalNormalReduction(a,d0,k,l,u,v,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(PolynomialReduction(p,u,v,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 primitive and irreducible;
#			 (4) d_{i1},...,d_{in} are pairwise shift inequivalent;
# 	 [p_1,...,p_m],  a list of polynomials 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: UpperBound
# Calling sequence:
#                  UpperBound(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[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*degree(P_i,y),i=1..m)
#		for mu running through the indices of tab.
################################################################################

UpperBound := 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*max(map(zz->zz[1],z[2]))*degree(z[1],y),tab[[lambda,mu]])));
  od;

  ord;
end:

################################################################################
# Name: LowerBound
# Calling sequence:
#                  LowerBound(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[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.
################################################################################

LowerBound := 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:

################################################################################
# Name: PrimitiveRNF
# Calling sequence:
#                  PrimitiveRNF(nf,df,y)
# Input: nf,df,  two coprime primitive polynomials in F[y];
#        y,  	 a variable name.
# Output: [u,v,s,t],  where
#           u,v,s,t - four 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.
################################################################################

PrimitiveRNF := proc(nf,df,y)
  local u,v,s,t,DS,pDS,nDS,i,k,p,pk,j,st;

  #-----------------------------------------
  # trivial case: nf or df is free of y
  #-----------------------------------------
  u := nf;
  v := df;
  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 := DispersionSet(u,v,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=y+k,v),y,'u');
    pk := subs(y=y-k,p);
    v := normal(v/pk);
    s := s*mul(subs(y=y-j,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=y+k,v),y,'u');
    pk := subs(y=y-k,p);
    v := normal(v/pk);
    t := t*mul(subs(y=y+j,p),j=0..-k-1);
  od;
  st := normal(s/t);
  
  [u,v,numer(st),denom(st)];
end:

################################################################################
# Name: IsIrredIntegerLinear
# Calling sequence:
#                  IsIrredIntegerLinear(p,x,y,'tp','po')
# Input: p,      an irreducible polynomial in QQ[x,y] with degree(p,y) > 0;
# 	 x,y,    two variable names;
#	 tp,po,  (optional) two unevaluated names.
# Output: true,   if p is 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[y] such that
#		                   p = P(lambda*x+mu*y);
# 	  false,  otherwise.
################################################################################

IsIrredIntegerLinear := proc(p,x,y)
  local d,hL,z_,tL,dL,n,i,P,m,h,k,c,sq,f,l,lambda,mu,cL,id,z,zz;

  #------------------------------------
  # trivial case: p is free of x
  #------------------------------------
  d := degree(p,x);
  if d = 0 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;
  
  #---------------------------------------------
  # check whether degrees in x and y are equal
  #---------------------------------------------
  if d <> degree(p,y) then return false; fi;

  #---------------------------------
  # compute homogeneous components
  #---------------------------------
  hL := [coeffs(collect(subs(x=z_*x,y=z_*y,p),z_),z_,'tL')];
  dL := map(degree,[tL],z_);
  n := nops(dL);
  hL := sort([seq([hL[i],dL[i]],i=1..n)],(z,zz)->z[2]<zz[2]);
  
  #-----------------------------------------------------------------
  # find the possible integer linear type
  # Note that p has at least one nontrivial homogeneous component 
  # since p is irreducible with degree(p,x) = degree(p,y) > 0.
  #------------------------------------------------------------------
  if hL[1][2] = 0 then
    P := hL[1][1];
    m := 2;
  else
    P := 0;
    m := 1;
  fi;
  h,k := op(hL[m]);
  c,sq := op(sqrfree(h,[x,y]));
  f,l := op(sq[1]);
  lambda,mu := op(map2(coeff,f,[x,y],1));

  if nops(sq) = 1 and k = l and type(lambda,integer) and lambda <> 0 and\
    type(mu,integer) and mu <> 0 then
    if mu < 0 then lambda,mu := -lambda,-mu; fi;
    P := P + c*y^k;
  else
    return false;
  fi;

  #----------------------------------------------
  # check whether each homogeneous component is
  # a power of (lambda*x+mu*y)
  #----------------------------------------------
  cL := table([]);
  for i from m+1 to n do
    h,k := op(hL[i]);
    if not RationalTools:-IsDivided(h,(lambda*x+mu*y)^k,y,'cL[k]') then
      return false;
    fi;
  od;
  
  if nargs = 4 then assign(args[4]=[lambda,mu]); fi;
  if nargs = 5 then
    id := [indices(cL,'nolist')];
    for i from 1 to nops(id) do
      k := id[i];
      P := P + cL[k]*y^k;
    od;
    assign(args[4]=[lambda,mu],args[5]=P);
  fi;
  
  true;
end:

################################################################################
# Name: IsShiftEquivalent
# Calling sequence:
#                  IsShiftEquivalent(a,b,y,'m');
# Input: a,b,  two primitive polynomials in F[y] of positive degrees in y;
#        y,    a variable name;
#	 m,    (optional) an unevaluated name.
# Output: true,   if a = Sy^k(b) for some integer k, and then the optional
# 	  	  4th argument m is assigned the integer k.
#         false,  otherwise.
################################################################################

IsShiftEquivalent := proc(a,b,y)
  local d,la,m;

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

  #--------------------------------
  # check the leading coefficient
  #--------------------------------
  la := lcoeff(a,y);
  if la <> lcoeff(b,y) then return false; fi;
  
  #-------------------------------
  # find and check the candidate
  #-------------------------------
  m := normal((coeff(a,y,d-1)-coeff(b,y,d-1))/d/la);
  if type(m,integer) and RationalTools:-NormalizeRatSum(a-subs(y=y+m,b)) = 0 then
      if nargs = 4 then assign(args[4]=m); 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:
