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,