10 ms·

Implementing Cosine in C from Scratch (2020)

220 points by alraj 3y ago

azhenley 3y ago
No need for archive.org, my website moved a few years ago: https://austinhenley.com/blog/cosine.html https://austinhenley.com/blog/cosine.html
- butterlesstoast 3y ago
  Thanks for taking the time to make such an incredible article. Made my day reading it this morning.
- dang 3y ago
  Fixed. Thanks!
- capitainenemo 3y ago
  BTW, WRT projects using a lookup table... https://hg.hedgewars.org/hedgewars/file/tip/hedgewars/uSinTable.pas https://hg.hedgewars.org/hedgewars/file/tip/hedgewars/uSinTa... I think there's even a github mirror someone made since that's where you were looking. Maybe a little out of date, but this is old code.
  - capitainenemo 3y ago
    oh... also, sine table isn't just a performance thing. Hedgewars does fixed point custom math for deterministic lockstep, so keeps the math simple and predictable on various compilers/hardware.
    - nitwit005 3y ago
      That's something I'd wondered about with games that require determinism. Seems fairly challenging to deal with for that sort of game.
      - capitainenemo 3y ago
        https://nicolas.brodu.net/en/programmation/streflop/index.html https://nicolas.brodu.net/en/programmation/streflop/index.ht... is apparently pretty reliable, but yeah, it's a big problem. Spring and 0AD both struggled with it too.
  - deleted 3y ago
    [deleted]
Helenarttr 3y ago
That was so amazing Information. https://www.ballsportspro.com/how-does-a-pickleball-ladder-work/ https://www.ballsportspro.com/how-does-a-pickleball-ladder-w...
xipix 3y ago
Nice. I'd love to see how this changes when you have SIMD and multiple cosines to compute in parallel. Also when you have to compute sine and cosine simultaneously which is often the case, and then you may be more interested in polar error than cartesian error.
- p0nce 3y ago
  You can do this with: https://github.com/RJVB/sse_mathfun https://github.com/RJVB/sse_mathfun Those functions do not have the same safety as libc, but even for scalar sin/cos it can be faster.
sampo 3y ago
> cosine repeats every 2pi > We can do better since the cosine values are equivalent every multiple of pi, except that the sign flips. There is one more step to take: Each π/2 long segment has identical shape, they are just pointing up or down (like you already noticed), and left or right. So you can reduce you basic domain down to not just 0...π, but to 0...π/2.
- amelius 3y ago
  But that takes extra clock cycles, since now you have to fiddle with the sign.
  - Someone 3y ago
    I think it’s likely that you get that more than back by having to do (at least) one iteration fewer. Those who know more about this: educate me.
    - adgjlsfhk1 3y ago
      Pretty much every implementation does reduction mod pi/2. Switching polynomials is basically free (and halving the domain significantly reduces the number of terms needed).
- IIAOPSW 3y ago
  And you can keep going with this logic, building cos from 0-x by recursively calculating cos from 0-x/2 then mirroring over some line if needed. Eventually you hit machine precision limit aka the base case cos(0) = 1.
- ajross 3y ago
  And still further: the function can be piecewise decomposed into separate polynomials, which is what the pervasive SunPro code in everyone's libc does. There's one "mostly parabolic" function to handle the curvy top, and another "mostly linear" one to handle the zero crossing region.
- gmiller123456 3y ago
jpfr 3y ago
The Taylor expansion locally fits a polynom based on the n first derivatives. If you want to find the "best" nth-degree polynom to approximate the sine function, functional analysis gives the tools for solving that optimization in closed form. By selecting an appropriate norm (in function space) you can minimize either the maximum error or the error integral over some range (e.g. the 0-\pi range). Here's a video on the subject. You might want to watch earlier ones also for more context. https://www.youtube.com/watch?v=tMlKZZf2Kac&list=PLdkTDauaUnQpzuOCZyUUZc0lxf4-PXNR5&index=28 https://www.youtube.com/watch?v=tMlKZZf2Kac&list=PLdkTDauaUn... Full disclosure, this is my university lecture on optimization that was recorded during Covid.
- stephencanon 3y ago
  It’s worth noting that for fixed-precision math libraries, we usually want the best approximation in the (possibly weighted) L-inf sense, which doesn’t generally have a closed-form solution, so we use iterative methods to generate approximations instead.
  - bigdict 3y ago
    What's a good textbook intro to these methods?
    - remcob 3y ago
      I wrote this intro for a friend and since then more people have found it useful: https://2π.com/22/approximation https://xn--2-umb.com/22/approximation If you want a proper textbook, I recommend L.N. Trefethen (2019). “Approximation Theory and Approximation Practice.” Which is available online here: https://epubs.siam.org/doi/book/10.1137/1.9781611975949 https://epubs.siam.org/doi/book/10.1137/1.9781611975949
      - stephencanon 3y ago
        That's a really nice overview of Remez, thanks for sharing!
femto 3y ago
If doing signal processing, an optimisation is to use an N-bit integer to represent the range 0 to 2.pi. Some example points in the mapping: 0->0, 2^(N-2)->pi/2, 2^(N-1)->pi, 2^N(wraps to 0)->2.pi(wraps to 0). If your lookup table has an M-bit index, the index to the lookup table is calculated with: index = (unsigned)theta>>(N-M), where theta in the N-bit integer representing the angle. The fractional part, which can be used for interpolation, is: theta & ((1<<(N-M))-1). If you choose M=N=word size (16 bits is often nice), the angle can be used directly as an index. With 16-bits accuracy is typically good enough without interpolation and the entire trig operation reduces to an array indexing operation. This minimises the overhead of converting an angle to a table index.
- rjmunro 3y ago
  Why just signal processing? Surely this is good for storing e.g. rotations of characters in games, longitude, in fact almost anywhere where you want an angle where a full rotation is the same as no rotation at all, and more precision as you approach 0° is not useful.
  - vikingerik 3y ago
    Accumulated error can matter, which is where you'd want that precision near 0°. If something is making tiny turns each frame for thousands of frames, floating-point inaccuracy can add up. (A smarter mathematical approach might remedy this, or might not; think of some game component that is supposed to jitter randomly each frame, maybe some particle system.) It's not common, but it's within hypothetical feasibility that the precision would be useful.
    - amitport 3y ago
      He is describing deterministic uniform round to nearest quantization. The most common smarter mathematical approach would be randomized rounding and it will solve this issue. (Issue being that the deterministic rounding is a biased estimator of the sin function)
      - 3y ago
Aardwolf 3y ago
" In all my time of coding, there has only been one situation in which I used cosine: games, games, games. function move(obj) { obj.x += obj.speed * Math.sin(obj.rotation); obj.y += obj.speed * Math.cos(obj.rotation); } " Why not store the velocity as a 2D vector instead? Then you still have to use cos/sin to compute this vector, but at least you don't need it every frame, plus often you don't need to use cos/sin to compute this vector either since forces that act on the velocity themselves can have an x and y component you can directly add to it
- messe 3y ago
  Sometimes it's just easier to work and think in polar coordinates. If that's the case, then (speed, rotation) is just as valid a choice as (vx, vy). For smaller games, the cost of computing sin and cos every frame is so astronomically small, that it's not even concerning yourself over.
  - Aardwolf 3y ago
    Hmm, I usually find polar coordinates harder to reason with due to the gimbal lock :) In 2D there's no gimbal lock of course, but there's still a number that behaves differently than the others, the angle wrapping around in range 0..2pi! And of course the article is about fast computation of sin/cos so the option of not computing it should be considered
    - thehappypm 3y ago
      Can you explain gimbal lock? I can’t wrap my head around it
      - Aardwolf 3y ago
        The main thing is how when representing the 3D looking orientation as multiple angles (pitch=up/down, yaw=compass direction), you lose control over the yaw axis when looking straight up.
aebtebeten 3y ago
Niklaus Wirth also has some cosine code (probably meant more for software floating point, or fp FPGA blocks); I don't know how they compare with these approximations but his seem to be within 1e-6 of python's math.cos ... https://people.inf.ethz.ch/wirth/ProjectOberon/Sources/Math.Mod.txt https://people.inf.ethz.ch/wirth/ProjectOberon/Sources/Math....
dm319 3y ago
The answer to Cos(1.57079632) (in radians) will give you a clue as to how your calculator does it. See here [0] [0] https://www.reddit.com/r/calculators/comments/126st95/cos157079632_is_a_bit_of_a_torture_test/ https://www.reddit.com/r/calculators/comments/126st95/cos157...
- the_third_wave 3y ago
  Quite well when you realise this comes down to cos(π/2) to which the thing answers correctly: cos(π/2) = 0. Using HiPER Calc on Android
  - dm319 3y ago
    The number is deliberately close but not the same as pi/2. See here: https://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv017.cgi?read=123814 https://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv017.c... And you can compare to wolfram alphas result.
- aidenn0 3y ago
  6.7948900e-9 which has only the first 5 digits correct (if you assume it would have to have been 6.79490 to get 6 digits correct since that is what the correct value rounds to)
JKCalhoun 3y ago
The author's Taylor Series looked salvageable to me. The range between -Pi/2 and Pi/2 looked fine. That useable range can be re-used to serve all other portions of the circle. Add a conditional that applies (1-result) for some angles and you're golden.
Technotroll 3y ago
Looks to me like you could compute the top half, and then just repeat the rest as a kind of mirror function, that repeats with some set translation. Am I wrong here?
- dh2022 3y ago
  I thought the same thing as well. When I saw the Taylor sum diverging at the "edges" I thought - well Taylor sum converges absolutely only in a compact interval. An easy solution to this limitation is to use the periodicity of the cos/sin function...
deleted 3y ago
[deleted]
amiga386 3y ago
One of the things I loved most about reading kernel and libc (or rather, libm) sources was the floating-point and FP emulation code. I thought it was immensely cool to see what others just expected an FPU instruction to do, someone had not only written out in C, but also wrote comments explaining the mathematics (rendered in ASCII), often with numerical analysis about error propagation, accuracy, etc. Some examples: https://sourceware.org/git/?p=newlib-cygwin.git;a=blob;f=newlib/libm/math/e_pow.c https://sourceware.org/git/?p=newlib-cygwin.git;a=blob;f=new... https://git.kernel.org/pub/scm/linux/kernel/git/torvalds/linux.git/tree/arch/x86/math-emu/poly_sin.c https://git.kernel.org/pub/scm/linux/kernel/git/torvalds/lin... https://git.kernel.org/pub/scm/linux/kernel/git/torvalds/linux.git/tree/arch/x86/math-emu/README https://git.kernel.org/pub/scm/linux/kernel/git/torvalds/lin...
throwawaaarrgh 3y ago
Strange. This article was easy to read, simple formating, not driven by trends or weird internet culture. The weren't any bombastic claims, aggrandizing statements or dramatic opinions. Just interesting information presented clearly without being verbose. Almost like it was written by an adult. How did this end up on HN?
- smlacy 3y ago
  Yeah, very much unlike your comment. Strange, that.
mistercow 3y ago
This is a fine introduction to how you even approach the problem of computing transcendental functions. It would have been nice to have some discussion of accuracy requirements rather than just eyeballing it and saying “more accuracy than I need”. This is a spot where inexperienced devs can easily get tangled up. An attitude like “ten digits? That’s so many! I’m only making a game, after all” is exactly the sort of thing that gets you into trouble if you start accumulating errors over time, and this is particularly easy to do with trig functions.
eschneider 3y ago
Trig functions are where accuracy and performance go to die. Accumulated error is a thing, so when optimizing always consider exactly how you're going to be using those functions so you make the 'right' tradeoffs for your application. One size definitely doesn't fit all here, so test and experiment.
rjmunro 3y ago
Can you get some better accuracy by noting that cos(x) === sin(π/2 - x) and using e.g. the taylor expansion for sin when π/4 < x < 3π/4?
sojuz151 3y ago
I would say the benchmarks with lookup tables are bad. In the benchmark cpu will keep the table in cache but in a real program this cache would have to be shared with rest of the code/data. This would either kill the performance of cosine or rest of the app
mvcalder 3y ago
If you like this sort of thing there’s a great book: Methods and Programs for Mathematical Functions by Stephen Moshier. It covers the computation of all sorts of special functions. The code is available in the cephes library but the book may be out of print.
- jonsen 3y ago
  Not even a single hit on bookfinder.com!
midjji 3y ago
A good idea is to not to compute the values of cos from 0-2pi, but further reduce the range, using cos(a) = cos(-a), and cos(2a) = 1-2cos(a), or cos(a+pi/4) =... So we really only ever need to be able to compute cos in the range 0-pi/4. Then for further accuracy we can do the taylor expansion around pi/8. (or other approximations) finally the number of terms for a fixed accuracy varies with the distance from pi/8,
- londons_explore 3y ago
  I think all real implementations use a lookup table... Small lookup tables (eg. 8 elements) all work out far smaller in silicon area than even one multiplier...
  - midjji 3y ago
    The math.h cos from the article counterexample is pretty trivial, so can you give an example? I doubt linear interpolation between the 8 values is enough, but cubic probably would be. This would also be close to e.g. 8 third order taylor approximations of the function, and a comparison would be interesting.
pacaro 3y ago
I went through the same exercise implementing trig functions for scheme in webassembly... It was a rabbit hole for sure https://github.com/PollRobots/scheme/blob/main/scheme.wasm/src/library/trig.wat#L369 https://github.com/PollRobots/scheme/blob/main/scheme.wasm/s... [Edit: above is just the range from 0-π/4, the branch cuts for cos are at https://github.com/PollRobots/scheme/blob/main/scheme.wasm/src/library/trig.wat#L266 https://github.com/PollRobots/scheme/blob/main/scheme.wasm/s... ]
dang 3y ago
Url changed from https://web.archive.org/web/20210513043002/http://web.eecs.utk.edu/~azh/blog/cosine.html https://web.archive.org/web/20210513043002/http://web.eecs.u..., which points to this.
dang 3y ago
Related: Implementing Cosine in C from Scratch (2020) - https://news.ycombinator.com/item?id=30844872 https://news.ycombinator.com/item?id=30844872 - March 2022 (134 comments) Implementing cosine in C from scratch - https://news.ycombinator.com/item?id=23893325 https://news.ycombinator.com/item?id=23893325 - July 2020 (20 comments)
xorvoid 3y ago
How good do you need it? Lol. This is the approximation that I used in for the animated sinwave example for SectorC: y ~= 100 + (x*(157 - x)) >> 7 https://github.com/xorvoid/sectorc/blob/main/examples/sinwave.c https://github.com/xorvoid/sectorc/blob/main/examples/sinwav...
- JKCalhoun 3y ago
  Tried entering it in Grapher (on the Mac) and got what looked more or less like a parabola. My algebra is rusty — but do you kind of have an x-squared in there?
londons_explore 3y ago
Note that you can combine the lookup table with the taylor series expansion... And you can even use the same lookup table for each. That means with 2 table lookups in a 32 entry table, a single multiply and add (and a few bit shifts), you can get ~9 bits of precision in your result, which is fine for most uses. It also probably makes a sin operation take ~1 clock cycle on most superscalar architectures, as long as your workload has sufficient instruction parallelism. Note that a smaller table typically works out faster because 32 entries fit in the cache, whereas repeated random entry into a 1024 entry table would probably kick a bunch of other stuff out of the cache that you wanted.
londons_explore 3y ago
OP's lookup table benchmarking is bad. The "modd(x, CONST_2PI);" at the top will dominate, by far, the runtime. Anyone who wants performance measures angles using n bit fixed point math, mapping 0 to 2*pi as 0 to (2^n)-1.
omgmajk 3y ago
>Maybe don't stare at that for too long... I sure did stare at that for way too long.
amadio 3y ago
Taylor expansion is usually not so good for this, you'll fare better with either Legendre or Chebyshev polynomials. Robin Green has some excellent material on the subject, which I am linking below. Faster Math Functions: https://basesandframes.files.wordpress.com/2016/05/fast-math-functions_p1.pdf https://basesandframes.files.wordpress.com/2016/05/fast-math... https://basesandframes.files.wordpress.com/2016/05/fast-math-functions_p2.pdf https://basesandframes.files.wordpress.com/2016/05/fast-math... Even faster math functions GDC 2020: https://gdcvault.com/play/1026734/Math-in-Game-Development-Summit https://gdcvault.com/play/1026734/Math-in-Game-Development-S...
aidenn0 3y ago
Marz's taylor series implementation as-is is significantly faster (~60% the runtime) than the glibc implementation at 6 terms on My Machine(tm), rather than about the same as in TFA. I haven't compared to LERP lookup table yet though.
charlieyu1 3y ago
Related: https://news.ycombinator.com/item?id=35381968 https://news.ycombinator.com/item?id=35381968 Cosine Implementation in C Much better approximation with only 7 terms
lntue 3y ago
So in the implementation of `cos_table__LERP`, you did technically 2 step range reduction: 1. Reduce x = x mod 2pi 2. Reduce index = 10^n * (x / 10^-n), and i - index = 10^n * (x mod 10^-n) With limited input range and required precision as in the tests, you can combine these 2 range reduction steps: 1. Choose the reduced range as power of 2 instead of power of 10 for cheaper modulus operation, let say `2^-N = 2^-7`. 2. Avoid the division in `modd(x, CONST_2PI)` by multiplying by `2^N / pi`. 3. Avoid the round trip `double -> int -> double` by using the `floor` function / instruction. Here is the updated version of `cos_table__LERP` which should have higher throughput and lower latency: ``` double cos_table_128_LERP(double x) { x = fabs(x); double prod = x TWO_TO_SEVEN_OVER_PI; } ```
lntue 3y ago
So in the implementation of cos_table_*_LERP, you did technically 2 step range reduction: 1. Reduce x = x mod 2*pi 2. Reduce index = floor(x / 10^-n), and i - index = 10^n * (x mod 10^-n) With limited input range and required precision as in the tests, you can combine these 2 range reduction steps: 1. Choose the reduced range as power of 2 instead of power of 10 for cheaper modulus operation, let say 2^-N = 2^-7. 2. Avoid the division in modd(x, CONST_2PI) by multiplying by 2^N / pi. 3. Avoid the round trip double -> int -> double by using the floor function / instruction. Here is the updated version of cos_table_*_LERP which should have higher throughput and lower latency: double cos_table_128_LERP(double x) { x = fabs(x); double prod = x * TWO_TO_SEVEN_OVER_PI; double id = floor(prod); double x = prod - id; /* after this step, 0 <= x < 1 */ int i = ((int)id) & 0xFF; /* i = id mod 2^8 */ return lerp(x, COS_TABLE[i], COS_TABLE[i + 1]); } You can also optimize lerp a bit more with the formula: lerp(w, v1, v2) = (1 - w) * v1 + w * v2 = w * (v2 - v1) + v1 We do employ this range reduction strategy in a more accurate way for trig functions in LLVM libc: - with FMA instructions: https://github.com/llvm/llvm-project/blob/main/libc/src/math/generic/range_reduction_fma.h - without FMA instructions: https://github.com/llvm/llvm-project/blob/main/libc/src/math/generic/range_reduction.h
fallingfrog 3y ago
Why not do a table, but also store a table of the 1st derivative (which, for cosine, you could use the same table again but shifted)? Then, you could do a 2nd order fit like a spline instead of a straight line between your table values. Betcha it would be crazy fast.
adeon 3y ago
I like how there's lots of replies showing different ways to do this, improve it, additional nuances and viewpoints etc. lots of smart people. If I ever want better solutions for some programming problem, I'll write a post about it and try to get it in frontpage of HN :)
fargle 3y ago
It's unusual, outside x86, for these functions to be hardware accelerated. So just about every libm has to do it in software. Nearly all libm implementations around are based off Sun's fdlibm from the 80's-90's, including a bunch of the variants cited below, *bsds, etc. They are basically updated slightly, but you can see their heritage in the structure of the code. The original is found on netlib these days: https://www.netlib.org/fdlibm/k_cos.c https://www.netlib.org/fdlibm/k_cos.c 14th order polynomial, but only uses ~7 terms. It's supposed to have error < 1 ULP. For fdlibm, it's pretty readable compared to some of the other fun ones. I seem to remember sqrt being a bit gnarly.
machina_ex_deus 3y ago
I would've done it using complex exponentiation. And the exponential function is very easy to estimate fast: instead of Taylor series, use: e^x=(1+x/n)^n pick n=2^N, and that's just a bit shift and N repeated multiplications. Probably much faster and accurate.
Const-me 3y ago
I once did that as well: https://github.com/Const-me/AvxMath/blob/master/AvxMath/AvxMathTrig.cpp https://github.com/Const-me/AvxMath/blob/master/AvxMath/AvxM... The method is different, and the OP hasn’t mentioned it — high-degree minimax polynomial approximation.

220 points by alraj 3y ago

No need for archive.org, my website moved a few years ago: https://austinhenley.com/blog/cosine.html https://austinhenley.com/blog/cosine.html

butterlesstoast 3y ago
Thanks for taking the time to make such an incredible article. Made my day reading it this morning.
dang 3y ago
Fixed. Thanks!
capitainenemo 3y ago
BTW, WRT projects using a lookup table... https://hg.hedgewars.org/hedgewars/file/tip/hedgewars/uSinTable.pas https://hg.hedgewars.org/hedgewars/file/tip/hedgewars/uSinTa... I think there's even a github mirror someone made since that's where you were looking. Maybe a little out of date, but this is old code.
- capitainenemo 3y ago
  oh... also, sine table isn't just a performance thing. Hedgewars does fixed point custom math for deterministic lockstep, so keeps the math simple and predictable on various compilers/hardware.
  - nitwit005 3y ago
    That's something I'd wondered about with games that require determinism. Seems fairly challenging to deal with for that sort of game.
    - capitainenemo 3y ago
      https://nicolas.brodu.net/en/programmation/streflop/index.html https://nicolas.brodu.net/en/programmation/streflop/index.ht... is apparently pretty reliable, but yeah, it's a big problem. Spring and 0AD both struggled with it too.
- deleted 3y ago
  [deleted]

Helenarttr 3y ago

That was so amazing Information. https://www.ballsportspro.com/how-does-a-pickleball-ladder-work/ https://www.ballsportspro.com/how-does-a-pickleball-ladder-w...

xipix 3y ago

Nice. I'd love to see how this changes when you have SIMD and multiple cosines to compute in parallel. Also when you have to compute sine and cosine simultaneously which is often the case, and then you may be more interested in polar error than cartesian error.

p0nce 3y ago
You can do this with: https://github.com/RJVB/sse_mathfun https://github.com/RJVB/sse_mathfun Those functions do not have the same safety as libc, but even for scalar sin/cos it can be faster.

sampo 3y ago

> cosine repeats every 2pi > We can do better since the cosine values are equivalent every multiple of pi, except that the sign flips. There is one more step to take: Each π/2 long segment has identical shape, they are just pointing up or down (like you already noticed), and left or right. So you can reduce you basic domain down to not just 0...π, but to 0...π/2.

amelius 3y ago
But that takes extra clock cycles, since now you have to fiddle with the sign.
- Someone 3y ago
  I think it’s likely that you get that more than back by having to do (at least) one iteration fewer. Those who know more about this: educate me.
  - adgjlsfhk1 3y ago
    Pretty much every implementation does reduction mod pi/2. Switching polynomials is basically free (and halving the domain significantly reduces the number of terms needed).
IIAOPSW 3y ago
And you can keep going with this logic, building cos from 0-x by recursively calculating cos from 0-x/2 then mirroring over some line if needed. Eventually you hit machine precision limit aka the base case cos(0) = 1.
ajross 3y ago
And still further: the function can be piecewise decomposed into separate polynomials, which is what the pervasive SunPro code in everyone's libc does. There's one "mostly parabolic" function to handle the curvy top, and another "mostly linear" one to handle the zero crossing region.
gmiller123456 3y ago

jpfr 3y ago

The Taylor expansion locally fits a polynom based on the n first derivatives. If you want to find the "best" nth-degree polynom to approximate the sine function, functional analysis gives the tools for solving that optimization in closed form. By selecting an appropriate norm (in function space) you can minimize either the maximum error or the error integral over some range (e.g. the 0-\pi range). Here's a video on the subject. You might want to watch earlier ones also for more context. https://www.youtube.com/watch?v=tMlKZZf2Kac&list=PLdkTDauaUnQpzuOCZyUUZc0lxf4-PXNR5&index=28 https://www.youtube.com/watch?v=tMlKZZf2Kac&list=PLdkTDauaUn... Full disclosure, this is my university lecture on optimization that was recorded during Covid.

stephencanon 3y ago
It’s worth noting that for fixed-precision math libraries, we usually want the best approximation in the (possibly weighted) L-inf sense, which doesn’t generally have a closed-form solution, so we use iterative methods to generate approximations instead.
- bigdict 3y ago
  What's a good textbook intro to these methods?
  - remcob 3y ago
    I wrote this intro for a friend and since then more people have found it useful: https://2π.com/22/approximation https://xn--2-umb.com/22/approximation If you want a proper textbook, I recommend L.N. Trefethen (2019). “Approximation Theory and Approximation Practice.” Which is available online here: https://epubs.siam.org/doi/book/10.1137/1.9781611975949 https://epubs.siam.org/doi/book/10.1137/1.9781611975949
    - stephencanon 3y ago
      That's a really nice overview of Remez, thanks for sharing!

femto 3y ago

If doing signal processing, an optimisation is to use an N-bit integer to represent the range 0 to 2.pi. Some example points in the mapping: 0->0, 2^(N-2)->pi/2, 2^(N-1)->pi, 2^N(wraps to 0)->2.pi(wraps to 0). If your lookup table has an M-bit index, the index to the lookup table is calculated with: index = (unsigned)theta>>(N-M), where theta in the N-bit integer representing the angle. The fractional part, which can be used for interpolation, is: theta & ((1<<(N-M))-1). If you choose M=N=word size (16 bits is often nice), the angle can be used directly as an index. With 16-bits accuracy is typically good enough without interpolation and the entire trig operation reduces to an array indexing operation. This minimises the overhead of converting an angle to a table index.

rjmunro 3y ago
Why just signal processing? Surely this is good for storing e.g. rotations of characters in games, longitude, in fact almost anywhere where you want an angle where a full rotation is the same as no rotation at all, and more precision as you approach 0° is not useful.
- vikingerik 3y ago
  Accumulated error can matter, which is where you'd want that precision near 0°. If something is making tiny turns each frame for thousands of frames, floating-point inaccuracy can add up. (A smarter mathematical approach might remedy this, or might not; think of some game component that is supposed to jitter randomly each frame, maybe some particle system.) It's not common, but it's within hypothetical feasibility that the precision would be useful.
  - amitport 3y ago
    He is describing deterministic uniform round to nearest quantization. The most common smarter mathematical approach would be randomized rounding and it will solve this issue. (Issue being that the deterministic rounding is a biased estimator of the sin function)
    - 3y ago

Aardwolf 3y ago

" In all my time of coding, there has only been one situation in which I used cosine: games, games, games. function move(obj) { obj.x += obj.speed * Math.sin(obj.rotation); obj.y += obj.speed * Math.cos(obj.rotation); } " Why not store the velocity as a 2D vector instead? Then you still have to use cos/sin to compute this vector, but at least you don't need it every frame, plus often you don't need to use cos/sin to compute this vector either since forces that act on the velocity themselves can have an x and y component you can directly add to it

messe 3y ago
Sometimes it's just easier to work and think in polar coordinates. If that's the case, then (speed, rotation) is just as valid a choice as (vx, vy). For smaller games, the cost of computing sin and cos every frame is so astronomically small, that it's not even concerning yourself over.
- Aardwolf 3y ago
  Hmm, I usually find polar coordinates harder to reason with due to the gimbal lock :) In 2D there's no gimbal lock of course, but there's still a number that behaves differently than the others, the angle wrapping around in range 0..2pi! And of course the article is about fast computation of sin/cos so the option of not computing it should be considered
  - thehappypm 3y ago
    Can you explain gimbal lock? I can’t wrap my head around it
    - Aardwolf 3y ago
      The main thing is how when representing the 3D looking orientation as multiple angles (pitch=up/down, yaw=compass direction), you lose control over the yaw axis when looking straight up.

aebtebeten 3y ago

Niklaus Wirth also has some cosine code (probably meant more for software floating point, or fp FPGA blocks); I don't know how they compare with these approximations but his seem to be within 1e-6 of python's math.cos ... https://people.inf.ethz.ch/wirth/ProjectOberon/Sources/Math.Mod.txt https://people.inf.ethz.ch/wirth/ProjectOberon/Sources/Math....

dm319 3y ago

The answer to Cos(1.57079632) (in radians) will give you a clue as to how your calculator does it. See here [0] [0] https://www.reddit.com/r/calculators/comments/126st95/cos157079632_is_a_bit_of_a_torture_test/ https://www.reddit.com/r/calculators/comments/126st95/cos157...

the_third_wave 3y ago
Quite well when you realise this comes down to cos(π/2) to which the thing answers correctly: cos(π/2) = 0. Using HiPER Calc on Android
- dm319 3y ago
  The number is deliberately close but not the same as pi/2. See here: https://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv017.cgi?read=123814 https://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv017.c... And you can compare to wolfram alphas result.
aidenn0 3y ago
6.7948900e-9 which has only the first 5 digits correct (if you assume it would have to have been 6.79490 to get 6 digits correct since that is what the correct value rounds to)

JKCalhoun 3y ago

The author's Taylor Series looked salvageable to me. The range between -Pi/2 and Pi/2 looked fine. That useable range can be re-used to serve all other portions of the circle. Add a conditional that applies (1-result) for some angles and you're golden.

Technotroll 3y ago

Looks to me like you could compute the top half, and then just repeat the rest as a kind of mirror function, that repeats with some set translation. Am I wrong here?

dh2022 3y ago
I thought the same thing as well. When I saw the Taylor sum diverging at the "edges" I thought - well Taylor sum converges absolutely only in a compact interval. An easy solution to this limitation is to use the periodicity of the cos/sin function...

deleted 3y ago

[deleted]

amiga386 3y ago

One of the things I loved most about reading kernel and libc (or rather, libm) sources was the floating-point and FP emulation code. I thought it was immensely cool to see what others just expected an FPU instruction to do, someone had not only written out in C, but also wrote comments explaining the mathematics (rendered in ASCII), often with numerical analysis about error propagation, accuracy, etc. Some examples: https://sourceware.org/git/?p=newlib-cygwin.git;a=blob;f=newlib/libm/math/e_pow.c https://sourceware.org/git/?p=newlib-cygwin.git;a=blob;f=new... https://git.kernel.org/pub/scm/linux/kernel/git/torvalds/linux.git/tree/arch/x86/math-emu/poly_sin.c https://git.kernel.org/pub/scm/linux/kernel/git/torvalds/lin... https://git.kernel.org/pub/scm/linux/kernel/git/torvalds/linux.git/tree/arch/x86/math-emu/README https://git.kernel.org/pub/scm/linux/kernel/git/torvalds/lin...

throwawaaarrgh 3y ago

Strange. This article was easy to read, simple formating, not driven by trends or weird internet culture. The weren't any bombastic claims, aggrandizing statements or dramatic opinions. Just interesting information presented clearly without being verbose. Almost like it was written by an adult. How did this end up on HN?

smlacy 3y ago
Yeah, very much unlike your comment. Strange, that.

mistercow 3y ago

This is a fine introduction to how you even approach the problem of computing transcendental functions. It would have been nice to have some discussion of accuracy requirements rather than just eyeballing it and saying “more accuracy than I need”. This is a spot where inexperienced devs can easily get tangled up. An attitude like “ten digits? That’s so many! I’m only making a game, after all” is exactly the sort of thing that gets you into trouble if you start accumulating errors over time, and this is particularly easy to do with trig functions.

eschneider 3y ago

Trig functions are where accuracy and performance go to die. Accumulated error is a thing, so when optimizing always consider exactly how you're going to be using those functions so you make the 'right' tradeoffs for your application. One size definitely doesn't fit all here, so test and experiment.

rjmunro 3y ago

Can you get some better accuracy by noting that cos(x) === sin(π/2 - x) and using e.g. the taylor expansion for sin when π/4 < x < 3π/4?

sojuz151 3y ago

I would say the benchmarks with lookup tables are bad. In the benchmark cpu will keep the table in cache but in a real program this cache would have to be shared with rest of the code/data. This would either kill the performance of cosine or rest of the app

mvcalder 3y ago

If you like this sort of thing there’s a great book: Methods and Programs for Mathematical Functions by Stephen Moshier. It covers the computation of all sorts of special functions. The code is available in the cephes library but the book may be out of print.

jonsen 3y ago
Not even a single hit on bookfinder.com!

midjji 3y ago

A good idea is to not to compute the values of cos from 0-2pi, but further reduce the range, using cos(a) = cos(-a), and cos(2a) = 1-2cos(a), or cos(a+pi/4) =... So we really only ever need to be able to compute cos in the range 0-pi/4. Then for further accuracy we can do the taylor expansion around pi/8. (or other approximations) finally the number of terms for a fixed accuracy varies with the distance from pi/8,

londons_explore 3y ago
I think all real implementations use a lookup table... Small lookup tables (eg. 8 elements) all work out far smaller in silicon area than even one multiplier...
- midjji 3y ago
  The math.h cos from the article counterexample is pretty trivial, so can you give an example? I doubt linear interpolation between the 8 values is enough, but cubic probably would be. This would also be close to e.g. 8 third order taylor approximations of the function, and a comparison would be interesting.

pacaro 3y ago

I went through the same exercise implementing trig functions for scheme in webassembly... It was a rabbit hole for sure https://github.com/PollRobots/scheme/blob/main/scheme.wasm/src/library/trig.wat#L369 https://github.com/PollRobots/scheme/blob/main/scheme.wasm/s... [Edit: above is just the range from 0-π/4, the branch cuts for cos are at https://github.com/PollRobots/scheme/blob/main/scheme.wasm/src/library/trig.wat#L266 https://github.com/PollRobots/scheme/blob/main/scheme.wasm/s... ]

dang 3y ago

Url changed from https://web.archive.org/web/20210513043002/http://web.eecs.utk.edu/~azh/blog/cosine.html https://web.archive.org/web/20210513043002/http://web.eecs.u..., which points to this.

dang 3y ago

Related: Implementing Cosine in C from Scratch (2020) - https://news.ycombinator.com/item?id=30844872 https://news.ycombinator.com/item?id=30844872 - March 2022 (134 comments) Implementing cosine in C from scratch - https://news.ycombinator.com/item?id=23893325 https://news.ycombinator.com/item?id=23893325 - July 2020 (20 comments)

xorvoid 3y ago

How good do you need it? Lol. This is the approximation that I used in for the animated sinwave example for SectorC: y ~= 100 + (x*(157 - x)) >> 7 https://github.com/xorvoid/sectorc/blob/main/examples/sinwave.c https://github.com/xorvoid/sectorc/blob/main/examples/sinwav...

JKCalhoun 3y ago
Tried entering it in Grapher (on the Mac) and got what looked more or less like a parabola. My algebra is rusty — but do you kind of have an x-squared in there?

londons_explore 3y ago

Note that you can combine the lookup table with the taylor series expansion... And you can even use the same lookup table for each. That means with 2 table lookups in a 32 entry table, a single multiply and add (and a few bit shifts), you can get ~9 bits of precision in your result, which is fine for most uses. It also probably makes a sin operation take ~1 clock cycle on most superscalar architectures, as long as your workload has sufficient instruction parallelism. Note that a smaller table typically works out faster because 32 entries fit in the cache, whereas repeated random entry into a 1024 entry table would probably kick a bunch of other stuff out of the cache that you wanted.

londons_explore 3y ago

OP's lookup table benchmarking is bad. The "modd(x, CONST_2PI);" at the top will dominate, by far, the runtime. Anyone who wants performance measures angles using n bit fixed point math, mapping 0 to 2*pi as 0 to (2^n)-1.

omgmajk 3y ago

>Maybe don't stare at that for too long... I sure did stare at that for way too long.

amadio 3y ago

Taylor expansion is usually not so good for this, you'll fare better with either Legendre or Chebyshev polynomials. Robin Green has some excellent material on the subject, which I am linking below. Faster Math Functions: https://basesandframes.files.wordpress.com/2016/05/fast-math-functions_p1.pdf https://basesandframes.files.wordpress.com/2016/05/fast-math... https://basesandframes.files.wordpress.com/2016/05/fast-math-functions_p2.pdf https://basesandframes.files.wordpress.com/2016/05/fast-math... Even faster math functions GDC 2020: https://gdcvault.com/play/1026734/Math-in-Game-Development-Summit https://gdcvault.com/play/1026734/Math-in-Game-Development-S...

aidenn0 3y ago

Marz's taylor series implementation as-is is significantly faster (~60% the runtime) than the glibc implementation at 6 terms on My Machine(tm), rather than about the same as in TFA. I haven't compared to LERP lookup table yet though.

charlieyu1 3y ago

Related: https://news.ycombinator.com/item?id=35381968 https://news.ycombinator.com/item?id=35381968 Cosine Implementation in C Much better approximation with only 7 terms

lntue 3y ago

So in the implementation of `cos_table__LERP`, you did technically 2 step range reduction: 1. Reduce x = x mod 2pi 2. Reduce index = 10^n * (x / 10^-n), and i - index = 10^n * (x mod 10^-n) With limited input range and required precision as in the tests, you can combine these 2 range reduction steps: 1. Choose the reduced range as power of 2 instead of power of 10 for cheaper modulus operation, let say `2^-N = 2^-7`. 2. Avoid the division in `modd(x, CONST_2PI)` by multiplying by `2^N / pi`. 3. Avoid the round trip `double -> int -> double` by using the `floor` function / instruction. Here is the updated version of `cos_table__LERP` which should have higher throughput and lower latency: ``` double cos_table_128_LERP(double x) { x = fabs(x); double prod = x TWO_TO_SEVEN_OVER_PI; } ```

lntue 3y ago

So in the implementation of cos_table_*_LERP, you did technically 2 step range reduction: 1. Reduce x = x mod 2*pi 2. Reduce index = floor(x / 10^-n), and i - index = 10^n * (x mod 10^-n) With limited input range and required precision as in the tests, you can combine these 2 range reduction steps: 1. Choose the reduced range as power of 2 instead of power of 10 for cheaper modulus operation, let say 2^-N = 2^-7. 2. Avoid the division in modd(x, CONST_2PI) by multiplying by 2^N / pi. 3. Avoid the round trip double -> int -> double by using the floor function / instruction. Here is the updated version of cos_table_*_LERP which should have higher throughput and lower latency: double cos_table_128_LERP(double x) { x = fabs(x); double prod = x * TWO_TO_SEVEN_OVER_PI; double id = floor(prod); double x = prod - id; /* after this step, 0 <= x < 1 */ int i = ((int)id) & 0xFF; /* i = id mod 2^8 */ return lerp(x, COS_TABLE[i], COS_TABLE[i + 1]); } You can also optimize lerp a bit more with the formula: lerp(w, v1, v2) = (1 - w) * v1 + w * v2 = w * (v2 - v1) + v1 We do employ this range reduction strategy in a more accurate way for trig functions in LLVM libc: - with FMA instructions: https://github.com/llvm/llvm-project/blob/main/libc/src/math/generic/range_reduction_fma.h - without FMA instructions: https://github.com/llvm/llvm-project/blob/main/libc/src/math/generic/range_reduction.h

fallingfrog 3y ago

Why not do a table, but also store a table of the 1st derivative (which, for cosine, you could use the same table again but shifted)? Then, you could do a 2nd order fit like a spline instead of a straight line between your table values. Betcha it would be crazy fast.

adeon 3y ago

I like how there's lots of replies showing different ways to do this, improve it, additional nuances and viewpoints etc. lots of smart people. If I ever want better solutions for some programming problem, I'll write a post about it and try to get it in frontpage of HN :)

fargle 3y ago

It's unusual, outside x86, for these functions to be hardware accelerated. So just about every libm has to do it in software. Nearly all libm implementations around are based off Sun's fdlibm from the 80's-90's, including a bunch of the variants cited below, *bsds, etc. They are basically updated slightly, but you can see their heritage in the structure of the code. The original is found on netlib these days: https://www.netlib.org/fdlibm/k_cos.c https://www.netlib.org/fdlibm/k_cos.c 14th order polynomial, but only uses ~7 terms. It's supposed to have error < 1 ULP. For fdlibm, it's pretty readable compared to some of the other fun ones. I seem to remember sqrt being a bit gnarly.

machina_ex_deus 3y ago

I would've done it using complex exponentiation. And the exponential function is very easy to estimate fast: instead of Taylor series, use: e^x=(1+x/n)^n pick n=2^N, and that's just a bit shift and N repeated multiplications. Probably much faster and accurate.

Const-me 3y ago

I once did that as well: https://github.com/Const-me/AvxMath/blob/master/AvxMath/AvxMathTrig.cpp https://github.com/Const-me/AvxMath/blob/master/AvxMath/AvxM... The method is different, and the OP hasn’t mentioned it — high-degree minimax polynomial approximation.