Boost logo

Boost :

Subject: [boost] [GSoC 2013] Draft Proposal for Bernoulli Numbers project
From: David Joy (david.joy_at_[hidden])
Date: 2013-05-01 08:39:59

Hi Chris,

Thanks for your reply. I hadn't looked too closely at the Multiprecision
project, and it was only after checking out the material on high-precision
integer arithmetic at MIT Open Courseware
that I realised I had seen something similar once before: a fascinating
by the (Fields Medallist) Timothy Gowers given at Gresham College.

My explorations of this are must wait, though; I feel much more at home in
the world of infinite sums, the Riemann Zeta function, the Gamma function
and the Chinese Remainder Theorem so have decided only to submit a proposal
for the Bernoulli numbers project. Here follows my first draft:


* Add support for Bernoulli numbers along with thread-safe caching.
* Add support for polygamma of positive integer order (requires Bernoulli
* Improve tgamma/lgamma for multiprecision types based on Stirling's approx.
* Optional: Add support for the Hurwitz zeta function.

* Main reference: (Richard P. Brent, David Harvey) Fast computation of
Bernoulli, Tangent and Secant numbers <> ,
plus formulae gleaned from Wikipedia pages (to be cross-checked for

* Bernoulli numbers can be found most useful in finite- or infinite sums,
rather than on their own, e.g.

  - polygamma:
  \psi^{(m)}(z) =
(-1)^{m+1}\sum_{k=0}^{\infty}\frac{(k+m-1)!}{k!}\frac{B_k}{z^{k+m}} \qquad m
\ge 1
  \psi^{(0)}(z) = \ln(z) - \sum_{k=1}^\infty \frac{B_k}{k z^k}

  - Stirling's approximation to (the natural logarithm of) the gamma
  \ln (\Gamma (z)) \sim \left(z-\tfrac{1}{2}\right)\ln(z) -z +
\tfrac{1}{2}\ln(2 \pi) + \sum_{n=1}^\infty \frac{B_{2n}} {2n(2n-1) z^{2n-1}}

  - Hurwitz zeta function (for negative n):
  \zeta(-n,x) = - \frac{B_{n+1}(x)}{n+1}
  where B_{n+1}(x) is the Bernoulli polynomial of degree n+1: B_n(x) =
\sum_{k=0}^n {n \choose k} B_{n-k} x^k
  (for positive n it is related to the polygamma function, above)

therefore the principal aim should be to return a suitable data structure
containing the first n Bernoulli numbers.

* The above paper shows that the first n Bernoulli numbers can be computed
with O( n^2 (log n)^{2 + o(1)} ) bit-operations and provides an algorithm,
FastTangentNumbers, by which the first n Tangent numbers can be computed
with that complexity, whence the first n Bernoulli numbers are obtained via
the relation:

T_n = (-1)^{n-1} 2^{2n} \frac{B_{2n}}{2n}

* The paper also provides another algorithm, TangentNumbers, similar to the
Akiyama-Tanigawa algorithm for Bernoulli numbers but only involving integer
arithmetic, and the suggestion that it may be more efficient than
FastTangentNumbers for n < 1000.

* Implement the algorithms TangentNumbers and FastTangentNumbers from the
reference, using the Boost.Multiprecision library.
* (Trivially) extend these to calculate the first n Bernoulli numbers (and
the first n Bernoulli polynomials?).
* Using m Bernoulli numbers (each of precision k - for some k and m to be
determined) implement the gamma function, polygamma function(s) and Hurwitz
zeta function, all to a given precision n, using the relations given above.

Miscelleous Thoughts
* It would make sense for the Bernoulli numbers to be returned as rational
numbers, using the Boost.Rational library, so that precision is not
needlessly lost by dividing.
* By extension, should the infinite sums needed in the implementation of the
other functions be returned as polynomials which can be evaluated similarly
to the Boost.Polynomial library? (So that we return a single object
representing the function, which the user can then evaluate at multiple

That's as far as I've got at present, but I'll be filling in the remaining
gaps before Friday, of course.



View this message in context:
Sent from the Boost - Dev mailing list archive at

Boost list run by bdawes at, gregod at, cpdaniel at, john at