Hacker News Board Hacker News, browsed like a classic message board
Live data from Hacker News

Implementing Cosine in C from Scratch (2020)

austinhenley.com

Post reply on HN
Print view View original

121–130 of 139 posts

Re: Implementing Cosine in C from Scratch (2020)

#121

Post by bertr4nd »

dancsi wrote: post #12

The 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 Cheby…

Are there libraries/tools that people use to do Remez/Chebyshev/etc. function expansions? I can do a basic Taylor series expansion by hand but I’m out of my depth with more sophisticated techniques.

Re: Implementing Cosine in C from Scratch (2020)

#122

Post by jeremysalwen »

dancsi wrote: post #12

The 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 Cheby…

Thanks for explaining this. I actually wrote a SIMD implementation of trig functions years ago, using the techniques you describe.

You can check it out: https://github.com/jeremysalwen/vectrig

I compared several different methods of generating polynomials of different sizes for speed and precision (spoilers: taylor series were the worst and minimax polynomials (Remez algorithm) were the best).

Another (surprising) thing which I learned during the project was that the range reduction was just as (if not more) important to the accuracy of the implementation than the polynomial. If you think about it, you will realize that it's actually pretty difficult to quickly and accurately compute the sin of large numbers like 2^50.

I also tried to directly optimize the coefficients for the accuracy of the polynomial on the required range, but that experiment was unsuccessful.

It's all there in the repository, the implementations, notes about the different polynomials used, and the accuracy/speed statistics for the different methods.

Re: Implementing Cosine in C from Scratch (2020)

#123

Post by adgjlsfhk1 »

snicker7 wrote: post #115
adgjlsfhk1 wrote: post #39

Earlier quoted context omitted.

Julia 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/03af78108afe6058f54c... for anyone interested.

In 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.

That's how we nerd snipe people into contributing to Base (only 10% joking).

Re: Implementing Cosine in C from Scratch (2020)

#125

Post by dietrichepp »

People 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/cal...

Re: Implementing Cosine in C from Scratch (2020)

#126

Post by eesmith »

aj7 wrote: post #64

So 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 loc…

Going 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

Re: Implementing Cosine in C from Scratch (2020)

#127

Post by stdbrouw »

drBonkers wrote: post #118
aaaronic wrote: post #18

Mandatory 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 ) It's actually kind of ridiculous how many times it comes up and works well enough in undergraduate Mechanics.

How small is small enough to use this approximation?

About 10°, depending on what kind of accuracy you'd like.

https://en.wikipedia.org/wiki/Small-angle_approximation#Erro...

Re: Implementing Cosine in C from Scratch (2020)

#128

Post by m1kal »

nwallin wrote: post #99
phonebucket wrote: post #78

This is a fun article! Alternative avenues that complement the approaches shown: - Padé Approximants ( https://en.wikipedia.org/wiki/Pad%C3%A9_approximant ) can be better than long Taylor Series for this kind of thing. - Bhaskara I's sin approximation ( https://en.wikipedia.org/wiki/Bhaskara_I%27s_sine_approximat... ) is easily adaptable to cosine, remarkably accurate for its simplicity and also fast to calculate.

I 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?

I learned from Bender, Orszag "Advanced Mathematical Methods for Scientists and Engineers I".

Re: Implementing Cosine in C from Scratch (2020)

#129

Post by gpderetta »

adgjlsfhk1 wrote: post #85
stncls wrote: post #63

Earlier quoted context omitted.

Yes 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!

My 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…

The trig call might just be a small component of a piece of cide that is called very often and wants the cache for its own needs.

I think that the grandparent concerns are very legitimate.

Re: Implementing Cosine in C from Scratch (2020)

#130

Post by ramshorns »

adgjlsfhk1 wrote: post #83
enriquto wrote: post #73

Earlier quoted context omitted.

> 0 correct digits when the output is close to 0 this is an amusing way to describe the precision of sub-normal floating point numbers

It's not just sub-normal numbers. As https://randomascii.wordpress.com/2014/10/09/intel-underesti... shows, fsin only uses 66 bits of pi, which means you have roughly no precision whenever abs(sin(x))<10^-16 which is way bigger than the biggest subnormal (10^-307 or so)

In that range, just returning x would be way better. Maybe even perfect actually - if x is less than 10^-16, then the error of x^3/6 is less than the machine precision for x.
Post reply on HN