9 ms·
Implementing Cosine in C from Scratch (2020)
- Gigachad 5y agoCan someone explain to me why calculating cosine seems to be one of the most complicated tasks? I have tried to read about it a few times and never really found an explanation for calculating it which makes sense.
- msk-lywenn 5y agoMy favorite implementation of a super fast but not that accurate cosinus is from a now defunct forum. Good thing we have webarchive ! https://web.archive.org/web/20111104163237/http://www.devmaster.net/forums/showthread.php?t=5784 https://web.archive.org/web/20111104163237/http://www.devmas...
- karma_fountain 5y agoYou can also use the Chebyshev polynomials to economise a higher degree taylor series with bit of a loss of accuracy.
- hatsunearu 5y agoSince cosine is even and periodic, you really only need like, up to x = 0 to pi/4 and the rest you can infer just from that data.
- mhh__ 5y agoThat is how it's often implemented i.e. shift back to the domain you actually need to be able to compute.
- dekhn 5y agolookup tables and lerp are the most common solution I've seen (for a wide range of functions). You can also memorize with an LRU cache.
- phkahler 5y ago>> lookup tables and lerp are the most common solution I've seen 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.
- LeifCarrotson 5y ago> Then you rotate that by the fine angle 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.
- marginalia_nu 5y agoRemember that cos(a+b)=cos(a)cos(b)-sin(a)sin(b)
- LeifCarrotson 5y agoAh, now it makes sense. Thanks for the reminder of 8th grade trig...it's been a while!
- phkahler 5y ago
- skulk 5y agoI learned how to write programs from Scratch[0]. Back in the pre-1.4 days (IIRC) there was no native sin/cos functions, so they had to be "implemented" somehow to translate arbitrary angles into delta (x,y). The way that I found people did this was as follows (for sin(A), cos(A)): Invisible sprite: [ when I receive "sine" ]: [ go to center ] [ turn A degrees ] [ move 1 step forward ] Now, anyone that wants to know sin(A), cos(A) just has to read that sprite's X and Y position. [0]: https://scratch.mit.edu https://scratch.mit.edu
- em3rgent0rdr 5y agoI guess scratch (or LOGO) is Turing complete, so someone could make a CPU built using a little turtle in a box where all functionality is built out of LOGO commands.
- Asraelite 5y agoIt's not just Turing complete but even practical. People have built entire CPU emulators that run pretty fast, e.g. https://scratch.mit.edu/projects/175052082/ https://scratch.mit.edu/projects/175052082/
- em3rgent0rdr 5y agoAmazing!
- deleted 5y ago[deleted]
- ThePhysicist 5y agoI remember doing this in Macromedia Flash as I wanted to make stuff go in circles and the language didn't have any higher math functions back then. Good ole times.
- dancsi 5y agoThe implementation of the __cos kernel in Musl is actually quite elegant. After reducing the input to the range [-pi/4, pi/4], it just applies the best degree-14 polynomial for approximating the cosine on this interval. It turns out that this suffices for having an error that is less than the machine precision. The coefficients of this polynomial can be computed with the Remez algorithm, but even truncating the Chebyshev expansion is going to yield much better results than any of the methods proposed by the author.
- l- 5y agohttps://github.com/ifduyue/musl/blob/master/src/math/__cos.c https://github.com/ifduyue/musl/blob/master/src/math/__cos.c
- llimllib 5y ago> Input y is the tail of x. What does "tail" mean in this context?
- lifthrasiir 5y agoIt is a double-double representation [1], where a logical fp number is represented with a sum of two machine fp numbers x (head, larger) and y (tail, smaller). This effectively doubles the mantissa. In the context of musl this representation is produced from the range reduction process [2]. [1] https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic https://en.wikipedia.org/wiki/Quadruple-precision_floating-p... [2] https://github.com/ifduyue/musl/blob/6d8a515796270eb6cec8a278cb353a078a10f09a/src/math/cos.c#L68-L76 https://github.com/ifduyue/musl/blob/6d8a515796270eb6cec8a27...
- jhgb 5y agoDoes it make sense to use a double-double input when you only have double output? Sine is Lispchitz-limited by 1 so I don't see how this makes a meaningful difference.
- aaaronic 5y agoI teach lookup tables specifically for sine/cosine (we just use one with a phase shift and modulus) on the Gameboy Advance. The hardware has no floating point unit support, so we do it in 8.8 fixed point with 1-degree resolution. 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.
- aaaronic 5y agoMandatory Physics "troll" comment: For small angles, sin(x) ~= x and cos(x) ~= 1 -- the ["small angle approximation"](https://en.wikipedia.org/wiki/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.
- jhgb 5y agoThis a perfectly nice Programming 101 example of a recursive function -- if an angle is too large, calculate its sine using the sine of the half-angle, otherwise return the angle (or, if you're fancy, some simple polynomial approximation of the sine). I'm sure everyone did this in school (we did, in fact).
- gugagore 5y agoIt's not obvious when you should switch over to an approximation (base case), so I'd say it's not a good example to introduce recursion. I have never seen it, and I did not do it in school.
- aaaronic 5y agoWhen the difference between iterations is less than whatever threshold of accuracy you're looking for, I'd assume. "Relaxation" algorithms for calculating fields work that way.
- jhgb 5y agoIt's not obvious but that makes it a nice example in parameter optimization as well.
- scythe 5y agoIt's "obvious" if you remember how the Taylor series for sin(x) works. The error between x and sin(x) is bounded above by x^3/6. So there you go.
- SnowHill9902 5y ago
- xioxox 5y agoHere's a solution which is vectorizable: http://gallium.inria.fr/blog/fast-vectorizable-math-approx/ http://gallium.inria.fr/blog/fast-vectorizable-math-approx/ I've also read that Chebyshev polynomials are far better for approximating functions than Taylor series (eg https://en.wikipedia.org/wiki/Approximation_theory https://en.wikipedia.org/wiki/Approximation_theory)
- freemint 5y agoThis matters less once you put things through the horner scheme. And if you are not you are leaving performance on the table.
- arthur2e5 5y agoHorner does not conflict with vectorization in this case. If you read the code, it actually uses Horner — the vectorizability applies to when you call the online function in a loop.
- waynecochran 5y agoYou are doing the polynomial evaluation wrong for Taylor Series. First of all you should be using Horner's Method to reduce the number of multiplications needed and then you should be using FMA to increase precision of the solution. https://momentsingraphics.de/FMA.html https://momentsingraphics.de/FMA.html
- naoqj 5y ago>I couldn't find any notable projects on GitHub that use a table for trig functions, but I'm sure they exist. Doom?
- freemint 5y agoJulia?
- adgjlsfhk1 5y agoJulia uses a very similar approach to musl here. Reduction mod 2pi followed by polynomial kernel. We use tables for `exp` and `log`, but the trig functions have good reductions so we use those instead. The implementation is here https://github.com/JuliaLang/julia/blob/03af78108afe6058f54cf766d8981d26a11c51f7/base/special/trig.jl#L29 https://github.com/JuliaLang/julia/blob/03af78108afe6058f54c... for anyone interested.
- snicker7 5y agoIn the Julia repl, just type '@less sin(pi)' to view the source code of sin. I don't know of any other programming language that makes it this easy to look at its internals.
- adgjlsfhk1 5y agoThat's how we nerd snipe people into contributing to Base (only 10% joking).
- zekica 5y agoDoom also uses the famous fast inverse square root method which although I know how it works still amazes me.
- naoqj 5y agoIt doesn't, to the best of my knowledge. It appeared in an id game for the first time in Quake 3.
- ok123456 5y agoA lot of Basic implementations, that had software support for floating point numbers, implemented it like this. Microsoft BASIC used a 6th degree polynomial iirc.
- nonrandomstring 5y agoAudio stuff uses different methods depending on the application. Additive oscillators need to be accurate and can use some closed forms and series, mostly people use Taylor and quadratic interpolation as the article shows. For tables remember you only really need to compute amd store one quadrant and shift it around. But for LFOs people use all manner of dirty and quick methods like Bhaskara. https://en.wikipedia.org/wiki/Bhaskara_I%27s_sine_approximation_formula https://en.wikipedia.org/wiki/Bhaskara_I%27s_sine_approximat...
- SeanLuke 5y agoA LERP lookup table is fast; but only somewhat slower, and probably quite a bit more accurate, would be a lookup table interpolated by a Catmull-Rom spline.
- dilyevsky 5y agoI remember researching cordic for some pet fpga project i had and finding some libs/discussions so it’s definitely used.
- blt 5y agoIf you want sinusoidal oscillation over time, then you can integrate the second-order differential equation x'' = -ω² x instead. This only costs a few arithmetic instructions per time step. On the other hand, it adds dependency on a state other than the current time, and you must be careful to use energy-conserving integration.
- duped 5y agoHere's a quick way of doing it with one matrix multiplication, you just need to precompute cos(ω) and sin(ω) with a bit of intuition: [ 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 There are two ways to look at this, from the systems perspective this is computing the impulse response of a critically stable filter with poles at e^(+/-jω). Geometrically, it's equivalent to starting at (1, 0) and rotating around the unit circle by ω radians each time step. You can offset to arbitrary phases. 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.
- em3rgent0rdr 5y agoSpeaking of lookup tables, I recall reading the DX-7 (https://www.righto.com/2021/11/reverse-engineering-yamaha-dx7.html https://www.righto.com/2021/11/reverse-engineering-yamaha-dx... shared a few months ago on HN), that DX-7 has a 4096-entry 14-bit sine ROM built into the chip itself.
- anonymousiam 5y agoI did the same thing for sqrt() seven years ago. It benchmarked at about five times faster than math.h/libm, but the trade was a reduction in the maximum value of the result. My version would not produce accurate results for inputs greater than 2^16. It did work very well for generating real-time laser pointing X/Y coordinates for an Archimedean Spiral. https://en.wikipedia.org/wiki/Archimedean_spiral https://en.wikipedia.org/wiki/Archimedean_spiral (This was on a MSP430 platform with a FPU that only did multiplication.)
- adgjlsfhk1 5y agoI'm surprised this was faster than an initial guess (half the exponent and mantisa) followed by newton's method.
- anonymousiam 5y agoNewton's method involves division, which is a problem on an embedded platform with a limited FPU that can only do multiplication. I tried several classic approaches and benchmarked all of them against the function in libm. I reviewed what I did and it turns out that my own memory is not that great. My version (nsqrt) is limited to input values of less than 2^14 and not 2^16 as I said above. Also, it's only 50% faster and not 5x as I said above. 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.
- deepspace 5y agoI was surprised that pre-calculating xx as double xx = xx made a difference. The author's compiler must not be great at optimization if this is the case.
- FartyMcFarter 5y agoI haven't checked, but off the top of my head this could be because this optimisation changes the precision of the result enough for it to be an illegal optimisation. 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.
- gweinberg 5y agoI think you should be able to get good results for precalculating the cosines for a set of angles and then using the formula for the cosine of the sum of 2 angles combined with the taylor series. The series converges very quickly for very small angles, right?
- marcodiego 5y agoI don't know if this benchmark is trustworthy. Using a lookup table pollutes cache. It definitely needs to be compared with interactions on other parts of the code to get a good performance measurement.
- stncls 5y agoThis is a very nice article, and I think the three initial objectives are fully attained (simple enough, accurate enough, fast enough). 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/#sm4fjqtvjyn36dyedncxh 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 Lookup table sizes: 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 bytes
- adgjlsfhk1 5y agoHow long does the cache thrashing loop take compared to 500 iterations? I'm not sure how much you can trash the cache while still being performance bottle-necked by a table-based implementation of a trig function.
- stncls 5y agoYes of course the cache thrashing takes a lot longer. But I modified the code to time only the cos() functions (see diff), using the rdtsc instruction. That's why it had to be every 500 iterations: If it was after each iteration, then rdtsc itself would become the bottleneck!
- adgjlsfhk1 5y agoMy point wasn't about how you were timing the code, what I meant is in the real world, I don't think there are programs where this will be an issue because if they are thrashing their caches too much, they won't be bottlenecked on trig. The advantage of table based functions is that they are really fast if you are calling a lot of them at once, but that's also the only case where they need to be fast. If only 1/10000 instructions is trig, it's OK if the trig is a little slower.
- anugrahit 5y agohttp://www.realitypixels.com/turk/computergraphics/FixedPointTrigonometry.pdf http://www.realitypixels.com/turk/computergraphics/FixedPoin... 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)
- aj7 5y agoSo 57 years ago, there weren’t any programming books. Well, there was one. McCracken and Dorn, Numerical Methods and Fortran Programming. A kindly childless couple, Ben and Bluma Goldin, who lived in the same apartment house in Brooklyn, bought me this book as a gift. An axiomatic lesson from the book was distrust of Taylor expansions, and the necessity of moving the expansion point of any function to a favorable location. Because harmonic functions are very well behaved, it should also be obvious that even a very coarse lookup table + interpolation should converge to high accuracy quickly. Finally, use of trigonometric identities to zero in on the desired angle is helpful. OK I’ll shut up now.
- chmaynard 5y agoWhy tell us the couple's names and tell us they were childless? Strange. I'm not inclined to trust anything you said after that.
- unethical_ban 5y agoIt's common in storytelling to add details. It's frustrating you call it bad faith because you don't get it.
- aj7 5y agoThey are long dead. It was sad. They doted on me and my sister because they had no kids. He was a salesman. He knew how to get books in Manhattan. He used to get my sister classics with the most beautiful color plate illustrations. In the sixties, these had little value. Today?
- eesmith 5y agoGoing on a tangent, as I understand it, the first programming book was "The Preparation of Programs for an Electronic Digital Computer: With special reference to the EDSAC and the Use of a Library of Subroutines", by Maurice Wilkes, David Wheller and Stanley Gill (1951). https://archive.org/details/programsforelect00wilk https://archive.org/details/programsforelect00wilk
- user_7832 5y agoSlightly off topic but kudos to the author for explaining the process extremely clearly. As someone without experience of low level programing this was quite understandable for me (albeit with uni-level maths knowledge).
- phonebucket 5y agoThis is a fun article! Alternative avenues that complement the approaches shown: - Padé Approximants (https://en.wikipedia.org/wiki/Pad%C3%A9_approximant 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_approximation_formula 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.
- nwallin 5y agoI have a few cases where I think a Padé approximation could really help me, but I've never been able to figure out, given a function, how to get the coefficients. Do you have any suggested reading?
- m1kal 5y agoI learned from Bender, Orszag "Advanced Mathematical Methods for Scientists and Engineers I".
- charlieyu1 5y agoI think off-zero Taylor series along with lookup table should be at least tried. You can approximate cos(x+a) using only cos(a) and sin(a) and the series converges really fast
- kevin_thibedeau 5y agoRemez always beats Taylor series.
- dheera 5y ago> So we are really only calculating cosine from 0 to pi. 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?
- dboreham 5y agoThis takes me back in time, to c. 1982 when I got curious about how computers calculated trig functions. Luckily I was using one of the first OSS systems (modulo AT&T lawyers): Unix V7. So I pulled up the source and took a look. You can still see it here: https://github.com/v7unix/v7unix/blob/master/v7/usr/src/libm/sin.c https://github.com/v7unix/v7unix/blob/master/v7/usr/src/libm... There's a comment referencing "Hart & Cheney", which is this book : https://www.google.com/books/edition/Computer_Approximations/jPVQAAAAMAAJ?hl=en https://www.google.com/books/edition/Computer_Approximations... Also luckily I had access to a library that had the book which I was able to read, at least in part. This taught me a few things : that the behavior of a computer can be a white box, not a black box; the value of having access to the source code; to use comments effectively; that for many difficult tasks someone probably had written a book on the subject. Also that polynomials can approximate transcendental functions.
- nwallin 5y agoIt's all well and good to have a library with fast sin/cos/tan functions, but always remember to cheat if you can. 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)> Instead, you might do: 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 you've replaced a sin/cos pair with 6 elementary float operations. 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.
- djmips 5y agoNice! Also see https://fgiesen.wordpress.com/2010/10/21/finish-your-derivations-please/ https://fgiesen.wordpress.com/2010/10/21/finish-your-derivat...
- todd8 5y agoThe article uses a game movement function as a possible motivation for needing the cosine. If object.rotation is constant as the spiral animation seems to suggest, there is an important optimization for rapidly computing the sequence: cos(theta), cos(theta + delta), cos(theta + 2delta), cos(theta + 3delta), ... 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)) The values for a and b must be calculated, but only once if delta is constant: a = 2 sin^2(delta/2) b = sin(delta) By using these relationships, it is only necessary to calculate cos(theta) and sin(theta) once for the first term in the series in order to calculate the entire series (until accumulated errors become a problem). [1] Press, William H. et. al., Numerical Recipes, Third Edition, Cambridge University Press, p. 219
- arthur2e5 5y agoIt’s an old article, so… guess I won’t beat a dead horse as much. * 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 https://stackoverflow.com/a/52841086 Others have already mentioned the Remez and the Horner stuff. Won’t repeat here.
- sylware 5y agoI have been looking for papers which are _very clean and serious_ mathematical proofs of accuracy for approximants of transcendental functions (remez/chebysheve/etc). Any pointers?
- userbinator 5y agoHowever, I can't find any references of whether it is used much today. CORDIC is still very common in embedded systems, especially with CPUs that have no hardware multiply/divide.
- helltone 5y agoWhat would it take to have a parameterized algorithm where you specify the precision you want and get a fast implementation for that precision?
- dietrichepp 5y agoPeople reach for the Taylor series because it’s what they remember from calculus, but there’s the Remez algorithm for finding the polynomial which minimizes the maximum error. 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/calc.py https://github.com/depp/ultrafxr/blob/master/math/coeffs/cal...
- 0x20cowboy 5y agoI haven’t seen this version mentioned in the thread - if you don’t need a lot of precision, here is a simple 4 line version[1] and here’s how it works[2]. Not sure who initially came up with it. [1] https://github.com/robrohan/wefx/blob/1a918cc2d5ad87402a383083c4af286d1fcd16f5/src/math.c#L63 https://github.com/robrohan/wefx/blob/1a918cc2d5ad87402a3830... [2] https://www.desmos.com/calculator/lo7cf60mjz https://www.desmos.com/calculator/lo7cf60mjz