Status: Accepted
Owner: skr...@gmail.com
Labels: Type-Enhancement Priority-Medium
New issue 2412 by skr...@gmail.com: Sketch of Sister Celine's algorithm
http://code.google.com/p/sympy/issues/detail?id=2412
Based on the algorithm explained in the book A=B (Chapter 4)
(http://www.math.upenn.edu/~wilf/AeqB.html), I implemented a sketch. At the
end of the chapter there are even some summations that can be solved with
the method, so we even have the test for the method.
There are some other summation algorithms explained on the book.
Things missing:
1) consider summation bounds, and special cases. The algorithm explained
does not consider this case, but I think is not difficult to get them
working.
2) Solve some issues with rsolve (see issue 2408)
3) Solve some issues with simplification and expansion of the Binomial
function.
4) Give the initial parameters to rsolve
The sketch is posted, since I can not dedicate much time to this and can
post to the git this sketch. I hope some one is wiling to take this code
and finish it. I'm open for discussion.
############################################
# SKETCH
############################################
from sympy import *
def findrecur(f, I, J):
i = Symbol("i")
j = Symbol("j")
n = Symbol("n")
a = Symbol("a")
k = Symbol("k")
F = Symbol("F")
y = 0
for ii in xrange(I+1):
for jj in xrange(J+1):
y += a(ii,jj)*(f(n-jj, k-ii)/f(n,k)).simplify()
#y = Sum(Sum(a(i,j)* (f(n-j, k-i)/f(n,k)), (i, 0, I)).doit(), (j, 0, J)
).doit().simplify()
t = together(y)
numer = t.as_numer_denom()[0]
z = collect(numer, k)
print z
#try:
lc = z.as_poly(k).all_coeffs()
#print lc
sols = solve(lc, a(0,0),a(1,0),a(1,1),a(0,1))
#print sols
yy = 0
for ii in xrange(I+1):
for jj in xrange(J+1):
if a(ii,jj) in sols:
yy += sols[a(ii,jj)]*(F(n-jj)).simplify()#, k-ii
else:
yy += a(ii,jj)*(F(n-jj)).simplify()#, k-ii
for ii in xrange(I+1):
for jj in xrange(J+1):
yy = yy.subs(a(ii,jj),1)
print "yy:",yy
return rsolve(yy,F(n))
#except:
# return None
#y = Sum(Sum(b(i,j)* (F(n-j, k-i)), (i, 0, I)).doit(), (j, 0, J)
).doit()
#print y
def f(n,k):
return factorial(n)/(factorial(k)*factorial(n-k))
def g(n,k):
return k*factorial(n)/(factorial(k)*factorial(n-k))
def fa(n,k):
a = Symbol("a")
return (-1)**k*binomial(n,k)*binomial(2*n-2*k,n+a)
def fb(n,k):
x = Symbol("x")
y = Symbol("y")
return binomial(x,k)*binomial(y,n-k)
def fc(n,k):
return k*binomial(2*n+1,2*k+1)
if __name__ == '__main__':
x,y = symbols("x y")
print findrecur(f,1,1)
print findrecur(binomial,1,1)
#error when solving the recurrence
#print findrecur(g,1,1)
#with 1,1 it does not found a recurrence, with 2,2 it gets stuck in
line y += a(ii,jj)*(f(n-jj, k-ii)/f(n,k)).simplify()
#print findrecur(fa,2,2)
#It does not expand binomial(y,n-k)
#print findrecur(fb,1,2)
print findrecur(fc,1,2)
--
You received this message because you are subscribed to the Google Groups
"sympy-issues" group.
To post to this group, send email to sympy-issues@googlegroups.com.
To unsubscribe from this group, send email to
sympy-issues+unsubscr...@googlegroups.com.
For more options, visit this group at
http://groups.google.com/group/sympy-issues?hl=en.
