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

lecture 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, 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