jhjourdan / SIMD-math-prims

Vectorizable implementations of some mathematical functions
MIT License
102 stars 12 forks source link

Source for the sin and cos formulas #5

Open Sahnvour opened 3 years ago

Sahnvour commented 3 years ago

Hello,

Thanks for your library. I was wondering if you could document how you determined the polynomials for sin and cos.
I tried using sollya with my limited knowledge to find some myself, and always ended up with terrible accuracy over [-pi;pi]. Other sources usually pick a more reduced range, which implies more work outside the approximation to reconstruct the final result.

Also, wouldn't it be better to use fpminimax instead of remez to find coefficients exactly representable in IEEE 754 ?

PS : your blog seems down since a few days https://gallium.inria.fr/blog/fast-vectorizable-math-approx/

jhjourdan commented 3 years ago

Thanks for your library. I was wondering if you could document how you determined the polynomials for sin and cos.

As documented in the source code:

  /* Generated in Sollya using:
     > f = remez(cos(x)-1, [|x*x, x*x*x*x, x*x*x*x*x*x, x*x*x*x*x*x*x*x|],
                            [0.000001, pi], 1, 1e-8);
     > plot(f-cos(x)+1, [0, pi]);
     > f+1
  */

And similar for sin.

I tried using sollya with my limited knowledge to find some myself, and always ended up with terrible accuracy over [-pi;pi].

Could you tell a bit more about what you did in Sollya exactly?

Other sources usually pick a more reduced range, which implies more work outside the approximation to reconstruct the final result.

Indeed, this is a possibility. But then you have to do extra computation, which either involve division (slow) or control flow (incompatible with SIMD, and sometimes hard to predict). Hence, increasing the degree of the polynomila might be a good idea instead.

PS : you blog seems down since a few days https://gallium.inria.fr/blog/fast-vectorizable-math-approx/

Right, there seems to be a problem in the HTTPS configuration of the server. This should be fixed soon. You can still access the blog using HTTP.

Sahnvour commented 3 years ago

Oh right, somehow I only looked at sin and missed it in cos, sorry.

Could you tell a bit more about what you did in Sollya exactly?

I tried to use fpminimax in [-pi;pi] but now I see you only used [0;pi]. I think I understand we can limit the approximation to that interval and mirror it thanks to the properties of even/odd functions, so that makes sense.

However I'm not sure why you use the -1 thing, if you could explain the trick.

jhjourdan commented 3 years ago

However I'm not sure why you use the -1 thing, if you could explain the trick.

That's completely orthogonal: one interesting property for a cos approximation is to be exactly equal to 1. at zero. The trick here is just a way to force the constant term of the polynomial is 1.