2939.sagews

# input: q(i) where q(i) is a polynomial
# output: p(n), where p(n)=\sum_{i=0}^n q(i)

i=var('i')
n=var('n')
LMAX=100

q=i^10
p=poly_sum(q)

print "q(i)=",q
print
print "p(n)=sum_{i=0}^n q(i)"
print
print "p(n)=", p, ", or,"
print
print "p(n)=", p.factor()

q(i)= i^10 p(n)=sum_{i=0}^n q(i) p(n)= 1/66*((n - 1)*((n - 2)*((n - 3)*((n - 4)*((n - 5)*((n - 6)*((n - 7)*((n - 8)*(3*(n - 9)*(2*n + 101) + 8470) + 98010) + 603306) + 1974357) + 3256836) + 2404875) + 627022) + 33759) + 66)*n , or, p(n)= 1/66*(n + 1)*(2*n + 1)*(n^2 + n - 1)*(3*n^6 + 9*n^5 + 2*n^4 - 11*n^3 + 3*n^2 + 10*n - 5)*n

# gives the closed form of \sum_{i=0}^n q(i)

def poly_sum(q):
    
    d=q.degree(i)
    R = PolynomialRing(QQ, 'n')
    
    for L in xrange(d,LMAX):
        points=[]
    
        for t in xrange(L):
            points.append((t,sum_q(q,t)))
      
        p=R.lagrange_polynomial(points)
        p2=p.subs(n=n-1)
        
        if (p-p2-q.subs(i=n)).is_zero()==1:
            return p

    return 0
    
# computing sum(t)

def sum_q(q,t):
    
    store=0
    for s in xrange(t+1):
        store+=q.subs(i=s)
    
    return store