[Edu-sig] Bernoulli Numbers
Kirby Urner
pdx4d@teleport.com
Sun, 20 May 2001 23:29:41 -0700
An algorithm for defining and deriving the Bernoulli numbers
involves starting with Pascal's Triangle for (b-1)^n and writing
a series of equations, each of which defines the next Bernoulli
number in the sequence, using all those that have come before.
b[2] = b[2] - 2*b[1] + 1
b[3] = b[3] - 3*b[2] + 3*b[1] - 1
b[4] = b[4] - 4*b[3] + 6*b[2] - 4*b[1] + 1
These may be re-expressed as follows:
>>> from bernoulli import *
>>> bernseries(7,1)
b[1] = (1/2.) * 1
b[2] = (1/3.) * (3*b[1]- 1)
b[3] = (1/4.) * (6*b[2] - 4*b[1]+ 1)
b[4] = (1/5.) * (10*b[3] - 10*b[2] + 5*b[1]- 1)
b[5] = (1/6.) * (15*b[4] - 20*b[3] + 15*b[2] - 6*b[1]+ 1)
b[6] = (1/7.) * (21*b[5] - 35*b[4] + 35*b[3] - 21*b[2] + 7*b[1]- 1)
b[7] = (1/8.) * (28*b[6] - 56*b[5] + 70*b[4] - 56*b[3] + 28*b[2] - 8*b[1]+ 1)
[1, (1/2.), (1/6.), 0, (-1/30.), 0, (1/42.), 0]
The last line gives a growing list of Bernoulli numbers,
starting with B[0]=1, B[1]=1/2. Notice that all subsequent
odd Bernoullis are 0. We can use another function to call
the above, but just slice off the last member of the list
to get the nth Bernoulli number, which will be rational, a
fraction:
>>> bernoulli(12)
(-691/2730.)
>>> bernoulli(14)
(7/6.)
>>> bernoulli(16)
(-3617/510.)
>>> bernoulli(20)
(-174611/330.)
You can check these with Mathematica, which has the Bernoullis
built in.
The above expressions were developed using regular polynomial
multiplication. In other words, we learn in algebra that
(b+1)*(b+1) = b**2 + 2*b + 1
and that
(b-1)*(b-1) = b**2 - 2*b + 1.
Then you can keep multiplying the latter product by (b-1), getting
polynomials of successively higher degree, and coefficients that
match those in the nth row of Pascal's Triangle (see my previous
post re Pascal's Triangle).
>>> p = Poly([1,-1],'b')
>>> p
b - 1
>>> p*p
b**2 - 2*b + 1
>>> p*p*p
b**3 - 3*b**2 + 3*b - 1
>>> p*p*p*p
b**4 - 4*b**3 + 6*b**2 - 4*b + 1
The algorithm for generating the Bernoullis involves changing
exponents to subscripts, i.e. b**6 or b**5 gets rewritten as b[6]
or b[5] respectively, with b by itself becoming b[1]. So for
example, starting with:
(b-1)**7 =
b**7 - 7*b**6 + 21*b**5 - 35*b**4 + 35*b**3 - 21*b**2 + 7*b - 1
we want that to convert this to something more like:
b[7] =
b[7] - 7*b[6] + 21*b[5] - 35*b[4] + 35*b[3] - 21*b[2] + 7*b[1]- 1
which further simplifies to the expression written below:
b[6] = (1/7.) * (21*b[5] - 35*b[4] + 35*b[3] - 21*b[2] + 7*b[1]- 1)
Consider again:
b[2] = b[2] - 2*b[1] + 1*b[0]
b[3] = b[3] - 3*b[2] + 3*b[1] - 1*b[0]
b[4] = b[4] - 4*b[3] + 6*b[2] - 4*b[1] + 1*b[0]
Taking the last equality, you see the b[4] terms cancel, and we can
bring 4*b[3] onto the left side, then divide by 4:
4*b[3] = 6*b[2] - 4*b[1] + 1*b[0]
b[3] = (1/4.) * (6*b[2] - 4*b[1] + 1*b[0])
Regular expressions proved useful for accomplishing this transformation
of exponents into subscripts e.g. consider the code below:
pat1 = re.compile('\*\*[0-9]*') # find exponents: **n
pat2 = re.compile('b ') # find b (no exponent)
def mksub(match):
return "["+''.join(list(match.group())[2:])+"]"
def bernseries(n,show=0):
"""
Extend Bernoulli series using successive powers
of (B-1)^n, with exponents demoted to subscripts
"""
b = [1] # b[0]=1
poly,term = Poly([1,-1]),Poly([1,-1])
for i in range(n):
poly *= term # get polynomial of next higher degree
bpoly = Poly(poly.coeffs[2:],'b')
bexpr = re.sub(pat1,mksub,str(bpoly)) # b**3 -> b[3]
bexpr = re.sub(pat2,"b[1]",bexpr) # b -> b[1]
if show: print ("b["+str(i+1)+"] = " +
str(Fraction(1,-poly.coeffs[1])) +
" * (" + bexpr)+")"
b.append(Fraction(eval(bexpr),-poly.coeffs[1]))
return b
In searching a string like '10*b**3 - 10*b**2 + 5*b- 1' with regular
expression '\*\*[0-9]*', we're looking for all occurances of **
followed by one or more digits, i.e. the exponent. The \*\* escapes
the double asterisk, to mean the literal '**', whereas the unescaped
asterisk following [0-9] says the pattern of digits 0 thru 9 may repeat
until it doesn't i.e. until whitespace is encountered.
So in '10*b**3 - 10*b**2 + 5*b- 1' we first find '**3', and what we
do next is give re.sub the name of a function for dealing with our
located string. How the substitution gets carried out depends on this
function, which automatically gets, as its only argument, the match
object returned when the pattern is found.
match.group() will disclose what exactly was found, and we need
that information, because we want to throw away the asterisks,
but *keep the digits*. "["+''.join(list(match.group())[2:])+"]"
is what does the job: stick square brackets at the start and
end, and slice off the stars (asterisks) in the middle. What
gets returned is [3] instead of **3, or [23] instead of **23 or
[444] in place of **444 etc. And re.sub does this work for the
whole length of the string, so after pat1 is through, we're left
with: '10*b[3] - 10*b[2] + 5*b- 1'
It's that final b that we want to deal with next, and pat2 does
this. Find 'b ' and replace it with 'b[1]' -- very simple.
Each time through the loop, we multiply our polynomial by
(x-1), bumping it up another degree. Then we skip the first two
coefficients and build a 2nd polynomial in terms of b, run the
regular expression substitions and then actually evaluate the
result as a fraction, in order to secure the next Bernoulli
number for our list (which list we consult at each pass, as
every Bernoulli is based on those that have come before).
I'm not claiming this is the best or fastest way to get the
Bernoullis. It just so happens that we have a polynomial object,
so successive multiplications by (x-1) is easily done. More usually,
the binomial theorem is used to generate the coefficients directly,
using n!/k!(n-k)!, where n is the degree of the polynomial, and k
the term. Here, we get those coefficents from the polynomial
objects.
So we're looking at mathematical relationships, using some open
source math objects already at our disposal, plus finding an
opportunity to use regular expressions in a not-too-complicated,
introductory, yet math-relevant way.
What good are Bernoulli numbers? Well, according to 'The Book of
Numbers' by Conway and Guy, it was actually a dude named Johann Faulhaber,
writing in 1631, who came up with the most familiar way in which the
so-called Bernoulli numbers are used: to generate a polynomial
giving the sum of the first n consecutive integers all raised to
the kth power. You get a polynomial of degree k+1.
In other words, suppose I want an expression for 1**10 + 2**10 +
3**10 + 4**10... + n**10, what would I write? My bernoulli.py module
will give the answer if you just invoke summa():
def choice(n,k):
return reduce(mul,
[1L]+range(n,n-k,-1))/reduce(mul,[1L]+range(2,k+1))
def summa(t):
terms = []
b = bernseries(t+1)
for k in range(t+1):
terms.append(Fraction(b[k]*choice(t+1,k),t+1))
return Poly(terms+[0])
[choice(n,k) is just n!/k!(n-k)! and figures into the coefficients, along
with the Bernoulli numbers, as per the above summa code]:
>>> summa(10)
(1/11.)*x**11 + (1/2.)*x**10 + (5/6.)*x**9 - x**7 + x**5 - (1/2.)*x**3
+ (5/66.)*x
The above polynomial will evaluate the the sum of the first x numbers
to the 10th power.
OK, so if I substitute 3 for x, I'll know the sum of 1**10 + 2**10 + 3**10.
Let's do that at the command line first:
>>> 1**10 + 2**10 + 3**10
60074
OK, now since summa() returns a Polynomial object, and I can evaluate
those simply be passing a value for x:
>>> p10 = summa(10)
>>> p10(3)
60074
How about the sum of the 10th powers of consecutive integers up to 1000
(inclusive). Conway and Guy give the answer on pg. 108 of 'The Book of
Numbers'. Let's see if our matches:
>>> p10(1000L)
91409924241424243424241924242500
That's it!
Kirby
Source code:
http://www.inetarena.com/~pdx4d/ocn/python/bernoulli.py
http://www.inetarena.com/~pdx4d/ocn/python/mathobjects.zip
mathobjects contains pyfraction, polynomial and simplematrix,
all of which have source code in ocn/python (as well as HTML
versions), but which are designed to operate in a subdirectory
named mathobjects with the included __init__.py i.e. as a
package.