Implementing Cosine in C from Scratch (2020)
austinhenley.com
austinhenley.com
What does "tail" mean in this context?
[1] https://en.wikipedia.org/wiki/Quadruple-precision_floating-p...
[2] https://github.com/ifduyue/musl/blob/6d8a515796270eb6cec8a27...
Also, in case of confusion, I was specifically commenting on the function over the [-pi/4, pi/4] domain in https://github.com/ifduyue/musl/blob/master/src/math/__cos.c , which the comment in https://news.ycombinator.com/item?id=30846546 was presumably about.
this is an amusing way to describe the precision of sub-normal floating point numbers
Do people do that in practice? It's on the FPU, which is basically legacy emulated these days, and it's inaccurate.
you'll pry my long doubles from my cold, dead hands!
The thing is, my first programming language was x86 assembler and the fpu was the funniest part. Spent weeks as a teenager writing almost pure 8087 code. I have a lot of emotional investment in that tiny rolling stack of extended precision floats.
You can check it out: https://github.com/jeremysalwen/vectrig
I compared several different methods of generating polynomials of different sizes for speed and precision (spoilers: taylor series were the worst and minimax polynomials (Remez algorithm) were the best).
Another (surprising) thing which I learned during the project was that the range reduction was just as (if not more) important to the accuracy of the implementation than the polynomial. If you think about it, you will realize that it's actually pretty difficult to quickly and accurately compute the sin of large numbers like 2^50.
I also tried to directly optimize the coefficients for the accuracy of the polynomial on the required range, but that experiment was unsuccessful.
It's all there in the repository, the implementations, notes about the different polynomials used, and the accuracy/speed statistics for the different methods.
I would have expected at least an LSQ approximation with a basis of Legendre polynomials thrown into the mix. I got that as a basic homework in my numerics class once, after we've shown to ourselves in the class that [1, x, x², x³...] is not a really good basis to project things onto.
Some of the performance claims have a caveat, though: Lookup tables and micro-benchmarks don't mix well. At all.
I just added a simple loop that thrashes the L3 cache every 500 iterations (diff here: https://goonlinetools.com/snapshot/code/#sm4fjqtvjyn36dyednc...). Now the method recommended at the end (cos_table 0_001 LERP) is slower than glibc's cos() (while still having an accuracy that is more than 10^8 times worse)!
Time benchmark output:
cos_table_1_LERP 0.0038976235167879
cos_table_0_1_LERP 0.0042602838286585
cos_table_0_01_LERP 0.0048867938030232
cos_table_0_001_LERP 0.0091254562801794
cos_table_0_0001_LERP 0.0139627164397000
cos_math_h 0.0089332715693581
cos_table_1_LERP 64 bytes
cos_table_0_1_LERP 512 bytes
cos_table_0_01_LERP 5040 bytes
cos_table_0_001_LERP 50280 bytes
cos_table_0_0001_LERP 502664 bytesI think that the grandparent concerns are very legitimate.
For instance, in a game, it's common to have a gun/spell or whatever that shoots enemies in a cone shape. Like a shotgun or burning hands. One way to code this is to calculate the angle between where you're pointing and the enemies, (using arccos) and if that angle is small enough, apply the damage or whatever.
A better way to do it to take the vector where the shotgun is pointing, and the vector to a candidate enemy is, and take the dot product of those two. Pre-compute the cosine of the angle of effect, and compare the dot product to that -- if the dot product is higher, it's a hit, if it's lower, it's a miss. (you can get away with very rough normalization in a game, for instance using the rqsrt instruction in sse2) You've taken a potentially slow arccos operation and turned it into a fast handful of basic float operations.
Or for instance, you're moving something in a circle over time. You might have something like
float angle = 0
for (...)
angle += 0.01
<something with sin(angle) cos(angle)>
const float sin_dx = sin(0.01)
const float cos_dx = cos(0.01)
float sin_angle = 0
float cos_angle = 1
for (...)
const float new_sin_angle = sin_angle * cos_dx + cos_angle * sin_dx
cos_angle = cos_angle * cos_dx - sin_angle * sin_dx
sin_angle = new_sin_angle
And there's the ever popular comparing squares of distances instead of comparing distances, saving yourself a square root.
In general, inner loops should never have a transcendental or exact square root in them. If you think you need it, there's almost always a way to hoist from an inner loop out into an outer loop.
Invisible sprite:
[ when I receive "sine" ]:
[ go to center ]
[ turn A degrees ]
[ move 1 step forward ]
For small angles, sin(x) ~= x and cos(x) ~= 1 -- the ["small angle approximation"](https://en.wikipedia.org/wiki/Small-angle_approximation)
It's actually kind of ridiculous how many times it comes up and works well enough in undergraduate Mechanics.
"Relaxation" algorithms for calculating fields work that way.
https://en.wikipedia.org/wiki/Small-angle_approximation#Erro...
This comes up in movement under constant rate of rotation as well as Fourier transforms. See [1] for details, but the basic idea uses the simple trig identities:
cos(theta + delta) = cos(theta) - (a*cos(theta) + b*sin(theta))
sin(theta + delta) = sin(theta) - (a*sin(theta) - b*cos(theta))
a = 2 sin^2(delta/2)
b = sin(delta)
[1] Press, William H. et. al., Numerical Recipes, Third Edition, Cambridge University Press, p. 219
It turns out the BIOS on the hardware can do sine/cosine, but you have to do that via inline assembly (swi calls) and it's no faster than the simple lookup table.
This is all just to get terms for the affine matrix when rotating sprites/backgrounds.
I've also read that Chebyshev polynomials are far better for approximating functions than Taylor series (eg https://en.wikipedia.org/wiki/Approximation_theory)
Alternative avenues that complement the approaches shown:
- Padé Approximants (https://en.wikipedia.org/wiki/Pad%C3%A9_approximant) can be better than long Taylor Series for this kind of thing.
- Bhaskara I's sin approximation (https://en.wikipedia.org/wiki/Bhaskara_I%27s_sine_approximat...) is easily adaptable to cosine, remarkably accurate for its simplicity and also fast to calculate.
Do you have any suggested reading?
If you’re measuring “how good is my cosine function” by calculating the maximum error, then it makes sense to use an optimization algorithm that minimizes that error.
The Remez algorithm is rather elegant. The basic idea is that you come up with a polynomial, figure out where the local error maximums are, and then calculate a new polynomial with better local error maximums in the same locations. Repeat. With the right conditions it converges quickly and to a global minimum.
For example, let’s say that you try to come up with a quarter-wave cosine approximation. You end up with a polynomial that has local error maximums at x=[0, 0.1, 0.3, 0.7, 1.4, 1.57], and various values of ∆y. You create a new polynomial and design it to have exact errors of ±∆y at the same x locations… the errors will alternate in sign as the polynomial goes above and below the target function. However, when you design the function, you’re just doing a polynomial fit through these (x,y) coordinates, so the new polynomial will have local error maximums at different x coordinates. Note that the boundaries are also error maximums. You end up with exactly the right number of degrees of freedom to solve the equation.
Each time through the loop you get a new set of x coordinates. Under the right conditions, these converge. Eventually the ∆y errors are all nearly equal to ± some global error maximum, with alternating signs.
I used this to approximate sine and exp functions for an audio synthesis tool I’m working on—the nice thing about polynomials over lookup functions is that they are easier to convert to SIMD.
Here is the NumPy code I used to find the coefficients. Note that this is not taken from any kind of numerical recipes cookbook and I can’t really vouch for the algorithm I used from a numerical standpoint—I only checked that it works empirically.
https://github.com/depp/ultrafxr/blob/master/math/coeffs/cal...
https://en.wikipedia.org/wiki/Bhaskara_I%27s_sine_approximat...
https://en.wikipedia.org/wiki/Archimedean_spiral
(This was on a MSP430 platform with a FPU that only did multiplication.)
Even a 50% increase in speed was well worth my effort though.
* Replacement sqrtf() function. This function runs about 50% faster than
* the one in the TI math library. It provides identical accuracy, but
* it is limited to input values less than 2^14. This shouldn't be a problem
* because the passed argument will always be in the range of
* 0-(max spiral scan period). The spiral scan period is presently about
* 140 seconds.
[ x(n + 1) ] = [cos(ω), -sin(ω)] [ x(n) ]
[ y(n + 1) ] [sin(ω), cos(ω)] [ y(n) ]
where
x(n) = cos(ω n)
y(n) = sin(ω n)
x(0) = 1
y(0) = 0
ω = frequency (radians) * time_step
It's not suitable for numerous cases (notably if the frequency needs to change in real time) but provided that the sum of the coefficients sum to `1.0` it's guaranteed stable and will not alias.
Not sure who initially came up with it.
[1] https://github.com/robrohan/wefx/blob/1a918cc2d5ad87402a3830...
* The range reduction is, well, probably worse than acos. Look for Payne–Hanek–Corbett. Not an issue for games, I assume.
* The LERP is nice and all, but Qt actually has an even better one: instead of doing a LERP, you use the local derivative also available from the lookup table. https://stackoverflow.com/a/52841086
Others have already mentioned the Remez and the Horner stuff. Won’t repeat here.
The author's compiler must not be great at optimization if this is the case.
Some compilers have options such as -ffast-math that result in ignoring some precision issues. This could be worth a try as an alternative to manual optimisation.
cordic beautifully explained; i have used this approach on fixed point dsp processors in another life; its simplicity, and elegance is beautiful (on controllers with no fp)
CORDIC is still very common in embedded systems, especially with CPUs that have no hardware multiply/divide.
The author missed one nice table method. When you need sin() and cos() split the angle into a coarse and fine angle. Coarse might be 1/256 of a circle and fine would be 1/65536 of a circle (whatever that angle is). You look up the sin/cos of the coarse angle in a 256 entry table. Then you rotate that by the fine angle which uses another 256 sin * cos entries. The rotation between coarse values is more accurate than linear interpolation and give almost the same accuracy as a 65536 entry table using only 768 entries. You can use larger tables or even another rotation by finer angle steps to get more accuracy.
Explain?
Are you trying to reduce the error in the linear interpolation by describing the convex curve between the endpoints?
The derivative of the cosine is the sine, and the derivative of the sine is the cosine...so I'd expect the required output of the fine angle table to interpolate 0.15 between 0.1 and 0.2 and the required output to interpolate 45.15 between 45.1 and 45.2 degrees will be very different.
Yeah I thought that might not be clear enough. Let's say you want 0.1 degree accuracy but don't want a 3600 entry table. Make one table for every 1 degree. You get to use the same table for sin and cos by changing the index. That is the coarse table with 360 entries. Then make another table with sin and cos values for every 0.1 degrees, or specifically 0.1*n for n going from 0-9. This is another 20 values, 10 sin and 10 cos for small angles.
Take your input angle and split it into a coarse (integer degrees) and fine (fractional degrees) angle. Now take the (sin(x),cos(x)) from the coarse table as a vector and rotate it using the sin & cos values of the fine angle using a standard 2x2 rotation matrix.
You can size these tables however you need. I would not use 360 and 10 for their sizes, but maybe 256 and 256. This can also be repeated with a 3rd table for "extra fine" angles.
That's probably not going to perform well for an operation this fast. The computation is faster than main memory read.
A common trick to make that faster is to not use doubles and radians to represent the phase. Instead, represent phase using an integer in some power-of-two range. That lets you truncate the phase to fit in a single period using a bit mask instead of the relatively slow modulo.
Doom?
You could go further and only to 0 to pi/2 as it is mirrored and flipped for pi/2 to pi.
Or you could go even further and do only 0 to pi/4 and use a simple trig identity and your presumably parallely-developed sin function for the cos of pi/4 to pi.
> One of the implementations is nearly 3x faster than math.h
Is this true even for -O3?