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

51–60 of 139 posts

Re: Implementing Cosine in C from Scratch (2020)

#51

Post by duped »

blt wrote: post #32

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

Here'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.

Re: Implementing Cosine in C from Scratch (2020)

#52

Post by jhgb »

lifthrasiir wrote: post #50
jhgb wrote: post #45

Earlier quoted context omitted.

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

The input might be double but the constant pi is not. Let f64(x) be a function from any real number to double, so that an ordinary expression `a + b` actually computes f64(a + b) and so on. Then in general f64(sin(x)) may differ from f64(sin(f64(x mod 2pi))); since you can't directly compute f64(sin(x mod 2pi)), you necessarily need more precision during argument reduction so that f64(sin(x)) = f64(sin(f64timeswhatev…

But am I correct in thinking that is at worst a 0.5 ulp error in this case? The lesser term in double-double can't be more than 0.5 ulp of the greater term and sensitivity of both sine and cosine to an error in the input will not be more than 1.

Also, in case of confusion, I was specifically commenting on the function over the [-pi/4, pi/4] domain in https://github.com/ifduyue/musl/blob/master/src/math/__cos.c , which the comment in https://news.ycombinator.com/item?id=30846546 was presumably about.

Re: Implementing Cosine in C from Scratch (2020)

#53

Post by adgjlsfhk1 »

jhgb wrote: post #52
lifthrasiir wrote: post #50

Earlier quoted context omitted.

The input might be double but the constant pi is not. Let f64(x) be a function from any real number to double, so that an ordinary expression `a + b` actually computes f64(a + b) and so on. Then in general f64(sin(x)) may differ from f64(sin(f64(x mod 2pi))); since you can't directly compute f64(sin(x mod 2pi)), you necessarily need more precision during argument reduction so that f64(sin(x)) = f64(sin(f64timeswhatev…

But am I correct in thinking that is at worst a 0.5 ulp error in this case? The lesser term in double-double can't be more than 0.5 ulp of the greater term and sensitivity of both sine and cosine to an error in the input will not be more than 1. Also, in case of confusion, I was specifically commenting on the function over the [-pi/4, pi/4] domain in https://github.com/ifduyue/musl/blob/master/src/math/__cos.c , whic…

Double rounding can still bite you. You are forced to incur up to half an ulp of error from your polynomial, so taking another half ulp in your reduction can lead to a total error of about 1 ulp.

Re: Implementing Cosine in C from Scratch (2020)

#54

Post by adgjlsfhk1 »

anonymousiam wrote: post #35

I 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 (This was on a…

I'm surprised this was faster than an initial guess (half the exponent and mantisa) followed by newton's method.

Re: Implementing Cosine in C from Scratch (2020)

#55

Post by jhgb »

gugagore wrote: post #47
jhgb wrote: post #20

Earlier quoted context omitted.

This 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).

It'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.

It's not obvious but that makes it a nice example in parameter optimization as well.

Re: Implementing Cosine in C from Scratch (2020)

#56

Post by lifthrasiir »

jhgb wrote: post #52
lifthrasiir wrote: post #50

Earlier quoted context omitted.

The input might be double but the constant pi is not. Let f64(x) be a function from any real number to double, so that an ordinary expression `a + b` actually computes f64(a + b) and so on. Then in general f64(sin(x)) may differ from f64(sin(f64(x mod 2pi))); since you can't directly compute f64(sin(x mod 2pi)), you necessarily need more precision during argument reduction so that f64(sin(x)) = f64(sin(f64timeswhatev…

But am I correct in thinking that is at worst a 0.5 ulp error in this case? The lesser term in double-double can't be more than 0.5 ulp of the greater term and sensitivity of both sine and cosine to an error in the input will not be more than 1. Also, in case of confusion, I was specifically commenting on the function over the [-pi/4, pi/4] domain in https://github.com/ifduyue/musl/blob/master/src/math/__cos.c , whic…

Yeah, sine and cosine are not as sensitive (but note that many libms target 1 or 1.5 ulp error for them, so a 0.5 ulp error might still be significant). For tangent however you definitely need more accurate range reduction.

Re: Implementing Cosine in C from Scratch (2020)

#57

Post by stncls »

This 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/#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

Re: Implementing Cosine in C from Scratch (2020)

#58

Post by enriquto »

brandmeyer wrote: post #44
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…

This has been the standard algorithm used by every libm for decades. Its not special to Musl.

But isn't this code rarely called in practice? I guess on intel architectures the compiler just calls the fsin instruction of the cpu.

Re: Implementing Cosine in C from Scratch (2020)

#59

Post by ufo »

l- wrote: post #15
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…

https://github.com/ifduyue/musl/blob/master/src/math/__cos.c

I think the polynomial calculation in the end looks interesting. It doesn't use Horner's rule.

Re: Implementing Cosine in C from Scratch (2020)

#60

Post by adgjlsfhk1 »

stncls wrote: post #57

This 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/#sm4fjqtvjyn36dyednc... ). Now the method re…

How 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.
Post reply on HN