{{ message }}

Instantly share code, notes, and snippets.

# goulu/faulhaber.py

Last active Jul 15, 2016
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 if m==1 else B# (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 commented Jul 1, 2016 • edited

 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 commented Jul 15, 2016 • edited

 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/