# Let # # D_{n,k} := Q(n,k)[ Sn, Sk, Sn^(-1), Sk^(-1) ] # # where Sn sends f(n,k) to f(n+1,k) and Sk sends f(n,k) to f(n,k+1). H0 := binomial(n,k)^7; # H0 is a hypergeometric term, which is equivalent to saying that # # Omega := Q(n,k) * H0 # # is a D_{n,k}-module. Let # # M := Omega / Delta_k(Omega) # # where Delta_k = Sk - 1. This M is a D-module where # # D := Q(n)[Sn, Sn^(-1)]. # # Note: since we can always multiply on the left with a power of Sn, we may assume # without loss of generality that an non-zero element L in D is actually in Q(n)[Sn] # and that its coefficient of Sn^0 is non-zero. H := H0/(2*n+3*k); # H is an element of Omega. Let # # Hbar = H + Delta_k(Omega) # # be the image of H in M. # Definition: # # The telescoper of H is the minimal operator L in D for which L( Hbar ) = 0. # # As mentioned in the above note, we may assume that L is in Q(n)[Sn]. The condition # L( Hbar ) = 0 is equivalent to saying that L(H) is in Delta_k(Omega), in other # words, L(H) equals Delta_k(c) for some c in Omega, called the certificate. # Zeilberger's algorithm computes both L and the certificate, as follows: # # with(SumTools): with(Hypergeometric): # v := Zeilberger(H, n, k, Sn): # Then v[1] = L and v[2] is the certificate. # Let D[k] = Q(n)[k, Sn, Sn^(-1)] and let # # N = D[k] * H0bar where H0bar = H0 + Delta_k(Omega). # # This N is a D-submodule of M. It consists of all elements of M with no k in the # denominator (other than those introduced by applying Sn or Sn^(-1) to H0). # # Footnote: remember that we're working modulo Delta_k(Omega), and some denominators # can be reduced away by subtracting suitable elements of Delta_k(Omega). However, # the denominator 2*n+3*k seen in H cannot be reduced away mod Delta_k(Omega) # so Hbar is not in N. # Our goal is L with L(Hbar) = 0 in M. Then L(Htilde) = 0 in M/N where # # Htilde = Hbar + N (an element of M/N). # # L should annihilate Hbar and hence Htilde, so the minimal annihilator of Htilde # must be a right-factor of L. # Annihilating Htilde means bringing Hbar into N, which simply means getting # rid of that denominator, the one that prevented it from being in N. # # To eliminate this denominator, we have to subtract functions with the same # denominator and make sure that the residues at that denominator cancel out. # # This function has the same denominator as H: H_shift := subs(n=n+3, k=k-2, H); # Now lets compute the quotient of their residues at the denominator 2*n + 3*k, # in other words, at k = -2/3 * n. RatFunction := simplify(convert(H_shift / H, GAMMA)): r := factor(subs(k = -2/3 * n, RatFunction)); # Now the denominator of # # H_shift - r * H # # cancels out; the terms have the same denominator 2*n + 3*k, and the factor r was # computed to make sure that their residues cancel out. # Note that Sk acts trivially on M. After all, we're working modulo Delta_k(Omega) # and Delta_k = Sk - 1. # # So the image of H_shift in M is equal to Sn^3( Hbar ). Therefore, the # denominator 2*n + 3*k in Sn^3( Hbar ) - r * Hbar cancels out, in other words # # R( Hbar ) is in N, or equivalently, R( Htilde ) = 0 in M/N. # # where R := Sn^3 - r; R := collect(primpart(R, Sn), Sn, factor); # to clear out denominators # R( Htilde) = 0 in M/N. This cannot happen for any R' with order < 3; the denominators # of H, Sn(H), Sn^2(H) are not in the same k-shift-equivalence class, and so they cannot # cancel mod Delta_k(Omega). So R is the minimal annihilator of Htilde, and must therefore # right-divide the minimal annihilator L of Hbar. # # The corresponding left-factor, call it L', is the telescoper of R(H), in other words, # the minimal annihilator of R( Hbar ). # # Now R( Hbar ) is in N, and as the above computation showed, submodules can lead # to factors of the telescoper. And this can be very efficient. The computation of # the order-3 right-factor R of L only took 0.027 seconds! This is very much faster # than computing L with Maple's Zeilberger algorithm (and then factoring it). # # This motivates us to look for more submodules. Consider the automorphism: # # phi: M --> M which sends k to n - k (*) # # In general, the substitution (*) does not need to be an automorphism of M, but in # our example it is, because H0(n, k) = H0(n, n-k) and H0 generates Omega. # # Let phi+ = (1 + phi)/2 and let phi- = (1-phi)/2. # # These are idempotents (e is an idempotent if e^2 = e). # Their sum is 1 and their product is 0 (use phi^2 = 1). # # In general, idempotents correspond to a direct sum decomposition: # # N = N+ oplus N- # # where N+ is the image of phi+ (equals the kernel of phi-) # and N- is the image of phi- (equals the kernel of phi+). # # Next, lets find a basis of N. Standard techniques used in telescoping algorithms # to reduce expressions modulo Delta_k( Omega ) can be used to show that every # element of N can be written as: # # a0 + a1 * k + ... + a6 * k^6 + Delta_k( Omega ). # # Now N+ consists of the submodule invariant under phi: k --> n-k. # The subring of Q(n)[k] invariant under that is Q(n)[ u ] where u := k*(n-k). # So a basis of N+ will be: # # B+ = {u^0, u^1, u^2, u^3} ( mod Delta_k(Omega) ) # # Now N- is anti-invariant under phi, so a basis of N- is # # B- = {v, u * v, u^2 * v} # # where v := k - (n-k), which is anti-invariant under phi, i.e. phi(v) = -v. # # Now B := B+ union B- is a basis of N. Later in this file we will write R( Hbar ) # as a Q(n)-linear combination of B. That way we can easily read off the components # of R( Hbar ) in N+ and in N-. # # Let L_plus and L_minus be the minimal annihilators of the components of R( Hbar ) # in N+ and in N-, then L' := LCLM( L_plus, L_minus ) will be minimal annihilator # of R( Hbar ), i.e. L' will be the telescoper of R(H), and then L' * R will be # the telescoper of H. # # To compute L_plus and L_minus, first we write R( Hbar ) in terms of the basis B, # and after that, we compute the action of Sn on B+ and B-. # Computing R( Hbar ), i.e., computing R(H) and reducing it mod Delta_k( Omega ). # Recall that k-shifts act trivially when we work mod Delta_k( Omega ), so the # shift k -> k-2 is harmless: RH := lcoeff(R, Sn) * subs(n=n+3, k=k-2, H) + tcoeff(R,Sn) * H: # Divide R(H) by H0 to get a rational function: RHd := normal(simplify(convert(subs(k=n-k, RH/H0), GAMMA))): # Remark (***) The k -> n-k in the above line makes it come out a bit nicer; I don't have # automated reduction code and this k -> n-k helps to make the reduction in the next lines # a bit simpler. But a reduction implementation should equally well work without this k -> n-k # substitution. # To compute with RHd, keep in mind that since we dropped the factor H0, so if we apply Sk # then we should remember to also multiply by Sk(H0)/H0 = ((n-k)/(k+1))^7. # # That leads to the following ad-hoc reduction of R(H), as represented by RHd = R(H)/H0, # modulo elements G in Delta_k( Omega ). for j from 4 to 0 by -1 do G := add(c[i]*k^i, i=0..6)/(k+j)^`if`(j=0,0,7); # in Omega with the factor H0 omitted G := subs(k=k+1,G) * ((n-k)/(k+1))^7 - G; # in Delta_k(Omega) with H0 omitted eq := {coeffs(rem(numer(normal( (RHd - G)*(k+j+1)^7 )), (k+j+1)^7, k),k)}; # Solve eq to find out what the c[i] in G should be, then subtract G to reduce RHd: RHd := normal(RHd - subs(solve(eq,{seq(c[i],i=0..6)}), G)); od: u := k*(n-k); v := k - (n-k); {coeffs(collect(RHd - add(c[i] * u^i, i=0..3) + v * add(d[i] * u^i, i=0..2), k),k)}: Decomp := factor(solve(%, {seq(c[i],i=0..3), seq(d[i],i=0..2)})): # # This Decomp gives us the projection of RH on N+ and N- # Specifically, the c[i] in Decomp give the weights for the B+ basis elements and the d[i] # give the weights of the B- basis elements. # Here we combine the basis elements in B+ into one object by putting variables c[i] in # front of them and then taking the sum: BP := add(c[i] * u^i, i=0..3): # This way we can apply Sn to all elements of B+ at once (not that this is faster than # doing this one basis element at a time, just saving a bit of programming): SnBP := subs(n = n+1, BP) * ((n+1)/(n-k+1))^7: SnBP := subs(k = n-k, SnBP): # See (***) above # Take a generic element of Delta_k( ... ) with the factor H0 removed: G := add(e[i]*k^i, i=0..6): G := subs(k=k+1, % ) * ( (n-k)/(k+1) )^7 - %: sol := solve({coeffs(rem(normal( (SnBP - G) * (k+1)^7 ), (k+1)^7, k), k)}, indets(G) minus {k,n}): SnBP := normal(SnBP - subs(sol, G)): # Reduction mod Delta_k( Omega ). SnBP := evala(subs(k = RootOf( u - U, k), SnBP)): # Write SnBP in terms of u instead of k. M := Matrix( [seq([seq(factor(coeff(coeff(SnBP, c[i]),U,j)),i=0..3)], j=0..3)] ); # This matrix describes the action of Sn on the basis B+ # The projection of R(H) on N+ written in terms of basis B+ is given by this vector V[0]: V[0] := [seq(subs(Decomp, c[i]), i=0..3)]: # Apply Sn to this 4 times, using the matrix M that describes the action of Sn on B+ for i to 4 do V[i] := map(factor, convert(M . Vector(subs(n=n+1, V[i-1])),list)) od: # Computing a linear dependence between V[0] .. V[4] gives L_plus: L_plus := subs(solve({op(add(c[i] *~ V[i], i=0..4))}, {seq(c[i],i=0..4)}), add(c[i]*Sn^i,i=0..4)): L_plus := collect(primpart(L_plus, Sn), Sn): a := proc(N::posint) global H,k,n; add(eval(H, n=N), k=0..N) end; seq(a(i), i=1..10); # The component N- contributes zero to the sequence, and therefore, L_plus * R will # annihilate the sequence a(1), a(2), ... # Indeed, it is the minimal recurrence for this sequence. # # So for computing a recurrence for the sequence, we do not need L_minus. # However, strictly speaking L_plus * R is not a telescoper for H because # it does not send Hbar to zero, due to the fact that the N- component is not # annihilated by L_plus * R. # The "full" telescoper, the minimal annihilator of Hbar, is equal to: # # LCLM( L_plus, L_minus ) * R. # # While L_minus is not useful if we want the minimal recurrence for the sequence, # lets compute it anyway so that we can compare with the "full" telescoper as computed # by Maple's Zeilberger algorithm. # Here we combine the basis B- into one object by putting variables d[i] in # front of them and then taking the sum: BM := v * add(d[i] * u^i, i=0..2): # Now apply Sn to BM: SnBM := subs(n = n+1, BM) * ((n+1)/(n-k+1))^7: SnBM := subs(k = n-k, SnBM): # See (***) above sol := solve({coeffs(rem(normal( (SnBM - G) * (k+1)^7 ), (k+1)^7, k), k)}, indets(G) minus {k,n}): SnBM := normal(SnBM - subs(sol, G)): # Reduction mod Delta_k( Omega ). SnBM := evala(subs(k = RootOf( u - U, k), -SnBM/v)): # Write SnBM in terms of B- M := Matrix( [seq([seq(factor(coeff(coeff(SnBM, d[i]),U,j)),i=0..2)], j=0..2)] ); # Describes the action of Sn on the basis B- V[0] := [seq(subs(Decomp, d[i]), i=0..2)]: # Projection of R(H) on N- written in terms of basis B- for i to 3 do V[i] := map(factor, convert(M . Vector(subs(n=n+1, V[i-1])),list)) od: L_minus := subs(solve({op(add(d[i] *~ V[i], i=0..3))}, {seq(d[i],i=0..3)}), add(d[i]*Sn^i,i=0..3)): L_minus := collect(primpart(L_minus, Sn), Sn): # The "full" telescoper is now LCLM( L_plus, L_minus ) * R. lprint("Total CPU time for the full telescoper in factored form", time()): with(LREtools): _Env_LRE_tau := Sn; _Env_LRE_x := n; L_full_expanded := MultiplyOperators(LCLM( L_plus, L_minus ), R): lprint("Total CPU time for the full telescoper in expanded form", time()): # Now compare with Maple's Zeilberger: _MAXORDER := 20; L_Zeilberger := SumTools:-Hypergeometric:-Zeilberger(H, n, k, Sn)[1]: evalb( normal(L_full_expanded - L_Zeilberger) = 0 ); # Found the same telescoper L.