Instantly share code, notes, and snippets.

# goulu/faulhaber.py

Last active July 15, 2016 11:51
Show Gist options
• Save goulu/5bbf24a3e2e25070904b79f49020448f to your computer and use it in GitHub Desktop.
Faulhaber formula to calculate sum of powers of integers using a generator for Bernouilli numbers
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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 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...)