#################################################################################
#                                                                               #
# MinimalAnnihilator                                                            #
# Compute a minimal annihilating operator of a bivariate hypergeometric term    #
# Input:  T  -  hypergeometric term of n and k                                  #
#         n , k  -  names                                                       #
#         E  -  name of a shift operator (  Ef(n) = f(n+1) )                    #
# Output: difference operator A = a_d*E^d + ... + a_1*E + a_0 such that         #
#         AT(n,k) = 0 and the order d is minimal                                #
#                                                                               #
#                                                Stanislav Polyakov, 2008       #
#                                                                               #
#################################################################################

MinimalAnnihilator := module()
   export ModuleApply;
   local HasAnnihilator, AnnihilatorBasis, Annihilator, Adjust;

   ModuleApply := proc( T , n , k , E )
     local is_hyp, s, p, b, R, A;

     # Type checking
     is_hyp := SumTools:-Hypergeometric:-IsHypergeometricTerm(T,n,'s','up');
     if not is_hyp then
        error "The input is expected to be a hypergeometric term in %1",
               n;
     end if;

     if type(s,freeof(k)) then
        A := denom(s)*E - numer(s);
        return A;
     elif HasAnnihilator(s,n,k,'p') then
        b := AnnihilatorBasis(p,n,k);
        R := s*p/eval(p,n=n+1);
     else
        return FAIL;
     end if;

     A := collect( Annihilator(b,n,E,R) , E );
     A;

   end proc:


   HasAnnihilator := 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
            a := numer(Normalizer(psol));
            true;
         end if;
      end if;
   end proc:


   AnnihilatorBasis := 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;

      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;

      L := map( normal , L );
      numer( convert(L,`+`) );

   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:


end module: