Compute a Zeilberger pair of a bivariate hypergeometric sequence
  
  Calling sequence:
    Zeilberger( T , n , k , E );
  
  Input:
    T - hypergeometric sequence in n and k
    n - name
    k - name
    E - name of a shift operator ( Ef(n) = f(n+1) )    

  Output:
    [ L , G ]
    L - difference operator of the form
      L = a[d] * E^d + ... + a[1] * E + a[0]
    where a[0], ..., a[d] are polynomials free of k, 
    G - hypergeometric sequence (or zero) such that
      G( n , k+1 ) - G( n , k ) = L T( n , k )
    and the order d is minimal.

  If the global variables _MINORDER and _MAXORDER are specified, the procedure
  will only be looking for operators L of the order d within a 
  _MINORDER.._MAXORDER range.

  In case of homogeneous recurrences, i.e. when L T( n , k ) = 0, an algorithm
  computing minimal annihilators is used.

  References:
    M. Petkovsek, H. Wilf, D. Zeilberger. A=B.
    A.K.Peters, Wellesley, Massachusetts, 1996, Chapter 6.

    Polyakov S.P. On homogeneous Zeilberger recurrences.
    Advances in Applied Mathematics, 2008, vol. 40, no. 1, pp 1-7.