suggestions on how to do this

[bgs]

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