Efficient implementation of the Hardy-Ramanujan-Rademacher formula

We describe how the Hardy-Ramanujan-Rademacher formula can be implemented to allow the partition function $p(n)$ to be computed with softly optimal complexity $O(n^{1/2+o(1)})$ and very little overhead. A new implementation based on these techniques achieves speedups in excess of a factor 500 over previously published software and has been used by the author to calculate $p(10^{19})$, an exponent twice as large as in previously reported computations. We also investigate performance for multi-evaluation of $p(n)$, where our implementation of the Hardy-Ramanujan-Rademacher formula becomes superior to power series methods on far denser sets of indices than previous implementations. As an application, we determine over 22 billion new congruences for the partition function, extending Weaver's tabulation of 76,065 congruences.Comment: updated version containing an unconditional complexity proof; accepted for publication in LMS Journal of Computation and Mathematic

The complexity of class polynomial computation via floating point approximations

We analyse the complexity of computing class polynomials, that are an important ingredient for CM constructions of elliptic curves, via complex floating point approximations of their roots. The heart of the algorithm is the evaluation of modular functions in several arguments. The fastest one of the presented approaches uses a technique devised by Dupont to evaluate modular functions by Newton iterations on an expression involving the arithmetic-geometric mean. It runs in time $O (|D| \log^5 |D| \log \log |D|) = O (|D|^{1 + \epsilon}) = O (h^{2 + \epsilon})$ for any $\epsilon > 0$, where $D$ is the CM discriminant and $h$ is the degree of the class polynomial. Another fast algorithm uses multipoint evaluation techniques known from symbolic computation; its asymptotic complexity is worse by a factor of $\log |D|$. Up to logarithmic factors, this running time matches the size of the constructed polynomials. The estimate also relies on a new result concerning the complexity of enumerating the class group of an imaginary-quadratic order and on a rigorously proven upper bound for the height of class polynomials

Computing π(N)\pi(N): An elementary approach in O~(N)\tilde{O}(\sqrt{N}) time

We present an efficient and elementary algorithm for computing the number of primes up to $N$ in $\tilde{O}(\sqrt N)$ time, improving upon the existing combinatorial methods that require $\tilde{O}(N ^ {2/3})$ time. Our method has a similar time complexity to the analytical approach to prime counting, while avoiding complex analysis and the use of arbitrary precision complex numbers. While the most time-efficient version of our algorithm requires $\tilde{O}(\sqrt N)$ space, we present a continuous space-time trade-off, showing, e.g., how to reduce the space complexity to $\tilde{O}(\sqrt[3]{N})$ while slightly increasing the time complexity to $\tilde{O}(N^{8/15})$. We apply our techniques to improve the state-of-the-art complexity of elementary algorithms for computing other number-theoretic functions, such as the the Mertens function (in $\tilde{O}(\sqrt N)$ time compared to the known $\tilde{O}(N^{0.6})$), summing Euler's totient function, counting square-free numbers and summing primes. Implementation code is provided

Fast computation of Bernoulli, Tangent and Secant numbers

We consider the computation of Bernoulli, Tangent (zag), and Secant (zig or Euler) numbers. In particular, we give asymptotically fast algorithms for computing the first n such numbers in O(n^2.(log n)^(2+o(1))) bit-operations. We also give very short in-place algorithms for computing the first n Tangent or Secant numbers in O(n^2) integer operations. These algorithms are extremely simple, and fast for moderate values of n. They are faster and use less space than the algorithms of Atkinson (for Tangent and Secant numbers) and Akiyama and Tanigawa (for Bernoulli numbers).Comment: 16 pages. To appear in Computational and Analytical Mathematics (associated with the May 2011 workshop in honour of Jonathan Borwein's 60th birthday). For further information, see http://maths.anu.edu.au/~brent/pub/pub242.htm

Sparse approaches for the exact distribution of patterns in long state sequences generated by a Markov source

We present two novel approaches for the computation of the exact distribution of a pattern in a long sequence. Both approaches take into account the sparse structure of the problem and are two-part algorithms. The first approach relies on a partial recursion after a fast computation of the second largest eigenvalue of the transition matrix of a Markov chain embedding. The second approach uses fast Taylor expansions of an exact bivariate rational reconstruction of the distribution. We illustrate the interest of both approaches on a simple toy-example and two biological applications: the transcription factors of the Human Chromosome 5 and the PROSITE signatures of functional motifs in proteins. On these example our methods demonstrate their complementarity and their hability to extend the domain of feasibility for exact computations in pattern problems to a new level

Computing fast and accurate convolutions

The analysis of data often models random components as a sum of in- dependent random variables (RVs). These RVs are often assumed to be lattice-valued, either implied by the problem or for computational efficiency. Thus, such analysis typically requires computing, or, more commonly, ap- proximating a portion of the distribution of that sum. Computing the underlying distribution without approximations falls un- der the area of exact tests. These are becoming more popular with continuing increases in both computing power and the size of data sets. For the RVs above, exactly computing the underlying distribution is done via a convolu- tion of their probability mass functions, which reduces to convolving pairs of non-negative vectors. This is conceptually simple, but practical implementations must care- fully consider both speed and accuracy. Such implementations fall prey to the round-off error inherent to floating point arithmetic, risking large rela- tive errors in computed results. There are two main existing algorithms for computing convolutions of vectors: naive convolution (NC) has small bounds on the relative error of each element of the result but has quadratic runtime; while Fast Fourier Transform-based convolution (FFT-C) has almost linear runtime but does not control the relative error of each element, due to the accumulation of round-off error. This essay introduces two novel algorithms for these problems: aFFT-C for computing convolution of two non-negative vectors, and sisFFT for com- puting p-values of sums of independent and identically-distributed lattice- valued RVs. Through a rigorous analysis of round-off error and its accumula- tion, both aFFT-C and sisFFT provide control of the relative error similar to NC, but are typically closer in speed to FFT-C by careful use of FFT-based convolutions and by aggressively discarding irrelevant values. Both accuracy and performance are demonstrated empirically with a variety of examples

Time- and Space-Efficient Evaluation of Some Hypergeometric Constants

The currently best known algorithms for the numerical evaluation of hypergeometric constants such as $\zeta(3)$ to $d$ decimal digits have time complexity $O(M(d) \log^2 d)$ and space complexity of $O(d \log d)$ or $O(d)$. Following work from Cheng, Gergel, Kim and Zima, we present a new algorithm with the same asymptotic complexity, but more efficient in practice. Our implementation of this algorithm improves slightly over existing programs for the computation of $\pi$, and we announce a new record of 2 billion digits for $\zeta(3)$