How to find the value of a function correct to $n$ decimal places as efficiently as possible?

Question

How to find the value of a function correct to $n$ decimal places as efficiently as possible?

For example, how to find $\sin(1)$ correct to six decimal places as efficiently as possible?

I computed $\sin(1)$ correct to six decimal places as follows:

$$|\sin(1) - (1 - \frac{1}{3!} + \frac{1}{5!} - \frac{1}{7!} + \frac{1}{9!} - \frac{1}{11!})| < \frac{1}{(2 \times 6 + 1)!} \fallingdotseq 1.6059 \times10^{-10}.$$
$$1 = 1.000000000.$$ $$-\frac{1}{3!} \fallingdotseq -0.166666667.$$ $$\frac{1}{5!} \fallingdotseq 0.008333333.$$ $$-\frac{1}{7!} \fallingdotseq -0.000198413.$$ $$\frac{1}{9!} \fallingdotseq 0.000002756.$$ $$-\frac{1}{11!} \fallingdotseq -0.000000025.$$ $$1 - \frac{1}{3!} + \frac{1}{5!} - \frac{1}{7!} + \frac{1}{9!} - \frac{1}{11!} \fallingdotseq 0.841470984.$$
$$|0.841470984 - (1 - \frac{1}{3!} + \frac{1}{5!} - \frac{1}{7!} + \frac{1}{9!} - \frac{1}{11!})| \leq 6 \times (5 \times 10^{-10}) = 30 \times 10^{-10}.$$

So,

$$|\sin(1) - 0.841470984| < 1.6059 \times10^{-10} + 30 \times 10^{-10} < 100 \times 10^{-10} = 10^{-8}.$$

So,

$$0.841470974 < \sin(1) < 0.841470994.$$

So, the answer is $\sin(1) \fallingdotseq 0.841470$.

By the way, how to find the value of $1 + \frac{1}{2} + \frac{1}{4} + \frac{1}{8} + \frac{1}{16} + \cdots$ correct to six decimal places in a similar way?

I guess we cannot compute the value of $1 + \frac{1}{2} + \frac{1}{4} + \frac{1}{8} + \frac{1}{16} + \cdots$ correct to six decimal places in a similar way.

Answer 1

I guess that you relied on the fact that you have an alternating series with decreasing terms. But

• these is no guarantee that stopping at the first term smaller than the desired accuracy is optimal (by the way, in your first example you can stop at $9!$),

• you can artificially turn a non-alternating series to an alternating one, for instance

$$1+\frac12+\frac14+\cdots\frac1{2^n}+\cdots=4-3+2-\frac32+1-\frac34+\frac12-\frac38+\cdots+\frac4{2^n}-\frac3{2^n}+\cdots$$

---

In general, you have two possible approaches:

• use an analytical expression of the remainder of the series. In the case of Taylor series, several formulas are available. https://en.wikipedia.org/wiki/Taylor%27s_theorem#Explicit_formulas_for_the_remainder

• find an upper bound on the tail of the summation by gross approximations of the terms, e.g.

$$\frac1{5!}+\frac1{6!}+\frac1{7!}+\cdots<\frac1{5!}+\frac1{5\cdot5!}+\frac1{5^2\cdot5!}+\cdots=\frac1{96}.$$

You can refine these expressions to get tighter bounds, but at some point, the cost of evaluating the remainder will exceed the cost of evaluating more terms of the series.

---

In practice, the important functions are evaluated in chosen intervals by means of precomputed approximation polynomials, which are optimal in a minmax or least-squares sense, with a chosen global accuracy.

Answer 2

There is no general answer. It depends very much on what you want to calculate and what tools you have available. Your calculation of $\sin(1)$ is a good one-it is justified by the alternating series theorem and is not too hard. It is basically using the Taylor series around $0$. The natural improvement would be to use the Taylor series around some closer point. That would take fewer terms, but they will be harder to calculate. You could say $\sin (1)=\cos(\frac \pi 2 -1)$ so $$\sin(1)\approx1-\frac1{2!}(\frac \pi 2-1)^2+\frac1{4!}(\frac \pi 2-1)^4-\frac 1{6!}(\frac \pi 2-1)^6$$ where the error bound is about $3\cdot 10^{-7}$. This has fewer terms than your expression, but each one is more complicated. Your expression was simplified by the fact that your argument was $1$. It would not be so nice if you asked for $\sin(.99343424)$. In that case my version would win easily.

Another subtlety is when you ask for correct to six places do you mean within $\pm 5\cdot 10^{-7}$ or do you mean you write down six decimal places that are correct. In the second case you might (but are unlikely to) be very close to a point where the rounding flips from one value to another and then need an extremely accurate value to get six good places.

For series like your $1 + \frac{1}{2} + \frac{1}{4} + \frac{1}{8} + \frac{1}{16} + \cdots$ where all the terms are positive you need to find another way of bounding the error. This one we sum analytically and get $2$.

Answer 3

As said, there is no general way. However, in some cases such as the problem of the sine function, there is an elegant solution.

We have $$\sin(x)=\sum_{n=0}^\infty \frac{(-1)^n}{(2n+1)!} x^{2n+1}=\sum_{n=0}^m \frac{(-1)^n}{(2n+1)!} x^{2n+1}+\sum_{n=m+1}^\infty \frac{(-1)^n}{(2n+1)!} x^{2n+1}$$ and, since it is an alternating series, you look for $m$ such that $$\frac {x^{2m+3}}{(2m+3)!} \leq 10^{-k}\implies (2m+3)! \geq x^{2m+3}\, 10^k$$ For the time being, let $p=2m+3$ and consider that we look first for the solution of $p!=x^p\, 10^k$.

If you look at this question of mine, you will find a magnificent approximation proposed by @robjohn, an eminent MSE user

$$p\sim ex\exp\left(W\left(\frac{1}{2 e x}\log \left(\frac{10^{2 k}}{2\pi x}\right)\right)\right)-\frac12$$

where appears Lambert function.

From this, we have for the sine function $$\color{blue}{m\sim \frac { e x} 2\exp\left(W\left(\frac{1}{2 e x}\log \left(\frac{10^{2 k}}{2\pi x}\right)\right)\right)-\frac74}$$ and for sure, you need to use $\lceil m \rceil$.

Applied to the case $k=6$ and $x=1$, this would give $m=3.22$ that is to say $m=4$. Checking $$\frac 1 {11!}=2.50 \times 10^{-8}\quad \text{(OK)} \qquad \text{while}\qquad \frac 1 {8!}=2.48 \times 10^{-5}\quad\text{(not OK)}$$ and $m=4$ was your answer.

To give you an idea about the accuracy, considering $m$ as a real, the exact solution of the equation would be $m=3.22280$ while @robjohn's approximation gives $m=3.22189$.

For a shortcut evaluation of $W(t)$ when $t$ is large, remember the approximation (have a look here) $$W(t) \sim L_1-L_2+\frac{L_2}{L_1}+\frac{L_2(L_2-2)}{2L_1^2}+\frac{L_2(6-9L_2+2L_2^2)}{6L_1^3}+\cdots$$ where $L_1=\log(t)$ and $L_2=\log(L_1)$. Applied to your case, this would give $W(t)=1.30344$ while the exact value would be $W(t)=1.29695$.