How to calculate a "rational" sine/cosine to specific digits of precision?

Question

I have a toy calculator I've been tasked to create for a class project. That's great - I have most of the normal calculator stuff done (order of operations, variables, precision settings, predefined functions). The issue is that, since I'm using rational arithmetic, there's no sin/cos function built-in, so I'll have to make my own.

And the issue with that is that sines are usually irrational, so I'll have to limit it to some precision or it'll run forever. I've successfully implemented the sine Maclaurin series in code. For reference, the relevant equation from the Wikipedia article is this:

$$\sin x = \sum_{n=0}^{\infty}{\frac{(-1)^n}{(2n + 1)!}x^{2n+1}} = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} - \cdots$$

My program runs only until this condition is true:

$$\frac{x^i}{i!}<\frac{1}{10^p}$$

where $i = 2n + 1$ and $p$ is the desired precision in digits, but it gets very slow quite quickly (at least using a C++ interpreter, but I want it to run quickly there before I make it faster by compiling). There are some optimizations such as keeping the last numerator/denominator, but it's still not that fast.

Calculating the sine of $5\pi$ to 15 digits of precision (bad example, but it's to prove a point) using my program took 5 seconds across 34 iterations, since many later iterations end up multiplying integers greater than $1 \times 10^{900}$. $\pi$ has to be represented as a fraction, so there were very large integers involved.

I'm not sure where to ask for a better algorithm but here (or Stack Overflow, from which I've been banned for over a year, so that's not going to happen...). Even though I'm using rational arithmetic, I'm wondering if there's a more computationally efficient way to approximate a $\sin$ and $\cos$ to some amount of digits of precision. I know this isn't a programming site, but surely someone can help with the algorithm side? (preferably explained in plain English)

Note: I could offload the work onto some arbitrary-precision floating point library and then convert that back to a "rational", but that seems like cheating to me... it'll be more of a last resort.

Answer 1

Probably a typo since $$\sin (x) = \sum_{n=0}^{\infty}{\frac{(-1)^n}{(2n + 1)!}x^{2n\color{red}{+}1}} = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} - \cdots$$

Then, for $|x| > \pi$, I should start adding or subtracting as many $\pi$ as possible writing $x=y+k \pi$ in order to stay in the range $0 \leq y \leq \pi$.

Finally, I would avoid any test (this is time consuming) and determine in advance the number of terms to be used for the required accuracy.

You want to approximate

$$\sin(y)= \sum_{n=0}^{p}{\frac{(-1)^n}{(2n + 1)!}y^{2n+1}}$$ such that $$\frac{y^{2p+3}}{(2p+3)!} \leq 10^{-k}\implies(2p+3)! \geq y^{2p+3}\, 10^k$$

If you look at this question of mine, you will see a magnificent approximation proposed by @robjohn. $$m!=y^m\, 10^k\implies m\sim ey\exp\left(\operatorname{W}\left(\frac k{ey}\log(10)-\frac1{2ey}\log(2\pi y)\right)\right)-\frac12$$ where $W(.)$ is Lambert function.

Just make $m=2p+3$ and get $p$.

Let us try with $x=123.456$ and $k=15$; we then need to subtract $39\pi$ to get $y=0.933887$. Using the formula, this will give $p=7$. You can easily check that this is correct.

Answer 2

I would suggest the following tips:

1. Since you are using rational arithmetic and you need to limit the accuracy because sine is typically irrational, the easiest way to do so is to work consistently with fractions of the form $\frac{a}{B^n}$ where $a$ is integer (of course), and $B$ is some suitable base - i.e. terminating base-$B$ fractions, and you limit $n$ as needed. One might want to choose $B = 10$, but for computer arithmetic $B = 2$ is probably best since then you can maximize efficiency using bit shifts to do a fair bit of the fraction arithmetic. (In that case, you will probably want to limit $n$ to 53, which is the number of bits needed to get around $10^{-16}$ residual as in IEEE 754 double precision floating point arithmetic.) You can truncate this representation as desired simply by shifting the numerator and denominator. Just don't reduce the fractions the usual way until perhaps the very end (if you're "allowed" to do that, i.e. how much control do you have over the rational arithmetic package?). Basically you'll be implementing fixed/floating (here really, fixed) point arithmetic with a fraction representation.

2. The second bit I'd want to point out is the Taylor series for sine really only works best when the input is in $(-1, 1)$. That is because outside this range, the "power of $x$" part grows up, while within it, it shrinks down. When it grows up, it's in effect "fighting" the factorials in the denominator which are trying to shrink the terms. Since the factorial eventually "wins" for suitably large term index, the series still converges, but you lose precision if you're doing this in floating-point, or will rack up huge fractions in rational, and moreover will have to sum a great deal of terms (in fact FP may not converge at all due to precision loss if the $x$ is large enough). Within this interval, however, because the power part is shrinking down, it actually works with the factorial and the series converges even faster than that for $e$ - that is to say, pretty damn fast. Taking it a bit more strictly as $(-\frac{1}{2}, \frac{1}{2})$ would be ideal and relatively simple. To ensure the argument is within this range, you can use trigonometric reduction identities like $$\sin\left(\frac{\theta}{2}\right) = \pm \sqrt{\frac{1 - \cos(\theta)}{2}}$$ $$\cos\left(\frac{\theta}{2}\right) = \pm \sqrt{\frac{1 + \cos(\theta)}{2}}$$ where you will need to mind your quadrant for the sign of the square root. For reduction of a much large angle, subtracting/adding $2\pi$ (i.e. $\tau$) repeatedly is possible, but actually just halving all the way down has the better asymptotic complexity.