[boost] [GSoC 2013] Draft Proposal for Bernoulli Numbers project

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
<http://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-006-introduction-to-algorithms-fall-2011/lecture-videos/lecture-11-integer-arithmetic-karatsuba-multiplication/>  
that I realised I had seen something similar once before: a fascinating 
lecture
<http://www.gresham.ac.uk/lectures-and-events/multiplying-and-dividing-whole-numbers-why-it-is-more-difficult-than-you-might/>  
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:

PROJECT PROPOSAL

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

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

* 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
function: 
  \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.

Intentions
* 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
points.)

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

Regards,

David

--
View this message in context: http://boost.2283326.n4.nabble.com/GSoC-2013-Bernoulli-Numbers-tp4646219p4646396.html
Sent from the Boost - Dev mailing list archive at Nabble.com.