[Numpy-discussion] sinpi/cospi trigonometric functions
Hi all, (I am very new to this mail list so please cut me some slack) trigonometric functions like sin(x) are usually implemented as: 1. Some very complicated function that does bit twiddling and basically computes the reminder of x by pi/2. An example in http://www.netlib.org/fdlibm/e_rem_pio2.c (that calls http://www.netlib.org/fdlibm/k_rem_pio2.c ). i.e. ~500 lines of branching C code. The complexity arises in part because for big values of x the subtraction becomes more and more ill defined, due to x being represented in binary base to which an irrational number has to subtracted and consecutive floating point values being more and more apart for higher absolute values. 2. A Taylor series for the small values of x, 3. Plus some manipulation to get the correct branch, deal with subnormal numbers, deal with -0, etc... If we used a function like sinpi(x) = sin(pi*x) part (1) can be greatly simplified, since it becomes trivial to separate the reminder of the division by pi/2. There are gains both in the accuracy and the performance. In large parts of the code anyways there is a pi inside the argument of sin since it is common to compute something like sin(2*pi*f*t) etc. So I wonder if it is feasible to implement those functions in numpy. To strengthen my argument I'll note that the IEEE standard, too, defines ( https://irem.univ-reunion.fr/IMG/pdf/ieee-754-2008.pdf ) the functions sinPi, cosPi, tanPi, atanPi, atan2Pi. And there are existing implementations, for example, in Julia ( https://github.com/JuliaLang/julia/blob/6aaedecc447e3d8226d5027fb13d0c3cbfbfea2a/base/special/trig.jl#L741-L745 ) and the Boost C++ Math library ( https://www.boost.org/doc/libs/1_54_0/boost/math/special_functions/sin_pi.hpp ) And that issue caused by apparently inexact calculations have been raised in the past in various forums ( https://stackoverflow.com/questions/20903384/numpy-sinpi-returns-negative-value https://stackoverflow.com/questions/51425732/how-to-make-sinpi-and-cospi-2-zero https://www.reddit.com/r/Python/comments/2g99wa/why_does_python_not_make_sinpi_0_just_really/ ... ) PS: to be nitpicky I see that most implementation implement sinpi as sin(pi*x) for small values of x, i.e. they multiply x by pi and then use the same coefficients for the Taylor series as the canonical sin. A multiply instruction could be spared, in my opinion, by storing different Taylor expansion number coefficients tailored for the sinpi function. It is not clear to me if it is not done because the performance gain is small, because I am wrong about something, or because those 6 constants of the Taylor expansion have a "sacred aura" about them and nobody wants to enter deeply into this. PPS: I am aware that it could be seen as rude to request a feature from an open source project but I am asking if there is a point in providing these functions in the first place. I could try to provide implementations for them in some time if it is indeed a worthwhile effort Yours, Tom. ___ NumPy-Discussion mailing list NumPy-Discussion@python.org https://mail.python.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] sinpi/cospi trigonometric functions
On Wed, Jul 14, 2021 at 11:39 AM Tom Programming wrote: > Hi all, > > (I am very new to this mail list so please cut me some slack) > > trigonometric functions like sin(x) are usually implemented as: > > 1. Some very complicated function that does bit twiddling and basically > computes the reminder of x by pi/2. An example in > http://www.netlib.org/fdlibm/e_rem_pio2.c (that calls > http://www.netlib.org/fdlibm/k_rem_pio2.c ). i.e. ~500 lines of branching > C code. The complexity arises in part because for big values of x the > subtraction becomes more and more ill defined, due to x being represented > in binary base to which an irrational number has to subtracted and > consecutive floating point values being more and more apart for higher > absolute values. > 2. A Taylor series for the small values of x, > 3. Plus some manipulation to get the correct branch, deal with subnormal > numbers, deal with -0, etc... > > If we used a function like sinpi(x) = sin(pi*x) part (1) can be greatly > simplified, since it becomes trivial to separate the reminder of the > division by pi/2. There are gains both in the accuracy and the performance. > > In large parts of the code anyways there is a pi inside the argument of > sin since it is common to compute something like sin(2*pi*f*t) etc. So I > wonder if it is feasible to implement those functions in numpy. > > To strengthen my argument I'll note that the IEEE standard, too, defines ( > https://irem.univ-reunion.fr/IMG/pdf/ieee-754-2008.pdf ) the functions > sinPi, cosPi, tanPi, atanPi, atan2Pi. And there are existing > implementations, for example, in Julia ( > https://github.com/JuliaLang/julia/blob/6aaedecc447e3d8226d5027fb13d0c3cbfbfea2a/base/special/trig.jl#L741-L745 > ) and the Boost C++ Math library ( > https://www.boost.org/doc/libs/1_54_0/boost/math/special_functions/sin_pi.hpp > ) > > And that issue caused by apparently inexact calculations have been raised > in the past in various forums ( > https://stackoverflow.com/questions/20903384/numpy-sinpi-returns-negative-value > https://stackoverflow.com/questions/51425732/how-to-make-sinpi-and-cospi-2-zero > https://www.reddit.com/r/Python/comments/2g99wa/why_does_python_not_make_sinpi_0_just_really/ > ... ) > > PS: to be nitpicky I see that most implementation implement sinpi as > sin(pi*x) for small values of x, i.e. they multiply x by pi and then use > the same coefficients for the Taylor series as the canonical sin. A > multiply instruction could be spared, in my opinion, by storing different > Taylor expansion number coefficients tailored for the sinpi function. It is > not clear to me if it is not done because the performance gain is small, > because I am wrong about something, or because those 6 constants of the > Taylor expansion have a "sacred aura" about them and nobody wants to enter > deeply into this. > The main value of the sinpi(x) formulation is that you can do the reduction on x more accurately than on pi*x (reduce-then-multiply rather than multiply-then-reduce) for people who particularly care about the special locations of half-integer x. sin() and cos() are often not implemented in software, but by CPU instructions, so you don't want to reimplement them. There is likely not a large accuracy win by removing the final multiplication. We do have sindg(), cosdg(), and tandg() in scipy.special that work similarly for inputs in degrees rather than radians. They also follow the reduce-then-multiply strategy. scipy.special would be a good place for sinpi() and friends. -- Robert Kern ___ NumPy-Discussion mailing list NumPy-Discussion@python.org https://mail.python.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] sinpi/cospi trigonometric functions
I'll note that SciPy actually does have `sinpi` and `cospi`-they just happen to be private: https://github.com/scipy/scipy/blob/master/scipy/special/functions.json#L58 https://github.com/scipy/scipy/blob/master/scipy/special/functions.json#L12 They are used extensively inside the module though as helpers in other special functions and have extensive tests of their own: https://github.com/scipy/scipy/blob/master/scipy/special/tests/test_mpmath.py#L533 https://github.com/scipy/scipy/blob/master/scipy/special/tests/test_mpmath.py#L547 https://github.com/scipy/scipy/blob/master/scipy/special/tests/test_mpmath.py#L960 https://github.com/scipy/scipy/blob/master/scipy/special/tests/test_mpmath.py#L1741 I have no objections to making them public; the PR is as simple as removing the underscores and adding a docstring. On Wed, Jul 14, 2021 at 9:17 AM Robert Kern wrote: > > On Wed, Jul 14, 2021 at 11:39 AM Tom Programming > wrote: >> >> Hi all, >> >> (I am very new to this mail list so please cut me some slack) >> >> trigonometric functions like sin(x) are usually implemented as: >> >> 1. Some very complicated function that does bit twiddling and basically >> computes the reminder of x by pi/2. An example in >> http://www.netlib.org/fdlibm/e_rem_pio2.c (that calls >> http://www.netlib.org/fdlibm/k_rem_pio2.c ). i.e. ~500 lines of branching C >> code. The complexity arises in part because for big values of x the >> subtraction becomes more and more ill defined, due to x being represented in >> binary base to which an irrational number has to subtracted and consecutive >> floating point values being more and more apart for higher absolute values. >> 2. A Taylor series for the small values of x, >> 3. Plus some manipulation to get the correct branch, deal with subnormal >> numbers, deal with -0, etc... >> >> If we used a function like sinpi(x) = sin(pi*x) part (1) can be greatly >> simplified, since it becomes trivial to separate the reminder of the >> division by pi/2. There are gains both in the accuracy and the performance. >> >> In large parts of the code anyways there is a pi inside the argument of sin >> since it is common to compute something like sin(2*pi*f*t) etc. So I wonder >> if it is feasible to implement those functions in numpy. >> >> To strengthen my argument I'll note that the IEEE standard, too, defines ( >> https://irem.univ-reunion.fr/IMG/pdf/ieee-754-2008.pdf ) the functions >> sinPi, cosPi, tanPi, atanPi, atan2Pi. And there are existing >> implementations, for example, in Julia ( >> https://github.com/JuliaLang/julia/blob/6aaedecc447e3d8226d5027fb13d0c3cbfbfea2a/base/special/trig.jl#L741-L745 >> ) and the Boost C++ Math library ( >> https://www.boost.org/doc/libs/1_54_0/boost/math/special_functions/sin_pi.hpp >> ) >> >> And that issue caused by apparently inexact calculations have been raised in >> the past in various forums ( >> https://stackoverflow.com/questions/20903384/numpy-sinpi-returns-negative-value >> >> https://stackoverflow.com/questions/51425732/how-to-make-sinpi-and-cospi-2-zero >> >> https://www.reddit.com/r/Python/comments/2g99wa/why_does_python_not_make_sinpi_0_just_really/ >> ... ) >> >> PS: to be nitpicky I see that most implementation implement sinpi as >> sin(pi*x) for small values of x, i.e. they multiply x by pi and then use the >> same coefficients for the Taylor series as the canonical sin. A multiply >> instruction could be spared, in my opinion, by storing different Taylor >> expansion number coefficients tailored for the sinpi function. It is not >> clear to me if it is not done because the performance gain is small, because >> I am wrong about something, or because those 6 constants of the Taylor >> expansion have a "sacred aura" about them and nobody wants to enter deeply >> into this. > > > The main value of the sinpi(x) formulation is that you can do the reduction > on x more accurately than on pi*x (reduce-then-multiply rather than > multiply-then-reduce) for people who particularly care about the special > locations of half-integer x. sin() and cos() are often not implemented in > software, but by CPU instructions, so you don't want to reimplement them. > There is likely not a large accuracy win by removing the final multiplication. > > We do have sindg(), cosdg(), and tandg() in scipy.special that work similarly > for inputs in degrees rather than radians. They also follow the > reduce-then-multiply strategy. scipy.special would be a good place for > sinpi() and friends. > > -- > Robert Kern > ___ > NumPy-Discussion mailing list > NumPy-Discussion@python.org > https://mail.python.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@python.org https://mail.python.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] sinpi/cospi trigonometric functions
On Wed, Jul 14, 2021 at 12:26 PM Joshua Wilson wrote: > I'll note that SciPy actually does have `sinpi` and `cospi`-they just > happen to be private: > > https://github.com/scipy/scipy/blob/master/scipy/special/functions.json#L58 > https://github.com/scipy/scipy/blob/master/scipy/special/functions.json#L12 > > They are used extensively inside the module though as helpers in other > special functions and have extensive tests of their own: > > > https://github.com/scipy/scipy/blob/master/scipy/special/tests/test_mpmath.py#L533 > > https://github.com/scipy/scipy/blob/master/scipy/special/tests/test_mpmath.py#L547 > > https://github.com/scipy/scipy/blob/master/scipy/special/tests/test_mpmath.py#L960 > > https://github.com/scipy/scipy/blob/master/scipy/special/tests/test_mpmath.py#L1741 > > I have no objections to making them public; the PR is as simple as > removing the underscores and adding a docstring. > Delightful! -- Robert Kern ___ NumPy-Discussion mailing list NumPy-Discussion@python.org https://mail.python.org/mailman/listinfo/numpy-discussion
