# recurrence-rsolve-check

I am was testing some ideas and don't understand where this went wrong.

from sympy import Function, rsolve
from sympy.abc import n
y = Function('y');z=Function('z')
f = y(n)-y(n-1)-(-1)^n*binomial(x,n);g=z(n)-z(n-1)-binomial(x,n)


f is the sume of alternating binomial coefficients which is known.

g is just the summation of binomial coefficients (purportedly nonexistent)

frec = rsolve(f,y(n),{y(0):1});grec = rsolve(g,z(n),{z(0):1})
frec


--

zrec


Now this is correct, the alternating sum of binomials

fdiff = frec-(frec.subs({n : n-1}))
simplify(fdiff)


Whereas this is not correct?

gdiff = grec-(grec.subs({n : n-1}))
simplify(gdiff)


Am I missing something in rsolve usage?

edit retag close merge delete

Sort by » oldest newest most voted

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,

more

Guess I will file a bug report! And... see if I can find my way through the code (: I have programmed in >10 languages; and the result is, I detest programming. Sigh...

( 2021-03-03 22:47:49 +0100 )edit

The bug is in sympy (that can be checked by comparing SSf and its Sage translation Sf) ; therefore, the bug should be filed as an issue against sympy. A Sage ticket could be useful to point it.

An alternative could be an independent implementation of rsolve, which is a much larger task (at least an order of magnitude...).

( 2021-03-03 23:40:26 +0100 )edit

Thanks, unfortunately it doesn't terminate for n not positive integer; which I was originally looking for. But I am pretty sure I have a straightforward way to finding a terminating Generalized Hypergeometric form. Still accumulates terms, but the number grows by n, so It's not too bad. I will look at rsolve though. I am doing these things as an amateur hobbyist.

( 2021-03-05 01:06:38 +0100 )edit