[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.