[Edu-sig] Re: quadratic functions

Kirby Urner urner@alumni.princeton.edu
Fri, 12 Oct 2001 22:12:20 -0700


On Fri, 12 Oct 2001 12:08:13 -0700, I wrote:

>And is there a more general expression for
>
>               n
>   f(n,c) =3D  SIGMA k^c =3D some polynomial of degree c+1 ?
>             k =3D 1
>
>
>The answer to the 2nd question is yes, and the Bernoulli numbers
>give us the coefficients of such a c+1 degree polynomial. =20

Just to spell that out at the command line, what we're=20
looking for is a closed form expression for the sum of
consecutive integers all to the cth power.  What we=20
expect is some polynomial of degree (c+1) with rational
coefficients in which the Bernoulli numbers figure (as
fractions).

So, for example, what is 1^4 + 2^4 + 3^4 + 4^4... n^4?

An easy function just does the arithmetic in a loop:

  >>> from operator import add
  >>> def sum(n,c):
	  return reduce(add,[i**c for i in range(1,n+1)])

Note that range(1,n) generates the list [1,2...n-1], hence
the parameter n+1.  [i**c for i in range(1,n+1)] is the=20
sequence of terms in list form, and reduce(add, list)=20
simply sums them up.

  >>> sum(10,4)
  25333

That's the sum of the first 10 natural numbers to the 4th power.

Bernoulli got his numbers from Johann Faulhaber, who wrote=20
about this series in 1631.[1]  With this series, Bernoulli=20
could brag about how quickly he could get the sum of the 10th
powers of the first 1000 natural numbers. =20

Using our sum function:

  >>> sum(1000,10)
  91409924241424243424241924242500L

(the L means 'long integer' -- it'll be dropped in later=20
versions of Python, as the integer types become more seamless).

But Bernoulli was using a polynomial:

  >>> 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 coefficients above aren't naked Bernoulli numbers; they're=20
folded into products which also involve choice(n+1,k), where
k is the kth term. =20

  kth coefficient =3D b[k]*choice(n+1,k)   (b[k] =3D kth bernoulli)
                    ------------------
                         n+1

Given x=3D1000, the powers are easy to compute.  summa(10)=20
returns a polynomial object into which we can substitute a=20
value using p() syntax i.e.:

  >>> p =3D summa(10)
  >>> p(1000)
  91409924241424243424241924242500

Let's look at the first 10 bernoullis:

  >>> bernseries(10)
  [1, (1/2), (1/6), 0, (-1/30), 0, (1/42), 0, (-1/30), 0, (5/66)]

Not so bad.  The odd ones become 0 after the 0th, 1st, 2nd.

The Bernoulli numbers get pretty crazy (see below) and difficult
to compute by hand after the first few.  But our computer=20
algorithm, using the Fraction object, is up to the task.=20

Let's ask for the 100th Bernoulli number:

  >>> bern(100)
  (-94598037819122125295227433069493721872702841533066936
  133385696204311395415197247711/33330)

Wow, that took almost 10 seconds.  And it agrees with the result=20
posted here (yay):

ftp://ftp.cdrom.com/pub/gutenberg/etext01/brnll10.txt

We calculate successive bernoullis using Pascal's Triangle.
The bernseries algorithm will show how this works if passed=20
a non-zero "show" parameter (2nd argument):

  >>> bernseries(4,1)
  b[1] =3D (1/2) * (1)
  b[2] =3D (1/3) * (3*b[1]- 1)
  b[3] =3D (1/4) * (6*b[2] - 4*b[1]+ 1)
  b[4] =3D (1/5) * (10*b[3] - 10*b[2] + 5*b[1]- 1)
  [1, (1/2), (1/6), 0, (-1/30)]

Those'd be the bernoullis used to compute the sum of=20
consecutive integers to the 4th power.  The corresponding
5th degree polynomial:

  >>> summa(4)
  (1/5)*x**5 + (1/2)*x**4 + (1/3)*x**3 - (1/30)*x

And the result for x =3D 10?

  >>> p =3D summa(4)
  >>> p(10)
  25333

That agrees with our earlier answer.
                                          n
And note we get the earlier formula for SIGMA k^2
                                        k =3D 1

>Well, you have metal cluster experiments, written up in=20
>Scientific American (Dec 1998), about agglomerations of atoms
                          ^^^^
                          1997

>in icosahedral form transforming into cuboctahedra.[1]  How many
>atoms in all?  These formula will tell you.

So how do we do the total number of spheres in a cuboctahedron
again?  I posted quite a bit about this over on math-teach=20
awhile ago.

=46or each layer, you've got the six squares and 8 triangles,=20
but then you have shared edges and vertices, so adding all=20
the squares and triangles involves double-counting edge=20
spheres and quadruple counting the vertices.  The number of
edge spheres is 24 * (n-2) and the number of vertices is 12
(see http://www.inetarena.com/~pdx4d/ocn/urner.html for a
graphic).  So 8*n(n+1)/2 + 4*n^2 - 24*(n-2) - 3*12 should
be it.  That simplifies to:  10*n^2 - 20*n + 12.  This=20
works for n>=3D2, i.e. a cuboctahedron with 2 spheres along
an edge has 12 spheres in that layer, with 3 spheres along
an edge has 42, then 92, then 162 and so on (icosahedral
shells the same, as per the jitterbug transformation:
http://www.inetarena.com/~pdx4d/ocn/numeracy0.html ).

What we want is a cumulative number of spheres in all=20
layers.  In other words:

        n
  1 + SIGMA (10k^2 - 20k + 12)
       k=3D2

the 1 being the nuclear sphere.

       n
  10 SIGMA k^2 =3D 10 [(1/3)*n^3 + (1/2)*n^2 + (1/6)*n) - 1]
      k=3D2

       n
  20 SIGMA k =3D  20 [n(n+1)/2 - 1]
      k=3D2

       n
  12 SIGMA 1 =3D 12[n - 1]
      k=3D2

Combining and simplifying gives:

  Cubocta(n) =3D (10/3)*n^3 - 5*n^2 + (11/3)*n - 1

(remember to add 1 in front of all the SIGMA terms).

Making this a function:

  >>> def cubocta(n):
  	  return int((10/3)*n**3 - 5*n**2 + (11/3)*n - 1)

  >>> map(cubocta,range(1,11))
  [1, 13, 55, 147, 309, 561, 923, 1415, 2057, 2869]

That's the first 10 cuboctahedral numbers.

Citation:
http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi=
?Anum=3DA005902

=3D=3D=3D=3D=3D=3D=3D=3D=3D

ID Number: A005902 (Formerly M4898)
Sequence:  1,13,55,147,309,561,923,1415,2057,2869,3871,5083,6525,8217,
           10179,12431,14993,17885,21127,24739,28741,33153,37995,43287,
           49049,55301,62063,69355,77197,85609,94611,104223,114465,
           125357,136919,149171
Name:      Centered icosahedral (or cuboctahedral) numbers, also crystal =
ball
              sequence for f.c.c. lattice.
Comments:  Called "magic numbers" in some chemical contexts.
References S. Bjornholm, Clusters..., Contemp. Phys. 31 1990 pp. 309-324.
           T. P. Martin, Shells of atoms, Phys. Reports, 273 (1996), =
199-241,
eq.
              (2).
           B. K. Teo and N. J. A. Sloane, Magic numbers in polygonal and
polyhedral
              clusters, Inorgan. Chem. 24 (1985), 4545-4558.
Links:     J. H. Conway and N. J. A. Sloane, Low-Dimensional Lattices =
VII:
Coordination Sequences, Proc. Royal Soc. London, A453 (1997), 2369-2389 =
(pdf,
ps).
=46ormula:   (2*n+1)*(5*n^2+5*n+3)/3.
Maple:     A005902:=3D n -> (2*n+1)*(5*n^2+5*n+3)/3;

=3D=3D=3D=3D=3D=3D=3D=3D=3D

Kirby

PS:  Source code for bernoulli.py is at the Oregon Tips website,
along with the mathobjects package upon which it depends:
http://www.oregon-tips.org/

[1]  See John H. Conway, Richard K. Guy, 'The Book of Numbers'=20
(Springer-Verlag, 1996) p. 106-109 re Faulhaber and Bernoulli=20
numbers.