By Claude:
The Zeilberger package (Gosper, GosperSum, Zeilberger, parGosper) and simplify_sum, which uses it, return wrong results when the shift quotient of the summand contains factors that don't depend on the summation variable k and differ by a nonnegative integer, like p and p-1, or 3 and 2.
(%i1) display2d:false$
(%i2) load(simplify_sum)$
(%i3) simplify_sum(sum(binomial(n,k)*p^k*(p-1)^(n-k),k,0,n));
(%o3) (2*p-2)^n
(%i4) subst([n = 2,p = 3],%);
(%o4) 16
(%i5) sum(binomial(2,k)*3^k*2^(2-k),k,0,2);
(%o5) 25
(%i6) simplify_sum(sum(binomial(n,k)*3^k*2^(n-k),k,0,n));
(%o6) 'sum(2^(n-k)*3^k*binomial(n,k),k,0,n)
(%i7) Zeilberger(binomial(n,k)*p^k*(p-1)^(n-k),k,n);
(%o7) [[(k*(p-1))/(n-k+1),[2*(p-1),-1]]]
(%i8) GosperSum(3^k/2^k,k,0,m);
(%o8) 2^(-m-1)*3^(m+1)*(m+1)
(%i9) subst(m = 0,%);
(%o9) 3/2
By the binomial theorem the first sum is (p+(p-1))^n = (2*p-1)^n, so at n = 2, p = 3 it is 25, not 16. The second sum is 5^n, which simplify_sum doesn't find. The recurrence from Zeilberger should have the coefficient 2*p-1, not 2*(p-1). GosperSum(3^k/2^k,k,0,0) is the single term 1, not 3/2. Also Gosper(3^k/2^k,k) returns k, while the certificate is 2, and simplify_sum(sum(binomial(n,k)*(p+1)^k*p^(n-k),k,0,n)) returns 2^n*p^n instead of (2*p+1)^n.
The factors involved don't depend on k, so they can't be integer shifts of each other in k, but the Gosper form computation treats them as such.
If I remember correctly,
nusumis another implementation of Gosper's method, so we could try to usenusumto check the result ofGosper.