add sinpi cospi functions

[simonbyrne]

Adds functions for accurately computing sin(pi*x) and cos(pi*x) using range reduction: the main benefit is that sinpi(0.5) == 1.0 and sinpi(1.0) == 0.0. Used internally in sinc function and several other places.

[stevengj]

+1

(In FFTW, we had to implement something much like this, albeit only for rational arguments.)

[simonbyrne]

Ah, I didn't think of rationals: I'll just tweak it a bit for those.

[simonbyrne]

Unfortunately I can't get sinpi(1//6) == 0.5, due to the rounding in the multiplication. It's possible to fix this (at least for this case) using the method I outlined in #4114, but it might not be worth the additional overhead.

[StefanKarpinski]

I'm quite sympathetic to having the correct expensive versions of these sorts of things and working hard to make sure we can get them as fast as possible. In general there tends to be an obvious imperfect but close and fast version – e.g. just using sin(2pi*x/360) in place if sind – the hard part is getting it absolutely right and as fast as possible. Eight flops sounds like a lot but if used in a loop on a super scalar processor (i.e. all modern ones), it actually often ends up only being a couple of cycles per iteration, and outside of loops, you'd almost always rather have accuracy than speed.

[simonbyrne]

The range reduction for sind and sinpi is the low hanging fruit here, as it considerably improves relative accuracy for operations like sinpi(1.0) (all cases I've tried have given relative accuracy to 2eps). The overhead seems to be about an extra 50% vs the simple multiplication, but it might be possible to optimise this further.

The rest is tricky, because not only do we have the rounding of the multiplication, we also have the rounding of sin itself, which can interact in complicated ways. For instance:

julia> sin(pi/6)
0.49999999999999994

julia> float64(sin(big(pi)/6))
0.5

julia> sin(float64(big(pi)/6))
0.5

Yet if we round 1/6 first:

julia> sin(pi*(1/6))
0.49999999999999994

julia> float64(sin(big(pi)*(1/6)))
0.5

julia> sin(float64(big(pi)*(1/6)))
0.49999999999999994

[StefanKarpinski]

Well, just let me know and I'll merge this if you feel it's ready.

[simonbyrne]

I'm happy with it as it stands, but others may have opinions as to whether it's worth inserting a clause to handle the sinpi(1//6) == sind(30) == 0.5 cases exactly (as these have exact Float64 answers).

The only other thing that could be worth thinking about is defining Degree and PiMultiple types to allow overloading of sin and cos instead of using extra functions, but I feel that's probably an argument for another day.

[StefanKarpinski]

Ok, merging.

[stevengj]

Can we please stop merging patches that don't update the documentation as needed? It would be good to institute this as a general rule.

[StefanKarpinski]

Yes, that's a fair point. @simonbyrne, would you please make another PR with docs?