What is the time complexity of computing sin(x) to t bits of precision?

Question

Short version of the question: Presumably, it's poly$(t)$. But what polynomial, and could you provide a reference?

Long version of the question: I'm sort of surprised to be asking this, because it's such an extremely basic sounding question. Here are some variants on it:

1. How much time does it take to compute $\pi$ to $t$ bits of precision?
2. How much time does it take to compute $\sin(x)$ to $t$ bits of precision?
3. How much time does it take to compute $e^x$ to $t$ bits of precision?
4. How much time does it take to compute $\mathrm{erf}(x)$ to $t$ bits of precision?
5. How much time and how many random bits does it take to generate a (discrete) random variable $X$ such that there is a coupling of $X$ with a standard Gaussian $Z \sim N(0,1)$ for which $|X - Z| < \delta$ except with probability at most $\epsilon$?

In my area of theory of computer science, no one seems to pay much attention to such questions; an algorithm description might typically read

"Generate two Gaussian random variables $Z$ and $Z'$ and check if $\sin(Z \cdot Z') > 1/\pi$"

or some such thing. But technically, one should worry about the time complexity here.

One colleague of mine who's more of an expert on these things assured me that all such "calculator functions" take time at most poly$(t)$. I well believe that this is true. But again, what polynomial (out of curiosity, at least), and what is the reference?

I kind of assumed that the answers would be in every single numerical analysis textbook, but I couldn't find them there. It seems (perhaps reasonably) that numerical analysis cares mainly about getting the answers to within a fixed precision like 32 or 64 bits or whatever. But presumably somebody has thought about getting the results to arbitrary precision, since you can type

Digits := 5000; erf(1.0);

into Maple and it'll give you an answer right away. But it seemed hard to find. After much searching, I hit upon the key phrase "unrestricted algorithm" which led me to the paper "An unrestricted algorithm for the exponential function", Clenshaw-Olver-1980. It's pretty hard to read, analyzing the time complexity for $e^x$ in terms of eight (??!) parameters, but its equation (4.55) seems to give some answers: perhaps $\tilde{O}(t^2)$ assuming $|x|$ is constant?

And really, all that work for little old $e^x$? As for erf$(x)$, I found the paper "The functions erf and erfc computed with arbitrary precision" by Chevillard in 2009. It was easier to read, but it would still take me some time to extract the answer; my first impression was $\tilde{O}(t^{3/2})$. But again, surely this question was not first investigated in 2009, was it?!

(By the way, question #5 is the one for which I really want to know the answer, but I can probably work it out from the answer to question #4.)

Answer 1

For state-of-the-art arithmetic algorithms, I'd recommend this book (a work in progress), available online, from Brent and Zimmermann: http://www.loria.fr/~zimmerma/mca/pub226.html See chapter 4.

As Steve points out, log, exp, and trig functions are $O(M(n) \log n)$ (in fact they're all calculated from log), where $M(n)$ is the cost of multiplication (theoretically $O(n\log n 2^{\log^*n})$ by Furer's algorithm) and $n$ is the number of bits accuracy. Pi also falls into this complexity. This is just theoretically, however; for less than a billion digits other algorithms with worse asymptotics are faster.

Erf is apparently harder, the book gives some algorithms based on power series and continued fractions, but it avoids giving an explicit computational cost as there are different convergence speeds in different regions.

Answer 2

$\pi$ can be computed with the hexadecimal BBP series, though apparently there are faster known ways to compute all of the bits to some level.

Knuth attributes to Brent JACM 23, 242 (1976) the result that $\log$, $\exp$, and $\arctan$ can be computed to $n$ significant bits in $O(M(n)\log n)$ where $M(n)$ is the complexity of multiplication for $n$-bit numbers. For more recent information the citations of this are probably a good bet.

Answer 3

I'll let others reply with the theoretical results.

In practice however, production systems (like Maple, Mathematica, GMP, etc) all use polyalgorithms (sometimes called hybrid algorithms) which are a collection of methods. Taking erf as an example, for real input, asymptotic series are used for large (in absolute value) inputs, and a Chebyshev-Pade approximant for small values. For complex input, the space is similarly divided up in regions, and a variety of methods can be used -- Laurent series, Puiseux series, Chebyshev-Pade approximants, continued fractions, etc. Range-reduction and reflection techniques are also quite common. So the answer for "What is the time complexity of XXX?" is really "it depends". In other words, in practice, production systems can frequently compute answers that seem to defy the theoretical results because the theoretical results are based on uniform algorithms while the practical implementations are based on adaptive algorithms.

But it looks like some very recent work on numerical computations for D-finite functions may well turn all of that on its head. See NumGfun: a Package for Numerical and Analytic Computation with D-finite Functions by Marc Mezzarobba [the full code is available as part of the gfun library].

Answer 4

Take a look at Ker-I Ko's book or his survey. I don't have the book with me right now. You can find the following in Weirauch's book:

"On compact subspaces of their domains real functions exp, sin, cos, tan, and their inverse can be computed in time $t_m(k)\log(k)$."

Where $t_m(k)$ is the time needed to compute the first k digits of the multiplication (of two real numbers) which is $O(k^2)$. This is from page 229. In the introduction he mentions that multiplication, $\sin$, $\exp$, and $\log$ can be computed in time $O(k^2)$.

--

Ker-I Ko, "Polynomial-time computability in analysis", in "Handbook of Recursive Mathematics," Volume 2, Recursive Algebra, Analysis and Combinatorics, Yu. L. Ershov et al., Eds., 1998, pp. 1271-1317. http://www.cs.sunysb.edu/~keriko/survey3.ps

Ker-I Ko, "Computational Complexity of Real Functions", Birkhauser, 1991

Klaus Weihrauch, "Computable Analysis: An Introduction", (Texts in Theoretical Computer Science. An EATCS Series), Springer, 2000

Answer 5

Try chapter 5 of Elementary Functions, Algorithms and Implementation by Jean-Michel Muller.

Answer 6

Knuth gives several algorithms for sampling from a normal distribution using elementary functions, given samples from a uniform distribution. As above, the expected complexity is $O(M(n) \log{n})$. Chapter 7 of Numerical Recipes is a brief survey of related methods. It might help to consult MPFR, since they cannot assume $O(1)$ floating-point operations. GSL has various implementations.

Answer 7

If you assume that the cost of basic arithmetic operations is independent of precision, which, of course, is not true and you know in advance the level of precision desired and you pre-compute some tables. You can calculate sin(x) and cos(x) in linear time with respect to precision using the CORDIC algorithm. I.e. After doing some relatively cheap up-front work it is possible to calculate sin(x) and cos(x) for different values of x to the same precision very cheaply.