cos_pi and sin_pi return erroneous sign for some big values

[jtlap]

using boost 1.72 in special functions part of boost math in file cos_pi.hpp lines 36 and 37:

36 if(iconvert(rem, pol) & 1) 37 invert = !invert;

on the call cos_pi( 324812120769207.375) (for instance), rem is 324812120769207 and iconvert(rem,pol) is -2147483648 which is certainly not expected as it changes the parity and so the invert boolean value is not changed and produces a wrong sign for the cosine. (cos_pi(1.375) is computed with the good sign because iconvert of 1 is 1)

Possible correction:

As the only important thing here is the parity of the flint value rem, I suggest to replace the test line 36 by

36 if( floor(rem/2)*2 != rem)

as it does not use any integral conversion

related:

sin_pi.h has exactly the same problem lines 41 and 42

How I discovered the error:

I wanted to use the function to test an simd implementation of mine of cospi and wondered why for some big values my result was the opposite. Using sollya I convinced myself that my result was the good one and found the bug in the quite easy to read implementation in boost math.

I made the patch locally and it corrected all errors in my tests. This is perhaps also a problem to see in iconvert.

Thank you for your very useful library.

J.T. Lapresté ([email protected])

[NAThompson]

@jtlap : Do you mind posting the diff from your patch or making a PR?

Thanks!

[jtlap]

ok, I have a pr ready, but I don't know how to push it as I have no login for 'https://github.com/boostorg/math.git/'

[jtlap]

Is there any guest login available ?

[jtlap]

however here is a diff result:

iff --git a/include/boost/math/special_functions/cos_pi.hpp b/include/boost/math/special_functions/cos_pi.hpp index c565ac08b..d395d262b 100644 --- a/include/boost/math/special_functions/cos_pi.hpp +++ b/include/boost/math/special_functions/cos_pi.hpp @@ -33,7 +33,7 @@ T cos_pi_imp(T x, const Policy& pol) x = -x; } T rem = floor(x);

• if(iconvert(rem, pol) & 1)

• if(floor(rem/2)*2 != rem) invert = !invert; rem = x - rem; if(rem > 0.5f) diff --git a/include/boost/math/special_functions/sin_pi.hpp b/include/boost/math/special_functions/sin_pi.hpp index c55c4ff79..cf4fff9f4 100644 --- a/include/boost/math/special_functions/sin_pi.hpp +++ b/include/boost/math/special_functions/sin_pi.hpp @@ -38,7 +38,7 @@ inline T sin_pi_imp(T x, const Policy& pol) invert = false;

  T rem = floor(x);

• if(iconvert(rem, pol) & 1)

• if(floor(rem/2)*2 != rem) invert = !invert;

[NAThompson]

@jtlap : You have to fork the repo, then the PR must be "compared across forks".

Your diff is sufficient for me to get it working, but of course I'd like you to get a hang of the "fork->PR" workflow so we can continue benefiting from your fixes!

[jtlap]

thank you, I will let you make the fix. I will fork "next time", if needed

[jtlap]

Also, I will not investigate further, but I am quite sure there will be problems with iconvert