Skip to content

Instantly share code, notes, and snippets.

@goulu
Last active Jul 15, 2016
Embed
What would you like to do?
Faulhaber formula to calculate sum of powers of integers using a generator for Bernouilli numbers
from scipy.special import binom as binomial
def bernouilli_gen(init=1):
"""generator of Bernouilli numbers
:param init: int -1 or +1.
* -1 for "first Bernoulli numbers" with B1=-1/2
* +1 for "second Bernoulli numbers" with B1=+1/2
https://en.wikipedia.org/wiki/Bernoulli_number
https://rosettacode.org/wiki/Bernoulli_numbers#Python:_Optimised_task_algorithm
"""
from fractions import Fraction
B, m = [], 0
while True:
B.append(Fraction(1, m+1))
for j in range(m, 0, -1):
B[j-1] = j*(B[j-1] - B[j])
yield init*B[0] if m==1 else B[0]# (which is Bm)
m += 1
def faulhaber(n,p):
""" sum of the p-th powers of the first n positive integers
:return: 1^p + 2^p + 3^p + ... + n^p
https://en.wikipedia.org/wiki/Faulhaber%27s_formula
"""
s=0
for j,a in enumerate(bernouilli_gen()):
if j>p : break
s=s+binomial(p+1,j)*a*n**(p+1-j)
return s//(p+1)
@tdsymonds

This comment has been minimized.

Copy link

@tdsymonds tdsymonds commented Jul 1, 2016

Just thought I'd let you know that I think you've pasted the faulhaber function in the middle of the bernoulli_gen function, and you need to add the binomial numpy import :) otherwise it's great! Thanks

@goulu

This comment has been minimized.

Copy link
Owner Author

@goulu goulu commented Jul 15, 2016

yes tdsymonds, your're right. Thanks!

BTW this code is part of my Goulib.math2 library (in fact it will be in the forthcoming version...)
I also wrote a blog article about this (in french) here: http://www.drgoulu.com/2016/05/31/pyramides-et-sommes-de-puissances/

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment