6
$\begingroup$

I am writing a rational number library for performing fast math on integer-only microcontrollers. So far, I have all the basic operations, as well as a first-step estimation function for square root. It takes a positive integer input and returns a rational number approximating the square root (accurate to roughly 4 bits of precision; to be refined further by Newton-Raphson or similar method). Obviously, inverse (reciprocal) square root comes for free since the answer is a fraction that can be trivially inverted.

I read somewhere that other functions (e.g. trigonometric functions) can be calculated in terms of inverse square root (can't find the link now), but I am having a hard time finding any algorithms since search results are almost entirely about inverse trigonometric functions (mostly high school help stuff).

Does anyone have a handy resource where I can find efficient methods for calculating sin, cos, tan, atan, etc. in terms of square root or inverse square root? Is it possible or am I imagining things?

I know this is not a "code" StackExchange but I believe this is fundamentally a math question. If you're interested in my code, it's here. The sqrt_est() function is somewhat remarkable in that it is very short, uses only shifts and adds, no looping or iteration, or even multiplies or divides! (core implementation is on lines 39-45.)

$\endgroup$
3
  • $\begingroup$ I have used (and advocated for) inverse square roots to eliminate trig functions when you already have another related trig function. But not to compute trig functions of angles from scratch. Is that what you remember? $\endgroup$
    – Dan Piponi
    Jan 31, 2021 at 3:56
  • $\begingroup$ @DanPiponi perhaps... I could probably implement at least one reasonably fast trig function. Any pointers to examples of what you're describing would be appreciated! $\endgroup$ Jan 31, 2021 at 4:57
  • $\begingroup$ To satisfy my curiosity, I am going to reveal my ignorance. Consider as a starting point using Taylor Series, where you first convert $x$ to the interval $(-\pi,\pi].$ Are you saying that the corresponding computational performance here can be improved? $\endgroup$ Jan 31, 2021 at 5:52

2 Answers 2

2
$\begingroup$

I personally am not aware of any method utilizing the inverse square root for evaluating trig functions and without a reference, it's really hard to understand what you're referring to. With such implementations, the general idea is to go with ready-made, optimized and tested implementations, but that probably won't answer your question.

Taylor Series is certainly a good starting point and although it technically has a time complexity of $O(n)$ it converges really fast (you can achieve really good precision by calculating a few terms) so depending on the scope of your library it may actually be ok. Chebyshev approximation is usually prefered since it converges faster and I haven't seen many real-world applications favouring Taylor expansions over it.

Another "dirtier" option is to simply use a look-up table and either do closest-match or interpolate, losing some precision along the way. You don't need to create a particular big look-up table either, just create one for many steps along the interval $[0, x]$ and then partition $[0, \pi]$ with a step of x. Then you can calculate $\sin(y)$ and $\cos(y)$ by writing them as:

$$\sin(y_1 + \delta_1) = \sin(y_1)\cos(\delta_1) + \cos(y_1)\sin(\delta_1)$$

$$\cos(y_2 + \delta_2) = \cos(y_2)\cos(\delta_2) + \sin(y_2)\sin(\delta_2)$$

where $\delta_1,\ \delta_2 < x$.

For $\arctan$ the basic idea is quite similar. You again usually end up partitioning an interval $[0, N]$ (for $x>N$ you assume the answer to simply be $\frac{\pi}{2}$ or use some rough and cheap approximation method) and approximate your function in each part using something like a Chebyshev series. You again use trig identities to try to bring your argument down to $[0, N]$ and evaluate according to your polynomial in that range. The specifics of this can be quite tricky without an FPU so I don't know if you would want to bother with that.

An alternative solution, is to try to solve the integrals in each interval instead of evaluating the series. For example, we know that $\arctan(x) = \int_0^x\frac{1}{1+x^2}$. You can again create a LUT and then either use a numerical integration method (like Newton-Cotes) which would require several evaluations of $\frac{1}{1+x^2}$ depending on the target precision, which I'm not sure you can afford. You could also utilize the much simpler to calculate Simpson's Rule. I haven't seen those methods used though, so it's probably not worth it in terms of efficiency, ease of implementation and precision. Just throwing this in in case it might be worth considering.

Another cheap way is using the CORDIC algorithm which intuitively works by doing successive rotations around the imaginary unit circle. You can find a surprisingly in-depth explanation in Matlab's documentation

Of course this is all math. If you asked in StackOverflow they would let you know about a myriad other things that factor in your implementation like cache misses, utilizing SIMD intrinsic and other stuff this answer didn't touch on ;)

$\endgroup$
1
  • $\begingroup$ Thanks for the variety of information. Yes, at this point I am more interested in the mathematical options than the specifics of machine implementation. If you look at my code for sqrt_est() you can see it’s already remarkably optimal, which is why I’d like to leverage it as much as possible in computing other functions. What I am really looking for is what @Dan Piponi mentioned in his first comment on my original question. $\endgroup$ Jan 31, 2021 at 15:16
1
$\begingroup$

Square roots are key to all sorts of computations because they let us compute arbitrary powers. If you want to compute $2^{1/7}$, for example, first write $\frac17$ in binary as $0.001001001001001\dots$, which tells us that $$ 2^{1/7} = 2^{1/8} \cdot 2^{1/64} \cdot 2^{1/512} \cdot 2^{1/4096} \cdot 2^{1/32768} \cdot \dots $$ All of these are repeated square roots: $2^{1/2} = \sqrt 2$, $2^{1/4} = \sqrt{2^{1/2}}$, $2^{1/8} = \sqrt{2^{1/4}}$, and so on. So we can take as many square roots as we need to, and multiply them together.


Square roots are key to trig functions because those are just powers. If you want to compute $\cos \frac\pi7$, for example, then it's enough to compute $(-1)^{1/7} = \cos \frac\pi7 + i \sin \frac\pi7$. This is done in the same way as above: $$ (-1)^{\frac17} = (-1)^{1/8} \cdot (-1)^{1/64} \cdot (-1)^{1/512} \cdot (-1)^{1/4096} \cdot (-1)^{1/32768} \cdot \dots $$ Computing iterated square roots $(-1)^{1/2} = \sqrt{-1}$, $(-1)^{1/4} = \sqrt{(-1)^{1/2}}$, and so on is a bit trickier, but comes down to real square roots in the end.


But if the $(-1)^x$ argument seems a little bit disingenuous to you, then we can rephrase this in terms of trig identities. Start with a known $(\cos x, \sin x)$ pair: it's convenient to start with $(\cos \frac\pi2, \sin \frac\pi2) = (0,1)$, but we can start anywhere if we want to. From $(\cos x, \sin x)$ we can compute $(\cos \frac x2, \sin \frac x2)$ by the pair of identities $$ \cos \frac x2 = \pm \sqrt{\frac{1 + \cos x}{2}}, \qquad \sin \frac x2 = \pm \sqrt{\frac{1 - \cos x}{2}}. $$ The angles are determined by the quadrant, which in this case is always going to be positive for both $\cos$ and $\sin$. This lets us compute $(\cos \frac{x}{2^n}, \sin \frac{x}{2^n})$ for many $n \ge 0$, which we can actually save in a table.

Then, if you want to compute $\cos \theta$ for an arbitrary angle $\theta \in [0,2x]$, where $x$ is our reference angle, first write $\frac{\theta}{x}$ in binary. This is telling us $\theta$ as a sum of some subset of $\{x, \frac x2, \frac x4, \frac x8, \frac x{16}, \frac x{32}, \dots\}$, all of whose sines and cosines we know. Angle addition formulas tell us $\cos \theta$ and $\sin \theta$.

This is very similar in flavor to CORDIC algorithm also mentioned in the other answer, but that one is more efficient (and also more complicated).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.