#################################################################################
# Implementation of Zeilberger's algorithm                                      #
# Uses an algorithm constructing minimal annihilators in case of homogeneous    #
# Zeilberger recurrences.                                                       #
#                                                                               #
# Reference:                                                         		#
#   M. Petkovsek, H. Wilf, D. Zeilberger (1996):                     		#
#   A=B, A.K.Peters, Wellesley, Massachusetts. (Chapter 6)           		#
#   S. P. Polyakov (2008):                                                      #
#   On homogeneous Zeilberger recurrences, Advances in Applied Mathematics,     #
#   vol. 40, N 1, pp. 1-7                                                       #
#                                                                               #
#                                                        Ha Le, 2001-2004       #
#                                                Stanislav Polyakov, 2007       #
#                                                                               #
#################################################################################

   Zeilberger := proc(FF,n,k,N)
     local F, s, r, r1, r2, s1, s2, MaxOrder, MinOrder, J, p0, i, j,
           r_k, s_k, p1, p2, p3, b, L, G, p, vars, t_k, is_hyp, msg,
           `z/coefs`, cns, inS1, Z, d, P, EnHoH, A, Lc, p_i, ab, rm, R;
     description "Zeilberger's algorithm";

     # Type checking
     F := FF;
     is_hyp := IsHypergeometricTerm(F,n,'s','down');
     if not is_hyp then
        error "The input is expected to be a hypergeometric term in %1 and %2",
               n,k;
     end if;
     is_hyp := IsHypergeometricTerm(F,k,'r','up');
     if not is_hyp then
        error "The input is expected to be a hypergeometric term in %1 and %2",
               n,k;
     end if;

     # check if F is an element of S_0
     if type(s,freeof(k)) then
        try
           G := Gosper(F,k);
        catch:
           userinfo(3,{'SumTools','Zeilberger'},
              "the input is an element of S_0");
	   L := TelescoperS_0(eval(s,n=n+1),N);
	   return [L,1]
	end try;
        return [1,G]
     end if;

     Z := true;
     # check if F is an element of S
     P := eval(s,n=n+1);
     inS1 := IsInS1(P,n,k,'p');
     if inS1 then
        ab:=annihilbasis(p,n,k);
        rm:=nops(ab);
        R := radnormal( P*p/eval(p,n=n+1) );
        d := ComputeUpperBound(r,k);
print(d);
	if d < 0 then Z := false end if;
     else rm:=infinity; end if;

     if Z then

       # specify the range or the order of the difference operator L    
       if assigned(_MAXORDER) and
          (type(_MAXORDER,'nonnegint') or _MAXORDER = infinity) then
          MaxOrder := _MAXORDER
       else
          MaxOrder := 6
       end if;
       if assigned(_MINORDER) and type(_MINORDER,nonnegint) then
          MinOrder := _MINORDER
       else
          MinOrder := 0
       end if;

       s1 := numer(s); s2 := denom(s);
       r1 := numer(r); r2 := denom(r);

       for J from MinOrder to MaxOrder do
          if J>=rm then
             L:=annihilator(ab,n,N,R);
             return [L,0];
          end if; 
          msg := sprintf(`applying Zeilberger's algorithm for order %a`,J);
          userinfo(3,{'SumTools','Zeilberger'},msg);
          vars := [seq(`z/coefs`[J-j],j=0..J)]; j := 'j';
          p0 := add(`z/coefs`[j]*(mul(eval(s1,n=n+j-i),i=0..j-1)*
                    mul(eval(s2,n=n+r),r=j+1..J)),j=0..J);
          r := 'r'; r_k := r1*mul(eval(s2,n=n+r),r=1..J);
          r := 'r'; s_k := r2*mul(eval(eval(s2,n=n+r),k=k+1),r=1..J);
          (cns,p2,p3,p1) := PolynomialNormalForm(Normalizer(r_k/s_k),k);
          p2 := cns*p2;
          p := p0*p1;
          try
            (p2,p3,p) := reduce(p2,p3,p,k);
            catch:
          end try;
          p := collect(expand(p),k);

          try
            b := GosperStep3(p2,p3,p,k,vars);
            catch "no nonzero polynomial solution":
            next;
          end try;
          vars := b[2]; b := b[1];
          p := eval(p,vars);
          L := eval(`+`(seq(`z/coefs`[j]*N^j,j=0..J)),vars);
          if L = 0 then next end if;
          L := collect(primpart(L,N),N);
          t_k := ApplyLtoF(L,F,N,n);                                    
          # This is the case when L is also the minimal annihilator for
          # F(n,k).
          if t_k = 0 or b = 0 then G := 1
          else G := eval(p3,k=k-1)/p*b*t_k end if;
          userinfo(3,{'SumTools','Zeilberger'},
                   "Zeilberger \'s algorithm successful");
          return [L,G];
        end do;
        error "No recurrence of order %1 was found", MaxOrder;

     else
         L:=annihilator(ab,n,N,R);
         return [L,0];
     end if;

   end proc:

   IsInS1 := proc(P,n,k,a)
      local nu, de, g, g2, req, psol;

      nu := collect(numer(P),k); de := collect(denom(P),k);
      g := primpart(de,k); g2 := primpart(nu,k);
      if degree(g,n) <> degree(g2,n) then
         false
      else
         req := g*z(n+1) = g2*z(n);
	 psol := LREtools['polysols'](req,z(n),{},'output'='onesol');
	 if psol = NULL or psol = 0 then
	    false
	 else
	    if nargs=4 and type(args[4],'name') then
	       a := numer(normal(psol))
	    end if;
	    true
	 end if;
      end if;
   end proc:

   ComputeUpperBound := proc(P,k)
      local co,a,b,c;

      (co,a,b,c) := PolynomialNormalForm(P,k);
      a := co*a;
      UpperBound(a,b,c,k)

   end proc:

   TelescoperS_0 := proc(s,N)
      denom(s)*N - numer(s)
   end proc:

   # attempt to reduce the size of the linear system of equations
   # to be solved in step 3 of Gosper's algorithm
   reduce := proc(p2,p3,p,n)
     local a, b, c, g, g1;
     description "attempt to reduce the size of the linear system";

     a := p2; b := p3; c := p;
     g := gcd(a,c,'a','c');
     if g <> 1 then
       g1 := eval(g,n=n+1);
       gcd(b,g1,'b','g1');
       a := Normalizer(a*g1);
     end if;
     (a,b,c)
   end proc:

#################################################################################
# Construct the polynomial normal form (or Gosper-Petkovsek form)		#
#										#
# Reference: A=B, M. Petkovsek, H. Wilf, D. Zeilberger (1996),			#
#            A.K.Peters, Wellesley, Massachusetts.				#
#            Gosper's Algorithm (Step 2), page 80, Section 5.3			#
#										#
#                                       Ha Le, October 2000			#
#										#
#################################################################################

   PolynomialNormalForm := proc(F,n)
     local f,g,z1,z2,z,H,N,p,q,s,j,i,prod,param,pp,qq,monic_a,monic_b,coefa,coefb;
     description "construct the polynomial normal form";             

     userinfo(3,{'SumTools','PolynomialNormalForm'},
              "construct the polynomial normal form");

     if not type(F,ratpoly(anything,n)) then
       error "wrong type of argument"
     end if;

     # Let F = z * f/g where f, g are monic relatively
     # prime polynomials, and z is a constant.
     # Compute H = [h_1,...,h_n] which is a list of
     # nonnegative integers such that
     # gcd(f,g+h) <> 1 for every element in H.
     f := numer(F); g := denom(F);
     gcd(f,g,'f','g');
     z1 := lcoeff(f,n); z2 := lcoeff(g,n);
     z := z1/z2; param := true;
     if not hastype(z1,'name') and not hastype(z2,'name') then
       f := f/z1; g := g/z2;
       param := false;
     end if;
     if degree(f,n)=0 or degree(g,n)=0 then H := FAIL
     else
        H := LREtools['dispersion'](g,f,n);
     end if;
     if H = FAIL then # [a,b,c]                               
       if not param then
         return z,MakeMonic(f,n)[2],MakeMonic(g,n)[2],1
       else
         return z,f/z1,g/z2,1
       end if;
     end if;
     H := sort(convert(H,list));

     # Compute the PolynomialNormalForm form for the input rational
     # function F
     N := nops(H);
     p := array(0..N); q := array(0..N); s := array(0..N);
     p[0] := f; q[0] := g;
     for j to N do
       s[j] := gcd(p[j-1],eval(q[j-1],n=n+H[j]),'coefa','coefb');
       p[j] := coefa; q[j] := eval(coefb,n=n-H[j]);
     end do;
     i := 'i'; j := 'j';
     prod := mul( mul( eval(s[i],n=n-j), j=1..H[i]), i=1..N);
     pp := MakeMonic(p[N],n); qq := MakeMonic(q[N],n);                   
     if not param then
       (Normalizer(pp[1]/qq[1]*z),pp[2],qq[2],MakeMonic(prod,n)[2])
     else
       monic_a := MakeMonic(pp[2],n);
       monic_b := MakeMonic(qq[2],n);
       (Normalizer(pp[1]/qq[1]*monic_a[1]/monic_b[1]),
        monic_a[2],monic_b[2],MakeMonic(prod,n)[2])
     end if;
   end proc:

#################################################################################
#                                                                               #
# Submodule: Hypergeometric                                                     #
# Input: a hypergeometric term t w.r.t. n                            		#
# Output: a hypergeometric term, say z w.r.t. n, such that           		#
#                z(n+1) - z(n) = t(n), if one exists;                		#
#         otherwise, it returns the error message                    		#
#         "no solution found".                                       		#
# Reference:                                                         		#
#   M. Petkovsek, H. Wilf, D. Zeilberger (1996):                     		#
#   A=B, A.K.Peters, Wellesley, Massachusetts. (Section 5.2, p77)    		#
#                                                                    		#
#					Ha Le, November 2001			#
#										#
#################################################################################

   Gosper := proc(t,n::name)
     local a, b, c, r, x, chk, z;
     description "Gosper's algorithm";

     # step 1 of Gosper's algorithm
     # form the ratio r(n) = t(n+1)/t(n) which should be a rational
     # function of n; otherwise, return error message.
     userinfo(3,{'SumTools','Gosper'},"Step 1 of Gosper's algorithm");
     chk := IsHypergeometricTerm(t,n,'r');
     if not chk then
       error "the input is expected to be a hypergeometric term in %1",n;
     end if;

     # step 2 of Gosper's algorithm
     # write r(n) = a(n)/b(n) * c(n+1)/c(n) where a(n),b(n),c(n)
     # are polynomials satisfying
     #     gcd(a(n),b(n+h)) = 1, for all nonnegative integers h.
     userinfo(3,{'SumTools','Gosper'},"Step 2 of Gosper's algorithm");
     (z,a,b,c) := PolynomialNormalForm(r,n);
     a := z*a;

     # step 3 of Gosper's algorithm
     # find a nonzero polynomial solution x(n) of
     # a(n)*x(n+1) - b(n-1)*x(n) = c(n), if one exists;               
     # otherwise, return error message.
     userinfo(3,{'SumTools','Gosper'},"Step 3 of Gosper's algorithm");
     try
       x := GosperStep3(a,b,c,n,[]);
       catch:
         error "no solution found";
     end try;
     Normalizer(eval(b,n=n-1)*x/c)*t;
   end proc:

   UpperBound := proc(a,b,c,n)
     local a1, b1, c1, deg_a, deg_b, deg_c, lc_a, lc_b, DD, A, B;

     a1 := collect(a,n); b1 := collect(b,n); c1 := collect(c,n);      
    
     deg_a := degree(a1,n); deg_b := degree(b1,n); deg_c := degree(c1,n);
     lc_a := Normalizer(lcoeff(a1,n)); lc_b := Normalizer(lcoeff(b1,n));
     if deg_a <> deg_b or lc_a <> lc_b then
       DD := {deg_c - max(deg_a,deg_b)};
     else
       A := coeff(a1,n,deg_a-1); B := coeff(eval(b,n=n-1),n,deg_b-1);
       DD := {deg_c-deg_a+1,Normalizer((B-A)/lc_a)};
     end if;
     DD := select(type,DD,nonnegint);
     if DD = {} then -1
     else max(op(DD)) end if;

   end proc:


   # Step 3 of Gosper's algorithm.
   # Reference: A=B, M. Petkovsek, H. Wilf, D. Zeilberger (1996),
   #            A.K.Peters, Wellesley, Massachusetts.
   #            Gosper's Algorithm (Step 3), page 86, Section 5.4
   #
   GosperStep3 := proc(a,b,c,n,vars1)
     local d, x, sys, vars, sols, i, G, msg, `g/coes`;
     description "a variant of step 3 of Gosper's algorithm";

     userinfo(5,{'SumTools','Gosper'},"find non-zero polynomial solution");
     d := UpperBound(a,b,c,n);
     userinfo(5,{'SumTools','Gosper'},"upper bound", d);

     if d<0 then x:=0 else x := `+`(seq(`g/coes`[i]*n^i,i=0..d)) end if;
     sys :=
       {coeffs(collect(expand(a*eval(x,n=n+1)-eval(b,n=n-1)*x-c),n),n)};
     vars := [seq(`g/coes`[d-i],i=0..d),
                op(vars1)];

     msg := sprintf(`size of the system: %a equations, %a unknowns`,
            nops(sys),nops(vars));
     userinfo(5,{'SumTools','Gosper'},msg);

     sols := [SolveTools:-Linear(sys,{op(vars)})];

     if sols = [] then
       error "no nonzero polynomial solution"
     end if;
     vars := MinimalLinearSpecialization([op(op(sols))], vars);
     sols := map((x,y)->lhs(x)=eval(rhs(x),y),op(sols),{op(vars)});
     G := eval(x,sols);

     if nops(vars1) = 0 then
       G
     else
       [G,sols]
     end if;
   end proc:


   ########################################################################
   #Calling sequence
   #
   #      MinimalLinearSpecialization(S, Par)
   #
   # Input: S: a list of nontrivial solutions of SolveTools:-Linear 
   #        Par: a list of unknowns ordered decreasingly
   #
   # Output: T: a nonzero solution obtained by specialzing parameters in S.
   #            (try to make nonzero entries appear as late as possible).
   #
   #						Ziming Li
   #
   #########################################################################
   
   MinimalLinearSpecialization := proc(S, Par)
      local l, A, i, f, p, le, re, S1, par, a ,b, c, vars, T;
   
      #########################
      # initialize
      #########################
   
      l := nops(Par);
   
      ##########################
      # special case
      ##########################
   
      if l = 0 then return [] end if;
   
      
      ##########################
      # sort solutions
      ##########################
      
      A := Array(1..l);
      for i from 1 to l do 
          f := S[i];
          if member(lhs(f), Par, 'p') then 
             A[p] := f;
          end if;
      end do;
      
      le := lhs(A[1]); re := rhs(A[1]);
      vars := indets(rhs(A[1])) intersect {op(Par)};
   
      ###############################
      # CASE1 no specialization
      ###############################
   
      if nops(vars) = 0 then
         if re <> 0 then 
            par := map(x->x=0, {op(Par)});
            return [A[1], seq(lhs(A[i])=eval(rhs(A[i]), par), i=2..l)];
         else
            T := MinimalLinearSpecialization([seq(A[i], i=2..l)], 
                                              [op(2..l, Par)]);  
            return [A[1], op(1..l-1, T)];
         end if;
      end if;
   
      #################################
      # CASE2 specialization
      ################################
   
      a := op(1, vars); 
      b := eval(rhs(A[1]),a=0);
      c := -b/lcoeff(rhs(A[1]), a);
      S1 := [seq(lhs(A[i])=evala(Simplify(eval(rhs(A[i]), a=c))), i=2..l)];
   
      if nops(vars) = 1 then
         if l = 1 then 
            if b <> 0 then 
               return [ lhs(A[1]) = b];
            else
               return [ lhs(A[1]) = eval(rhs(A[1]), a=1)];
            end if;
         else
            if b <> 0 or {seq(rhs(S1[i]), i=1..l-1)} <> {0} then
               T := MinimalLinearSpecialization(S1, [op(2..l, Par)]);
               return [le=0, op(1..l-1, T)];
            else
               return [le=eval(re, a=1), 
                       seq(lhs(A[i])=evala(
                                 Simplify(eval(rhs(A[i]), a=1))), i=2..l)];
            end if;
         end if;
       else
         T := MinimalLinearSpecialization(S1, [op(2..l, Par)]);
         [le=0, op(1..l-1, T)];
       end if;
   
   end proc: 
     
#################################################################################
# Name: MakeMonic (used internally)						#
#										#
# Make each factors of f \in K[x] monic						#
#										#
#                                       Ha Le, October 2000                     #
#										#
#################################################################################

MakeMonic := proc(f,x)
    local lc, newf, l, i, oi, oi1;

    if type(f,freeof(x)) then return(f,1) end if;
    newf := f; lc := 1;
    if type(f,`+`) then
      lc := lcoeff(newf,x); newf := collect(newf/lc,x,'normal');
    elif type(f,`^`) then
      oi := op(1,newf); l := lcoeff(oi,x);
      lc := l^op(2,newf); newf := collect(oi/l,x,'normal')^op(2,newf);
    elif type(f,`*`) then
      lc := 1; newf := 1;
      for i to nops(f) do
        oi := op(i,f);
        if type(oi,`^`) then
          oi1 := op(1,oi);
          l := lcoeff(oi1,x);
          lc := lc * l^op(2,oi);
          newf := newf * collect(oi1/l,x,normal)^op(2,oi);
        else
          l := lcoeff(oi,x);
          lc := lc * l;
          newf := newf * collect(oi/l,x,normal);
        end if;
      end do;
    end if;
    lc,newf
end proc:

#################################################################################
# Name: IsHypergeometricTerm							#
#										#
# Check to see if a given expression h is a hypergeometric term of a particular #
# variable n. If the third optional argument RR, which is any name is given,    #
# and if h is a hypergeometric term, then RR will be assigned to the 		#
# certificate of h.								#
#										#
#                                       Ha Le, October 2000			#
#										#
# Allowed the optional argument 'up' or 'down'. If it is 'up' then the		#
# certificate is defined as h(n+1)/h(n); otherwise, h(n)/h(n-1). It is useful	#
# for the implementation of Zeilberger's algorithm. Note that this option is	#
# used internally.								#
#										#
#                                       Ha Le, November 2001			#
#										#
#################################################################################

IsHypergeometricTerm := proc(h,n,RR)
  description                                                  
    "Check to see if the given expression is a hypergeometric term";
  local R, dir, rat, i, a, b, func, ind;

  dir := 'up'; rat := false;
  for i from 3 to nargs do
    if member(args[i],{'up','down'}) then
      dir := args[i]
    elif type(args[i],'name') then
      rat := true
    end if;
  end do;

  if dir = 'up' then
    # add a smart
    if type(h,specfunc(anything,'Product')) then
       func := op(1,h);
       ind := op([2,1],h);
       a := op([2,2,1],h); b := op([2,2,2],h);
       if type(func,ratpoly(anything,ind)) and 
          type(a,integer) and evalb(b=n-1) then
          R := eval(func,ind=n);
       else
          return procname(value(h),args[2..nargs]);
       end if;
    else
       try
          R := eval(h,n=n+1)/h
       catch:
          return false
       end try;
    end if;
  else
    try
       R := h/eval(h,n=n-1)
    catch:
       return false
    end try;
  end if;

  if not type(R,ratpoly(anything,n)) then
    R := traperror(sumtools['simpcomb'](R));
    if R = lasterror then return false end if;
    if not type(R,ratpoly(anything,n)) then
      if hastype(R,specfunc(anything,'Product')) then
         # regular description
         R := value(R)
      end if;
      R := normal(R,'expanded');                           
      if not type(R,ratpoly(anything,n)) then
        return false;
      end if;
    end if;
  end if;
  if rat then
    RR := Normalizer(R)
  end if;
  true
end proc:

   ##############################################################################
   # Name: ApplyLToF (internally used)						#
   # Application of recurrence operator L to F                          	#
   ##############################################################################
   
   ApplyLtoF := proc(L1,F,N,n)
     local An,  d, L, i;
   
     An := OreTools:-SetOreRing(n,'shift');
     d := degree(L1,N);
     L := 'OrePoly'(seq(coeff(L1,N,i),i=0..d));
     OreTools:-ApplyOrePoly(L,F,An)
   
   end proc:



   ##############################################################################
   # Calculating minimal annihilators                                           #
   ##############################################################################

annihilbasis:=proc(p,n,k)
  local c, i, basis, degrees, e;
  c:=[seq(coeff(p,k,i),i=0..degree(p,k))];
  basis:=[]; degrees:=[];
  for e in c do
    while member(degree(e,n),degrees,'i') do
       e:=e-basis[i]*lcoeff(e,n)/lcoeff(basis[i],n);
    end do;
    if not Testzero(e) then 
      basis:=[basis[],e]; 
      degrees:=[degrees[],degree(e,n)];
    end if;   
  end do;
  basis;
end proc:


annihilator:=proc(bas,n,E,R)
  local r,i,j,a,vars,cns,s,L,t,A;
  r:=nops(bas);
  vars:=[seq(a[j],j=0..r)];
  cns:={seq(`+`(seq(vars[j+1]*eval(bas[i],n=n+j),j=0..r)),i=1..r)};
  s:=SolveTools:-Linear(cns,vars);
  s:=adjust(s,vars[r+1]);
  L:=eval([seq(vars[j+1]*E^j,j=0..r)],s);
  t:=R;
  for i from 1 to r do
     L[i+1]:=L[i+1]/t; t:=t*eval(R,n=n+i);
  end do;
  A := collect( numer( `+`( op(L) ) ) , E);
  A;
end proc:


adjust:=proc(sols,x)
  local freevars,s,t,xr,j,ev;
  description "binds free variables in solution keeping x nonzero";
  freevars:=[];
  s:=sols;
  for t in s do
    if lhs(t)=rhs(t) then freevars:=[freevars[],lhs(t)]; end if; 
    if lhs(t)=x then xr:=rhs(t); end if;
  end do;
  if Testzero(xr) then xr:=1; end if;
  for t in freevars do
    j:=0; 
    while Testzero(eval(xr,t=j)) do j:=j+1; end do;
    ev:=t=j;
    xr:=eval(xr,ev);
    s:=eval((s minus {t=t}),ev) union {ev};
  end do; 
  map(proc(y) lhs(y)=normal(rhs(y)); end proc,s);
end proc: