Working with the recurrence for a sum of squares.
Requirements: ore_algebra package
# Define the ring of shift operators to encode P-recurrence
from ore_algebra import *
Ind.<n> = PolynomialRing(QQ); Shift.<Sn> = OreAlgebra(Ind)
# Define P-recurrence for the sum of k^2 with k from 0 to n
rec_sum = (n+1)^2*Sn^2 - (2*n^2+6*n+5)*Sn + (n+2)^2
rec_sum
(n^2 + 2*n + 1)*Sn^2 + (-2*n^2 - 6*n - 5)*Sn + n^2 + 4*n + 4
# Define P-recurrence for n*(n+1)*(2n+1)/6
rec_closed = n*(n+1)*(n*2+1)*Sn - (n+1)*(n+2)*(2*n+3)
rec_closed
(2*n^3 + 3*n^2 + n)*Sn - 2*n^3 - 9*n^2 - 13*n - 6
# It's easy to verify the closed form satisfies the recurrence derived from the sum
closed = n*(n+1)*(2*n+1)/6
rec_sum(closed)
0
# The following computation shows that rec_closed is a factor of rec_sum, so both the sum and closed form satisfy rec_sum
# Because the sum and closed form are equal for their first two values we then see they are identically equal
Q, R = rec_sum.quo_rem(rec_closed)
print(Q)
print(Q * rec_closed - rec_sum)
((1/2*n + 1/2)/(n^2 + 7/2*n + 3))*Sn + (-1/2*n - 1)/(n^2 + 5/2*n + 3/2) 0
Instead of using the built-in functions of ore_algebra, we can solve an explicit system to derive a recurrence satisfied by the difference of the sum and the closed form, which we then show is zero.
(There is no point -- other than pedagogical completeness -- in doing this when ore_algebra will do it for you, and in less time!)
# Equation for s(n), which encodes the sum
function('s')(n)
eq_sum = rec_sum(s(n))
eq_sum
(n^2 + 2*n + 1)*s(n + 2) - (2*n^2 + 6*n + 5)*s(n + 1) + (n^2 + 4*n + 4)*s(n)
# Equation for c(n), which encodes the closed form
function('c')(n)
eq_closed = rec_closed(c(n))
eq_closed
(2*n^3 + 3*n^2 + n)*c(n + 1) - (2*n^3 + 9*n^2 + 13*n + 6)*c(n)
# We know from dimension counting that s(n) - c(n) satisfies a linear recurrence of order four
# Define this equation with unknown coefficients
var('A B C D')
eq_diff = A * (c(n)-s(n)) + B * (c(n+1)-s(n+1)) + C * (c(n+2)-s(n+2)) + D * (c(n+3)-s(n+3))
eq_diff
D*(c(n + 3) - s(n + 3)) + C*(c(n + 2) - s(n + 2)) + B*(c(n + 1) - s(n + 1)) + A*(c(n) - s(n))
# We re-write this equation, using the recurrences to represent everything in terms of s(n), s(n+1), and c(n)
subs_s = solve(eq_sum,s(n+2))[0]
subs_c = solve(eq_closed,c(n+1))[0]
eq = eq_diff
eq = eq.subs(subs_s.subs(n=n+1))
eq = eq.subs(subs_s)
eq = eq.subs(subs_c.subs(n=n+2))
eq = eq.subs(subs_c.subs(n=n+1))
eq = eq.subs(subs_c.subs(n=n))
eq
D*((2*(n + 2)^2 + 7*n + 20)*(2*(n + 1)^2 + 7*n + 13)*(2*n^2 + 7*n + 6)*c(n)/((2*(n + 2)^2 + n + 2)*(2*(n + 1)^2 + n + 1)*(2*n^2 + n)) + (((n + 1)^2 + 4*n + 8)*s(n + 1) - (2*(n + 1)^2 + 6*n + 11)*((2*n^2 + 6*n + 5)*s(n + 1) - (n^2 + 4*n + 4)*s(n))/(n^2 + 2*n + 1))/((n + 1)^2 + 2*n + 3)) + C*((2*(n + 1)^2 + 7*n + 13)*(2*n^2 + 7*n + 6)*c(n)/((2*(n + 1)^2 + n + 1)*(2*n^2 + n)) - ((2*n^2 + 6*n + 5)*s(n + 1) - (n^2 + 4*n + 4)*s(n))/(n^2 + 2*n + 1)) + B*((2*n^2 + 7*n + 6)*c(n)/(2*n^2 + n) - s(n + 1)) + A*(c(n) - s(n))
# We can solve this for the coefficients of c(n), s(n), and s(n+1)
# There is a family of solutions, giving recurrences satisfied by the difference s(n)-c(n)
sols = solve([eq.coefficient(c(n)), eq.coefficient(s(n)), eq.coefficient(s(n+1))],[A,B,C,D])
sols
[[A == (n^2*(2*r1 + r2) + 2*n*(5*r1 + 2*r2) + 13*r1 + 4*r2)/(n^2 + 2*n + 1), B == -(n^2*(3*r1 + 2*r2) + 6*n*(2*r1 + r2) + 14*r1 + 5*r2)/(n^2 + 2*n + 1), C == r2, D == r1]]
# As expected, one of these solutions is the recurrence rec_sum that both s(n) and c(n) satisfy
function('d')(n)
found_rec = (A * d(n) + B * d(n+1) + C * d(n+2) + D * d(n+3)).subs(sols[0]).subs(r1=0,r2=1).numerator()
print(found_rec)
n^2*d(n + 2) - 2*n^2*d(n + 1) + n^2*d(n) + 2*n*d(n + 2) - 6*n*d(n + 1) + 4*n*d(n) + d(n + 2) - 5*d(n + 1) + 4*d(n)
print((found_rec - rec_sum(d(n))).expand())
0