It turns out that rsolve gives pure and simple wrong results :
var("n,p,q")
f,g= (function(u) for u in ("f", "g"))
import sympy
DR=f(p)-f(p-1)-binomial(n,p)
AR=g(p)-g(p-1)-((-1)^p)*binomial(n,p)
SSf=sympy.rsolve(*map(sympy.sympify, (DR, f(p), {f(0):1})))
Sf=SSf._sage_().factor()
SSg=sympy.rsolve(*map(sympy.sympify, (AR, g(p), {g(0):1})))
Sg=SSg._sage_().factor()
At this point,
$$\mathtt{Sf}\,=\,\frac{{n \choose p + 1} + 1}{n + 1} \qquad
\mathtt{Sg}\,=\,\frac{\left(-1\right)^{p} {\left(p + 1\right)} {n \choose p + 1}}{n}$$
This can be compared with the results given by Mathematica :
hypergeometric2f1=function("hypergeometric2f1",
nargs=4,
eval_func=lambda self, *args:hypergeometric(args[:2], [args[2]], args[3]),
print_latex_func=lambda self, *args:latex(hypergeometric(args[:2], [args[2]], args[3])),
conversions={mathematica:"Hypergeometric2F1"})
MMf=mathematica.RSolve([DR==0, f(0)==1], f(p), p)[1][1][2]
Mf=MMf.sage(locals={'Hypergeometric2F1' : hypergeometric2f1})
MMg=mathematica.RSolve([AR==0, g(0)==1], g(p), p)[1][1][2]
Mg=MMg.sage(locals={'Hypergeometric2F1' : hypergeometric2f1})
which gives :
$$\mathtt{Mf}\,=\,-{n \choose p + 1} \,_2F_1\left(\begin{matrix} 1,-n + p + 1 \\ p + 2 \end{matrix} ; -1 \right) + 2^{n} \qquad \mathtt{Mg}\,=\,\frac{\left(-1\right)^{p} {\left(p + 1\right)} {n \choose p + 1}}{n}$$
Sage and Mathematica obtain the same result for the alternating sum, but not for the direct sum. A numerical comparison with the "brute force" sums show that Sage (i. e. Sympy) is wrong and suggests (but does not prove) that Mathematica is right :
sage: matrix(zip([sum(binomial(n,q),q,0,u).subs(n==5) for u in range(6)],[Mf.subs({n:5,p:u}).n() for u in range(6)],[Sf.subs({n:5,p:u}).n() for u in range(6)]))
[ 1 1.00000000000000 1.00000000000000]
[ 6 6.00000000000000 1.83333333333333]
[ 16 16.0000000000000 1.83333333333333]
[ 26 26.0000000000000 1.00000000000000]
[ 31 31.0000000000000 0.333333333333333]
[ 32 32.0000000000000 0.166666666666667]
An analog numerical comparison suggests that the results common to both CASes is right :
sage: matrix(zip([sum(((-1)^q)*binomial(n,q),q,0,u).subs(n==5) for u in range(6)],[Mg.subs({n:5,p:u}).n() for u in range(6)],[Sg.subs({n:5,p:u}).n() for u in range(6)]))
[ 1 1.00000000000000 1.00000000000000]
[ -4 -4.00000000000000 -4.00000000000000]
[ 6 6.00000000000000 6.00000000000000]
[ -4 -4.00000000000000 -4.00000000000000]
[ 1 1.00000000000000 1.00000000000000]
[ 0 0.000000000000000 0.000000000000000]
HTH,
