8 ms·
How do computers calculate sine?
- vardump 3y agoCORDIC is how it's usually done in hardware (and FPGAs). https://en.wikipedia.org/wiki/CORDIC https://en.wikipedia.org/wiki/CORDIC
- adgjlsfhk1 3y agono it's not. cordic has awful convergence of 1 bit per iteration. pretty much everyone uses power series.
- f1shy 3y agoThat is 64 iterations for a double, that is nothing!
- adgjlsfhk1 3y ago53, but that's still a lot more than the 5th degree polynomial that you need.
- rcxdude 3y agoyeah, but 52 adds can be a lot cheaper than a few multiplies, if you're making them out of shift registers and logic gates (or LUT). in a CPU or GPU, who cares, moving around the data is 100x more expensive than the ALU operation.
- Const-me 3y ago> in a CPU or GPU, who cares, moving around the data is 100x more expensive than the ALU operation Moving data is indeed expensive, but there’s another reason to not care. Modern CPUs take same time to add or multiply floats. For example, the computer I’m using, with AMD Zen3 CPU cores, takes 3 cycles to add or multiply numbers, which applies to both 32- and 64-bit flavors of floats. See addps, mulps, addpd, mulpd SSE instructions in that table: https://www.uops.info/table.html https://www.uops.info/table.html
- adgjlsfhk1 3y ago> moving around the data is 100x more expensive than the ALU operation. This is exactly the problem with CORDIC. 52 dependent adds requires moving data from a register to the ALU and back 52 times.
- rcxdude 3y agoIt's the problem with CORDIC in that context, yes!
- ajross 3y agoActually pretty much everyone implements double precision sin/cos using the same (IIRC) pair of 6th order polynomials. The same SunPro code exists unchnaged in essentially every C library everywehre. It's just a fitted curve, no fancy series definition beyond what appears in the output coefficients. One for the "mostly linear" segment where the line crosses the origin and another for the "mostly parabolic" peak of the curve.
- perihelions 3y agoCORDIC is pretty obsolete, AFAIK. Its advantage is that its hardware requirements are absolutely tiny: two (?) accumulator registers, and hardware adders and shift-ers—I think that's all. No multiplication needed, in particular. Very convenient if you're building things from discrete transistors, like the some of those earlier scientific calculators! (Also has a nice property, apparently, that CORDIC-like routines exist for a bunch of special functions and they're very similar to each other. Does anyone have a good resource for learning the details of those algorithms? They sound elegant).
- pclmulqdq 3y agoCORDIC still is used in tiny microcontrollers (smaller than Cortex-M0) and in FPGAs when you are very resource-constrained. Restricting the domain and using Chebyshev/Remez is the way to go pretty much everywhere.
- f1shy 3y agoMultiplication is pretty much needed in cordic! And is far from obsolete! It works perfectly fine, and dont have any of the problems said in the article.
- perihelions 3y agoIt uses multiplication by powers of two, which is a floating-point bit shift.
- kevin_thibedeau 3y agoCORDIC doesn't use multipliers. That's the whole appeal for low performance hardware since it's all shifts and adds. It can still be useful on more capable platforms when you want sin and cos in one operation since there is no extra cost.
- dilyevsky 3y agoI think still used out of necessity when hw floating point not available (like fpgas)
- 3y ago
- bigbillheck 3y agoWhy would you ever use CORDIC if you had any other option?
- vardump 3y agoIt's great for hardware implementations, because it's simple and you get good/excellent accuracy. I wouldn't be surprised if that's still how modern x86-64 CPUs compute sin, cos, etc. That said, last time I had to do that in software, I used Taylor series. Might not have been an optimal solution. EDIT: AMD's Zen 4 takes 50-200 cycles (latency) to compute sine. I think that strongly suggests AMD uses CORDIC. https://www.agner.org/optimize/instruction_tables.pdf https://www.agner.org/optimize/instruction_tables.pdf page 130. Same for Intel, Tiger Lake (Intel gen 11) has 60-120 cycles of latency. Page 353. I'd guess usually ~50 cycles for Zen 4 (and ~60 for Intel) for float32, float64/float80 datatype. Denormals might also cost more cycles.
- bigbillheck 3y agoThey switched away from CORDIC at one point: https://www.intel.com/content/www/us/en/developer/articles/technical/the-difference-between-x87-instructions-and-mathematical-functions.html?wapkw=fsin https://www.intel.com/content/www/us/en/developer/articles/t... (there doesn't seem to actually be a linked article there, just the summary)
- vardump 3y agoPretty weird Intel's sine computation latency hasn't changed all that much over the years. Latencies have been pretty similar for 20 years. EDIT: That's a paper for a software library, not the CPU's internal implementation. Which is probably still done with CORDIC.
- bigbillheck 3y ago> EDIT: That's a paper for a software library, not the CPU's internal implementation. Unless you're seeing something I'm not, it's talking about x87, which hasn't been anything other than 'internal' since they stopped selling the 80486sx.
- eska 3y agoYou might also find this video interesting: "Finding the BEST sine function for Nintendo 64" https://www.youtube.com/watch?v=xFKFoGiGlXQ https://www.youtube.com/watch?v=xFKFoGiGlXQ
- t-3 3y agoLink doesn't appear to be valid, but aren't these usually precalculated and stored in a lookup table?
- lupire 3y agoThink about how big that lookup table would be for an 64bit double float. They use a lookup table for nπ/16, and then a polynomial to approximate the difference to other values.
- Etherlord87 3y agoThat would be quite a big lookup table... Half of float numbers are in -1…+1 range. You don't need the negative part, so that's a quarter of float numbers, over one billion numbers. And then some to get from 1 to π/2. And that's only float32!
- paulpauper 3y agosine is easy because the series is globally convergent and fast converging
- ot 3y agoDid you read the article? It is specifically about how the series looks simple, but the error is actually very bad if you do things naively.
- paulpauper 3y agoThat still makes it easier compared to computing constants in which the series are not globally convergent, like inverse trig functions. Obviously, you would have to break it apart to speed convergence.
- zgs 3y agoIt would also be extremely inaccurate. The x^n numerators grow very quickly and digits get lost because unlimited precision isn't available. Likewise, the n! denominators also grow rapidly. Then the series is alternating which means cancellation is happening for every added term. If you don't believe me try for x=10.
- deleted 3y ago[deleted]
- planede 3y agoYou need around 26 terms for x=10 if you do it without reduction and you want an accuratish result for double precision. You wouldn't evaluate the terms naively from left-to-right. x - x^3/3! + x^5/5! - ... = x * (1 - x^2/(2*3) * (1 - x^2/(4*5) * ... ) ) I just checked in python and you get a result that is around 1000*machine epsilon off. Not great, not terrible.
- duped 3y ago1 - cos^2(x), obviously
- nh23423fefe 3y agoinstead, you could just double negate to optimize away the square root -(-(sin(x))
- ChainOfFools 3y agoclear, but too verbose. 1/csc is what you want. or for style points just -cos'
- tails4e 3y agoI've seen the third order Taylor series used, but with the coefficients calculated at various offsets for a quarter wave. So you lookuo where you are in the quarter wave, then look up the 3 or 4 cofficients. This keeps the error somewhat bounded as the size of X is a small so the series does not diverge too much.
- amelius 3y agoAnd arcsin?
- paulpauper 3y agoharder due to convergence issues
- convolvatron 3y agohttp://steve.hollasch.net/cgindex/math/inccos.html http://steve.hollasch.net/cgindex/math/inccos.html is a great technique if you need a fast integer approximation for some some arbitrary sampling interval (i.e. motor control)
- phkahler 3y agoI've been doing FoC motor control for a long time and I've settled on a neat little fixed point approximation for sin/cos. I haven't been able to find the blog I got the idea from. It's accurate to 9 bits, but is very symmetric and hits 1,0,-1 exactly. It's also smooth which usually makes it better than lookup tables.
- warpech 3y agoThis made me realize that trigonometric functions are not deterministic across different CPU architectures, OS, and programming languages (floating point precision aside). E.g. I would assume that Math.sin(x) returns the same thing in NodeJS on Windows and Mac/M1, but it turns out it is necessarily so. https://stackoverflow.com/questions/74074312/standard-math-functions-reproducibility-on-different-cpus https://stackoverflow.com/questions/74074312/standard-math-f...
- adgjlsfhk1 3y agosome languages (e.g. Julia) provide their own math library do that you get the same results across across operating systems.
- TylerE 3y agoSafer to assume that floats are never deterministic.
- jacobolus 3y agoFloats follow a clear specification which determines precisely how basic arithmetic should work. They should work the same on all popular modern platforms. (Whether specific software libraries are the same is a separate question.)
- microtherion 3y agoP.J. Plauger's _The Standard C Library_ provides an implementation for all functions in the (then) C standard: https://www.amazon.com/Standard-Library-P-J-Plauger/dp/0138380120?ref_=ast_author_dp&dib=eyJ2IjoiMSJ9.J2FZMVfhMnUjy3nORwaNfJ39GKhZvMa1t-YBXfQeaEgAuQf63AYkxWWCauQjjBeo9Z3_OrNF4PZDHa-l_tdvzR3ooYzJBGyfigUwXLxNHiszIfSYtPsgIzjEKpUmRBVOSMlrnXBoG26XMFM6WGfX6gSXoSJCch-ZygQNZo-OuqAfqvHlSY4NZBnBxp3ORQySq32fkwuWqv516zBqZMCmp6fPUdbJG1rJgG9Z5yJxIc8.Kz9qxgL-6ZhzUqEeAARs4sTCGFg8uFGlXxRP21iGrRc&dib_tag=AUTHOR https://www.amazon.com/Standard-Library-P-J-Plauger/dp/01383...
- dboreham 3y agoThis was my first use of open source, around 1978. I wondered how calculators and computers did trig functions, and was also using Unix V7. We had a large disk and kept the source on line. So I was able to find this: https://www.tuhs.org/cgi-bin/utree.pl?file=V7/usr/src/libm/sin.c https://www.tuhs.org/cgi-bin/utree.pl?file=V7/usr/src/libm/s... and from there this book: https://www.biblio.com/book/computer-approximations-john-f-hart-e/d/1504828236 https://www.biblio.com/book/computer-approximations-john-f-h...
- perihelions 3y agoIs this still current? The paper has a publication year of 1999.
- jxy 3y agoIt's much clearer if you read one of the source code of the libm. Plan 9: https://9p.io/sources/plan9/sys/src/libc/port/sin.c https://9p.io/sources/plan9/sys/src/libc/port/sin.c Freebsd: https://cgit.freebsd.org/src/tree/lib/msun/src/k_sin.c https://cgit.freebsd.org/src/tree/lib/msun/src/k_sin.c
- anthk 3y agoDoes 9front keep the same sin.c implementation?
- jxy 3y agohttp://git.9front.org/plan9front/plan9front/HEAD/sys/src/libc/port/sin.c/f.html http://git.9front.org/plan9front/plan9front/HEAD/sys/src/lib...
- lifthrasiir 3y agoFreeBSD code is missing the range reduction step (it's named a "kernel" for the reason): https://cgit.freebsd.org/src/tree/lib/msun/src/e_rem_pio2.c https://cgit.freebsd.org/src/tree/lib/msun/src/e_rem_pio2.c
- planede 3y agoAh, I always referred to the musl implmenetation, but I just now realized that they copied the Freebsd one.
- jxy 3y agoit says "Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved." and the directory name is literally: `msun`
- deepthaw 3y agoIgnorant question: Given the ridiculous number of transistors and so on we can use in CPUs and GPUs nowadays how feasible is a relatively huge trig lookup table burned into rom?
- zokier 3y agoThere is huge number of 64bit floats and huge portion of those are between 0..pi/2.
- bigbillheck 3y agoSeems like a terrible idea on latency grounds alone.
- aidenn0 3y agoA lookup table (for any function) that covered all values between 0 and 1 in single precision, would be ~4GB; there are approximately 1B values between 0 and 1, and the result of each value is 4 bytes. Such a table for double-precision would be much, much larger.
- azhenley 3y agoI blogged my adventure of implementing cosine from scratch and how others have done it: https://austinhenley.com/blog/cosine.html https://austinhenley.com/blog/cosine.html
- toolslive 3y agoAfter reducing the interval, you don't want to use the Taylor series as you're building an approximation that's really good in 0 but not so good moving away from 0. It's better to use an interpolating polynomial (Chebychev comes to mind) over the whole target interval.
- paulpauper 3y agoThere are many ways to do this. It's not a difficult problem unless memory is constrained.
- demondemidi 3y agoI thought modern CPUs since the late 1980's used a lookup table for trig/transcendental functions. Is the LUT just an expansion of the polynomial? I never really understood how FPUs worked...
- Someone 3y agoThat would take way too much room. A full lookup table would have 2^64 entries of 64 bits each, at 2^70 bits of ROM. For comparison: - Apple’s M2 Ultra has about 134 billion transistors. That’s about 2^38. - Avogadro’s number is about 2^79. Reducing the argument to a small range around zero decreases that a lot, but not enough by a far stretch. There are 2^52 doubles in [0.5, 1.0), 2^52 more in [0.25, 0.5], 2^52 more in [0.125, 0.25], etc. so you’d still easily need 2^52 entries or 2^58 bits (likely way, way more)
- demondemidi 3y agoThey DO use a lookup table because that’s what the FDIV but came from: https://en.m.wikipedia.org/wiki/Pentium_FDIV_bug https://en.m.wikipedia.org/wiki/Pentium_FDIV_bug “It is implemented using a programmable logic array with 2,048 cells, of which 1,066 cells should have been populated with one of five values: −2, −1, 0, +1, +2.” Not sure what you’re trying to demonstrate, they wouldn’t store every single float!! I hope don’t program. ;)
- lifthrasiir 3y agoOlder CPUs generally have used CORDIC (which does use LUT but that's only a part of the algorithm) due to its simplicity and compactness, while later CPUs with extensive microcode support would do the same thing as software implementations.
- aidenn0 3y agoTFA says they use a 32 entry LUT, then do some math on the result.
- eh_why_not 3y agoAnyone experienced with the Remez algorithm mentioned at the end of the article? The degree-9 polynomial, said to be a thousand times better than the original Taylor approximation in maximum error, also appears to be very close to the Taylor series in the first place. Rounding the Taylor coefficients to 6 digits after the decimal: 1/3! = 0.166667 1/5! = 0.008333 1/7! = 0.000198 1/9! = 0.000027(56) The first 2 are exact, the third is 5 digits only (so 0.000190), and the fourth is more different starting from the 6th digit (0.000026019). The delta in the 9-th order is expected if you were to truncate the Taylor series starting from the 11th order to infinity (+ x^11 / 11! - x^13/13! ...).
- stephencanon 3y agohttps://en.wikipedia.org/wiki/Remez_algorithm https://en.wikipedia.org/wiki/Remez_algorithm It’s a very simple iterative algorithm, essentially the dumbest thing that could possibly work (like most good algorithms). It fails to converge for functions that have poles nearby unless you have a very good initial guess (the Chebyshev or Carathéodory-Fejér approximants are ~always good starting points and easily computed). In practice you want to optimize a weighted L-inf norm rather than absolute, because floating-point errors are measured in a relative norm.
- olooney 3y agoI don't think the polynomial given in the article was calculated via Remez. Perhaps the author merely meant it as an illustration. Here is Wolfram Alpha plotting the error of equation from the article: https://www.wolframalpha.com/input?i=plot+sin%28x%29+-+P%28x%29+on+%5B0%2C+pi%2F2%5D+where+P%28x%29+%3D+x+-0.166667x%5E3+%2B+0.00833x%5E5+-0.00019x%5E7+%2B+2.6019%5Ccdot10%5E%7B-6%7Dx%5E9 https://www.wolframalpha.com/input?i=plot+sin%28x%29+-+P%28x.... You'll note that the error is small near zero, and large near pi/4, which is characteristic of a Taylor series around zero, and not the characteristic "level" oscillation characteristic of Remez approximations[1]. Note also that the polynomial only includes odd terms, which is not something I would expect from Remez unless it was run on the symmetric interval [-pi/4, pi/4]. I ran Remez for the problem and after 4 iterations obtained a degree 8 polynomial with error less than 1e-9, but it didn't look anything like the polynomial given in the article. 1.7209863008943345e-05x^8 -0.00024575124459624625x^7 +7.194649190849227e-05x^6 +0.008268794893899754x^5 +3.425379759410762e-05x^4 -0.16667692317020713x^3 +1.5422400957890642e-06x^2 +0.9999999106213322x +8.526477071446453e-10 Although of course the first few digits of the low-order matching terms will be very similar - any polynomial approximation method will agree there because they are fitting the same function, after all. But by the time we reach x^5 or x^7 the agreement is very loose, really only the first couple digits are the same. [1]: https://en.wikipedia.org/wiki/Approximation_theory https://en.wikipedia.org/wiki/Approximation_theory
- staplung 3y agoI recently learned how Doom was ported to the SNES. It's quite impressive. The SNES hardware was nowhere near fast enough to do all the trig calculations needed for the game but cartridge based games had a trick up their sleeve: they could include actual hardware inside the cart that the game code could make use of. It was more expensive but if you expected to sell a boatload of copies, it could be worth it. However, even using extra hardware wasn't enough in this case. So they pre-calculated lookup tables for sine, cosine, tangent etc. for every angle at the necessary precision. They were helped by the fact that the game resolution in this case was fairly low. If you're interested, you can peruse the C code that was used to generate the tables. Here's the file for sine/cosine: https://github.com/RandalLinden/DOOM-FX/blob/master/source/mksin.c https://github.com/RandalLinden/DOOM-FX/blob/master/source/m...
- pillusmany 3y agoGames targetting pre-Pentium PCs also used precomputed trig tables. Pentium was fast enough that it didn't matter as much. Just a few years later it was slower to read a trig precomputed table.
- xarope 3y agoin other words, for those of us who remember, they used the equivalent of a slide rule
- fodkodrasz 3y agoMore like a Trigonometry table, which predates even slide rules: https://en.wikipedia.org/wiki/Mathematical_table https://en.wikipedia.org/wiki/Mathematical_table
- BD103 3y agoYup, I remember watching a video about how the RAM bus is the bottleneck when running Super Mario 64 on the N64. The original implementation used trig lookup tables, but the person optimized it by instead using Taylor series (I think) and some negation / shifting.
- Solvency 3y agoWhy isn't this just done with an industry standard lookup table these days?
- simonblack 3y agoThe same way you can eat an elephant: one byte at a time. Any calculating job can be undertaken by a proper Turing machine. You just have to keep in mind the old triangle of usage: Cost, capability and speed. If a human can calculate a sine, so can any full-blown computer.
- mettamage 3y agoWhat’s the best way to calculate it by hand? I’m brushing up my math basics (I graduated CS while dodging the math requirements) and it frustrates me that in trig I need to remember values at all. The values such as sqrt(2)/2 make sense but how hard is it to calculate sin(5 degrees) by hand?
- HenryPrickett 3y agoUse a Taylor series with a four function calculator. 0 is a decent approximation of sin near 0. x is a better one. x - x^3/6 is an even better one. x - x^3/6 + x^5/120 ... Note that x here is in radians rather than degrees so convert (degrees * pi/180) first. Repeat until you're satisfied with how many stable digits you get
- empath-nirvana 3y agohttps://www.youtube.com/watch?v=3d6DsjIBzJ4 https://www.youtube.com/watch?v=3d6DsjIBzJ4
- _v7gu 3y agoFive is small enough that you can get away with sin(5 deg) = 5pi/180
- sema4hacker 3y agoI remember seeing the source code for a version of SpaceWar! running on an Adage Graphics Terminal (a one's complement machine) around 1970 that used a precomputed sine table. I wonder what the first program was that ever used a precomputed trig table.
- Tade0 3y agoI was truly amazed when my high school computer science teacher expanded the sine function into a Taylor series for us to implement in class. Couldn't wrap my head around the concept until the topic was revisited in college, but idea was there and helped me understand the tradeoffs brought by "fast math" libraries I was using during high school.
- aidenn0 3y agoJust using the 32 entry LUT and the small-angle approximations (sin(x) = x, cos(x) = 1-x^2/2) lets you calculate sin(x) within +/- 0.00015, which is rather impressive for something that can be done quickly by hand. If you use cos(x) = 1 then you are still accurate to within 1% [edit] I think I also found an error in TFA? It seems like picking the "best" N would allow r to be in the range -pi/32 < r < pi/32, which makes the 3rd order Taylor series have an error of 4e-12, significantly better than the error range for 0 < r < pi/16