You could probably use scipy.base.polynomial, but it's easy enough to
implement a polynomial yourself. Just use a dict-- each key represents
the power and each value the coefficient of the polynomial.
You didn't say exactly how efficient you need this. It takes only a
couple seconds to sum 100 of the b(k)'s using the implementation below.
This gets you roots out to about A=500.
Perhaps there are bugs in the below implementation. Either way, you
can compare your results and its results and if they're not the same,
well then there's a problem.
------------------------------------------
from rational import rational #get from bitconjurer.org/rational.py
def add(a, b, s=1):
c = {}
for k, v in b.items():
if not a.has_key(k):
c[k] = s*v
for k, v in a.items():
vb = b.get(k, 0)
c[k] = a[k] + s*vb
return c
def raise1(a):
b = {}
for k, v in a.items():
b[k+1] = v
return b
def scale(a, s):
b = {}
for k, v in a.items():
b[k] = v*s
return b
def eval(x, p):
s = 0.0
for k, v in p.items():
s = s + float(v.num)/float(v.den)*x**k
return s
b = {-3 : {}, -2 : {}, -1 : {}, 0 : {0:rational(1,1)}, 1 : {}}
N = 100
for k in range(2,N):
b[k] = scale(raise1(add(b[k-5],b[k-2],-1)),rational(1,k**2))
o = [b[0]]
for k in range(1, N):
o.append(add(o[-1], b[k]))
for x in range(0,800):
print x, eval(x, o[-3]), eval(x, o[-2]), eval(x, o[-1])
# o[-1],o[-2], and o[-3] start to split around x = 500
--
http://mail.python.org/mailman/listinfo/python-list
