I'm trying to write a reciusive function to compute $$S_p=1^p+2^p+3^p+\dots+n^p,$$ where both $n$ and $p$ are positive integers. My Sage code correctly returns $S_1$, but does not work for $p>1$. Any ideas?
I am not a big Sage user, but I am using it as a way to introduce pre-calculus students to CAS. Hopefully, I will be able to share this post with my students, and certainly, give credit to anyone that helps fix this mess.
Here's the code:
sage: reset()
sage: k,n,p=var('k,n,p')
sage: assume(k,'integer')
sage: assume(n,'integer')
sage: assume(p,'integer')
sage: def S(p):
....: if p == 1:
....: return n*(n+1)/2
....: else:
....: return ((n+1)^(p+1) - n - 1 - sum(binomial(p+1,k)*S(p-k+1),k,2,p))
....: /(p+1)
....:
sage: S(1)
1/2*(n + 1)*n
sage: S(2)
