[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

 >>> bernoulli(12)
 >>> bernoulli(14)
 >>> bernoulli(16)
 >>> bernoulli(20)

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)+")"
      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 

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,

  def summa(t):
      terms = []   
      b = bernseries(t+1)
      for k in range(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

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)

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)

That's it!


Source code:

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