
Boost : 
Subject: [boost] [GSoC 2013] Draft Proposal for Bernoulli Numbers project
From: David Joy (david.joy_at_[hidden])
Date: 20130501 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 highprecision
integer arithmetic at MIT Open Courseware
<http://ocw.mit.edu/courses/electricalengineeringandcomputerscience/6006introductiontoalgorithmsfall2011/lecturevideos/lecture11integerarithmetickaratsubamultiplication/>
that I realised I had seen something similar once before: a fascinating
lecture
<http://www.gresham.ac.uk/lecturesandevents/multiplyinganddividingwholenumberswhyitismoredifficultthanyoumight/>
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 threadsafe 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 crosschecked 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+m1)!}{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(2n1) z^{2n1}}
 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_{nk} 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)} ) bitoperations 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)^{n1} 2^{2n} \frac{B_{2n}}{2n}
* The paper also provides another algorithm, TangentNumbers, similar to the
AkiyamaTanigawa 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/GSoC2013BernoulliNumberstp4646219p4646396.html Sent from the Boost  Dev mailing list archive at Nabble.com.
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk