* Re: sin_pi and cos_pi [not found] <PsQhDk0sE8yDv7AdpVknc-L90U5KXCOImz_GJaw_iVm9q9Inq1YdQkE4g0BoWmIAvt_E_NoObbeLb8bbWgkLbrXm-S9J-E-2NnaJbYja63g=@protonmail.com> @ 2017-06-28 20:51 ` Gerard Jungman 2017-06-29 7:53 ` Patrick Alken 0 siblings, 1 reply; 3+ messages in thread From: Gerard Jungman @ 2017-06-28 20:51 UTC (permalink / raw) To: Konrad, gsl-discuss, alken On 06/28/2017 04:51 AM, Konrad wrote: > Dear Gerard and all, > > thanks for your comments, I wasn't aware of the caveats concerning > round() but your implementation takes care of that. Like I say, I do not know if what I did was the best choice, but I do know that round(), and friends like nearbyint(), seemed to be problematic. > I'll try to combine our codes if that's ok with you That's fine. Whatever works. Obviously, any improvements are welcome, since I would like to know what I missed. > in particular I'd like to provide the > error estimate since this is (for me) one the main reasons to use GSL. I will comment about error estimates below, since that is something that I have thought about in the intervening years. > There are a few questions remaining before I proceed: > > why are you reducing the argument to [-1,1] but not to [-0.5,0.5]? > This would avoid one of the special cases you treat explicitly, > parameterized by 1-z. That may be only a matter of style. You have an explicit test for fabs(fracx) > 0.5 which flips it around. I like to avoid the extra flip and treat the case directly. I find flips like that confusing. In the end, I think we have the same number of branches, by direct count. The branch patterns will be different though, so it is hard to compare. Code size is different too. > in the GSL implementation do we actually need the special case for > fabs(z) around 0.5 > As far as I know this special case is implemented > reliably in gsl_sf_sin_e and gsl_sf_cos_e which can be used for the > intermediate cases as well. It took me a while to come up with the right comment here. In the end, I think the problem I have is with the idea of deferring to the GSL implementations at all. There are a number of reasons not to do this. First, there is logic and machinery in gsl_sf_sin_e() and gsl_sf_cos_e() for handling large arguments. That effort is wasted, since we know that we have a reduced argument, so why not avoid the cost of invoking the full implementation. You also avoid the headache of introducing some possible problems that could be lurking in the GSL functions. When evaluating a function near a zero, or other special point, I like to use a method which manifestly converges to the exact answer, without reliance on the logic in some other function, even if I wrote that function myself... If I had more patience, I would replace the invocation of the library functions with direct evaluations, such as a rational approximation. But I did not get around to trying that. The series that I do use are kind of lame, but they do satisfy the minimum guarantee of a smooth transition. I could make some further technical comments about series, but that would lead in another direction. Second, the GSL functions can never compete with the library implementations. They were never meant to; instead they were built to produce some kind of error estimate, under some circumstances. I now consider that motivation to be wrong; see the discussion below about error estimates. Of course, the library functions expend unnecessary effort doing their argument reduction gymnastics too. Ideally, I would like to have a gold-plated implementation for arguments which are already reduced, which can be invoked in situations like this. But the gymnastics in the library functions are highly optimized, so they bother me less. > I'll check out the glibc implementation as well, but I guess it lacks > error estimates, so as a first step I'd like to stick with the reduction > to our GSL-trig functions and look into optimizing it with the glibc > methods afterwards. Let me know what you think of this strategy. This leads into the discussion of error estimates. Decades ago, there was a real concern about issues like argument reduction and the level of precision in library functions like sin() and cos(). There may still be some issues on weird platforms, or with certain compiler flags on some platforms. But people who work with such things are aware of the issues and act accordingly. For practical concerns, we can assume IEEE 754 or some equivalent. In this case, the problems for most (probably all) elementary functions were solved a long time ago. For example, the argument reduction problem was solved at Sun in the early 1990s. You can read the short research paper by Ng, "Argument Reduction: Good to the Last Bit". It's all over the web. Here is one link: https://www.csee.umbc.edu/~phatak/645/supl/Ng-ArgReduction.pdf The code in glibc is precisely the code written by this group at SunPro. It still carries the Sun copyright notice, 1993. It basically guarantees that sin() and cos() are correct to the last bit. The error in their private implementation of sin_pi() will be no worse. Obviously that case is much easier than reducing mod PI; as you see for yourself, you basically just need modf() for sin_pi(). So there is essentially no error in sin() and cos(), or sin_pi(). More precisely, you can assign an error which is equivalent to one ulp. So there is really not much point in computing an error. And certainly not the way that the GSL implementations do it. The culture of that implementation really stems from the funky Fortran world of the 1970s, when carpets were orange-colored shag, and computer consoles had built-in ashtrays. So my advice is to avoid the GSL trig implementations and use the standard library. If you need an error to propagate into something, you can assume the error on sin() and cos() is one ulp. In particular, since the GSL interfaces require returning an error, I suggest just using the library sin() and cos() and returning one ulp, basically a fractional error of DBL_EPSILON or so. My general feelings about error estimation are as follows. First, if a function is implemented in such a way as to guarantee precision to the last bit, then error estimates are just a waste. The elementary functions in the standard library can be assumed to have this behavior. Second, if a function cannot be made this precise without recourse to crazy machinery (like deferral to a multi-prec implementation), then that limitation is probably just part of the interface, and should be accepted as such, as long as the result produced is reasonable. For example, if you lose 10 bits in a complicated special function, I may not cry. Of course, such limitations should be carefully measured, circumscribed, and documented. This seems to be the case for some library implementations in certain contexts. For example, the complex functions, like csin() in the C99 standard, are not "guaranteed" to be good to the last bit. I think some vendors could not be bothered to figure out all the cases; there is a large range of behavior in a dopey little function like csin(). I am not sure of the status in glibc; I would be interested to know. I think there is explicit language somewhere in the glibc manual that says something like "well, we do our best, but there are no guarantees of bit-accuracy". Probably somewhere around the discussion of complex-valued functions or arithmetic. Yes, even complex arithmetic is a problem for some vendors... Third, if a function cannot be implemented without catastrophic loss of precision in some regime, then it probably should not be implemented at all. Fixing such a problem is a research project, but it should not be inflicted on unknowing users. GSL suffers somewhat from this last problem. It was partly a utility and partly a research project. As such, the quality of the special functions varies a great deal. Some of those functions should really not be in there at all. I can say this, because it is my fault. What I really wanted from the error estimates was a way to characterize the behavior of the implementations. In some sense, I wanted them to be able to cry for help when they were stretched too far. Ideally, once you knew that everything was ok for your application, you could then turn them off and get bare-metal performance. That second part of the dream never came close to realization. Nowadays, I think the right procedure is to characterize the functions from outside. You want to prove that the functions satisfy a minimum specification by a combination of scanning values and a little theorem-proving. This is also a very enlightening practice in general. You get to see where all the weird jumps and hiccups are, just by staring at a scatter plot of value scans. -- G. Jungman By the way, for some reason I seem to not be getting the gsl-discuss emails anymore. Not sure why. Maybe institutional spam filters or something. But if I miss some of this discussion, that may be the reason. ^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: sin_pi and cos_pi 2017-06-28 20:51 ` sin_pi and cos_pi Gerard Jungman @ 2017-06-29 7:53 ` Patrick Alken 0 siblings, 0 replies; 3+ messages in thread From: Patrick Alken @ 2017-06-29 7:53 UTC (permalink / raw) To: Gerard Jungman, Konrad, gsl-discuss Hi Gerard, Interesting post. In the recent release (2.4) I had tried to fix bug #45730 (https://savannah.gnu.org/bugs/index.php?45730), where the bessel_j2 function failed for large input arguments. This was found to be due to calls to the gsl_sf_sin and gsl_sf_cos functions with large inputs, so those functions definitely have bugs for large arguments. Anyway I switched the gsl_sf_sin/cos calls to libm sin/cos which fixed the issue, however I had no clue how to assign the correct error to the sin/cos calls in the j2 function, largely because I have no idea how the SF errors work anyway. I ended up setting the error to 0, which according to your post below seems to be the correct approach. But I remember there were still some issues with the tests, since the SF test suite does some complex tests on these error values, and it didn't like the j2 error result for large inputs. After some tinkering I eventually got the tests to pass ok. I wonder if we should just replace all gsl_sf_sin/cos with sin/cos ? Also what is "ulp" ? Patrick On 06/28/2017 10:51 PM, Gerard Jungman wrote: > On 06/28/2017 04:51 AM, Konrad wrote: >> Dear Gerard and all, >> >> thanks for your comments, I wasn't aware of the caveats concerning >> round() but your implementation takes care of that. > > Like I say, I do not know if what I did was the best choice, but I > do know that round(), and friends like nearbyint(), seemed to be > problematic. > >> I'll try to combine our codes if that's ok with you > > That's fine. Whatever works. Obviously, any improvements > are welcome, since I would like to know what I missed. > >> in particular I'd like to provide the error estimate since this is >> (for me) one the main reasons to use GSL. > > I will comment about error estimates below, since that is something > that I have thought about in the intervening years. > >> There are a few questions remaining before I proceed: >> >> why are you reducing the argument to [-1,1] but not to [-0.5,0.5]? >> This would avoid one of the special cases you treat explicitly, >> parameterized by 1-z. > > That may be only a matter of style. You have an explicit test > for fabs(fracx) > 0.5 which flips it around. I like to avoid > the extra flip and treat the case directly. I find flips > like that confusing. In the end, I think we have the same > number of branches, by direct count. The branch patterns > will be different though, so it is hard to compare. Code > size is different too. > >> in the GSL implementation do we actually need the special case for >> fabs(z) around 0.5 >> As far as I know this special case is implemented reliably in >> gsl_sf_sin_e and gsl_sf_cos_e which can be used for the intermediate >> cases as well. > > It took me a while to come up with the right comment here. In the > end, I think the problem I have is with the idea of deferring to > the GSL implementations at all. There are a number of reasons > not to do this. > > First, there is logic and machinery in gsl_sf_sin_e() and > gsl_sf_cos_e() for handling large arguments. That effort is > wasted, since we know that we have a reduced argument, so > why not avoid the cost of invoking the full implementation. > You also avoid the headache of introducing some possible > problems that could be lurking in the GSL functions. > > When evaluating a function near a zero, or other special point, > I like to use a method which manifestly converges to the exact > answer, without reliance on the logic in some other function, > even if I wrote that function myself... > > If I had more patience, I would replace the invocation of the > library functions with direct evaluations, such as a rational > approximation. But I did not get around to trying that. The > series that I do use are kind of lame, but they do satisfy > the minimum guarantee of a smooth transition. I could > make some further technical comments about series, but > that would lead in another direction. > > Second, the GSL functions can never compete with the library > implementations. They were never meant to; instead they were > built to produce some kind of error estimate, under some > circumstances. I now consider that motivation to be > wrong; see the discussion below about error estimates. > > Of course, the library functions expend unnecessary effort doing > their argument reduction gymnastics too. Ideally, I would like to > have a gold-plated implementation for arguments which are already > reduced, which can be invoked in situations like this. But the > gymnastics in the library functions are highly optimized, so > they bother me less. > > >> I'll check out the glibc implementation as well, but I guess it lacks >> error estimates, so as a first step I'd like to stick with the >> reduction to our GSL-trig functions and look into optimizing it with >> the glibc methods afterwards. Let me know what you think of this >> strategy. > This leads into the discussion of error estimates. Decades ago, there > was a real concern about issues like argument reduction and the level > of precision in library functions like sin() and cos(). There may > still be some issues on weird platforms, or with certain compiler > flags on some platforms. But people who work with such things are > aware of the issues and act accordingly. > > For practical concerns, we can assume IEEE 754 or some equivalent. In > this case, the problems for most (probably all) elementary functions > were solved a long time ago. For example, the argument reduction > problem was solved at Sun in the early 1990s. You can read the > short research paper by Ng, "Argument Reduction: Good to the Last Bit". > It's all over the web. Here is one link: > > https://www.csee.umbc.edu/~phatak/645/supl/Ng-ArgReduction.pdf > > The code in glibc is precisely the code written by this group at SunPro. > It still carries the Sun copyright notice, 1993. It basically guarantees > that sin() and cos() are correct to the last bit. The error in their > private implementation of sin_pi() will be no worse. Obviously that > case is much easier than reducing mod PI; as you see for yourself, > you basically just need modf() for sin_pi(). > > So there is essentially no error in sin() and cos(), or sin_pi(). > More precisely, you can assign an error which is equivalent to one ulp. > So there is really not much point in computing an error. And certainly > not the way that the GSL implementations do it. The culture of that > implementation really stems from the funky Fortran world of the 1970s, > when carpets were orange-colored shag, and computer consoles had > built-in ashtrays. > > So my advice is to avoid the GSL trig implementations and use the > standard library. If you need an error to propagate into something, > you can assume the error on sin() and cos() is one ulp. > > In particular, since the GSL interfaces require returning an error, > I suggest just using the library sin() and cos() and returning one > ulp, basically a fractional error of DBL_EPSILON or so. > > > My general feelings about error estimation are as follows. > First, if a function is implemented in such a way as to > guarantee precision to the last bit, then error estimates > are just a waste. The elementary functions in the standard > library can be assumed to have this behavior. > > Second, if a function cannot be made this precise without > recourse to crazy machinery (like deferral to a multi-prec > implementation), then that limitation is probably just part > of the interface, and should be accepted as such, as long > as the result produced is reasonable. For example, if you > lose 10 bits in a complicated special function, I may not > cry. Of course, such limitations should be carefully > measured, circumscribed, and documented. > > This seems to be the case for some library implementations in > certain contexts. For example, the complex functions, like csin() > in the C99 standard, are not "guaranteed" to be good to the last > bit. I think some vendors could not be bothered to figure out > all the cases; there is a large range of behavior in a dopey > little function like csin(). I am not sure of the status in > glibc; I would be interested to know. > > I think there is explicit language somewhere in the glibc manual > that says something like "well, we do our best, but there are > no guarantees of bit-accuracy". Probably somewhere around the > discussion of complex-valued functions or arithmetic. Yes, > even complex arithmetic is a problem for some vendors... > > Third, if a function cannot be implemented without catastrophic > loss of precision in some regime, then it probably should not > be implemented at all. Fixing such a problem is a research > project, but it should not be inflicted on unknowing users. > > GSL suffers somewhat from this last problem. It was partly a > utility and partly a research project. As such, the quality > of the special functions varies a great deal. Some of those > functions should really not be in there at all. I can say > this, because it is my fault. > > > What I really wanted from the error estimates was a way to > characterize the behavior of the implementations. In some > sense, I wanted them to be able to cry for help when they > were stretched too far. Ideally, once you knew that everything > was ok for your application, you could then turn them off and > get bare-metal performance. That second part of the dream > never came close to realization. > > Nowadays, I think the right procedure is to characterize the > functions from outside. You want to prove that the functions > satisfy a minimum specification by a combination of scanning > values and a little theorem-proving. This is also a very > enlightening practice in general. You get to see where > all the weird jumps and hiccups are, just by staring > at a scatter plot of value scans. > > -- > G. Jungman > > > By the way, for some reason I seem to not be getting > the gsl-discuss emails anymore. Not sure why. Maybe > institutional spam filters or something. But if I > miss some of this discussion, that may be the reason. > ^ permalink raw reply [flat|nested] 3+ messages in thread
* sin_pi and cos_pi @ 2017-06-27 19:22 Gerard Jungman 0 siblings, 0 replies; 3+ messages in thread From: Gerard Jungman @ 2017-06-27 19:22 UTC (permalink / raw) To: gsl-discuss, alken [-- Attachment #1: Type: text/plain, Size: 3392 bytes --] > After Gerards suggestion I've implemented (see attached) the routines > sin_pi and cos_pi. Essentially, modf and fmod do the trick. I'm a bit > surprised at how simple the solution is, so please check if you see > any problems, but the numerical tests all pass. > > Best, > Konrad Since this is really a design and implementation discussion, and not a bug discussion, I wanted to move it to gsl-discuss. Also, I have no idea how to interact with the bug forum... Anyway, I do have some comments about the method. modf() is almost certainly the right choice for the basic argument reduction. But the combination of round() and fmod() in the "sign" step is probably not ideal. Unfortunately, there are something like a dozen different functions for rounding in the standard library, and some of them have surprising consequences. It's a bit of a minefield. round() is one of the nasty surprise functions in the set. The specification requires it to round away from zero, regardless of the current rounding direction. That sounds like a nice and easy-to-understand choice, until you learn that the way it does this is by _always_ resetting the rounding direction before pushing the instructions into the pipeline. Typically (and certainly on x86-ish cpus) this requires flushing the pipeline before and after the operation. Very expensive. I was looking at the C99 standard rounding functions a few months ago, and some tests showed that round() can easily run twenty times slower than other choices in the set. Some of the functions are specified to not raise an FE_ type hardware exception. Again, this sounds ok, until you learn that this is done by disabling floating-point exceptions before executing the instructions. This operation itself is expensive, combined with the fact that it also has to flush the pipeline. Basically, any of the functions which require changing the state of the floating-point unit will probably have terrible performance implications. Amongst all the choices, testing showed that rint() was relatively harmless. modf() is also ok, as far as I can tell. Attached is an implementation that I was working on earlier this year. It is probably far from perfect. But maybe it has better performance. I would be happy if anybody could prove me wrong. But you'll need to get out your measurement tools and read the docs for all those dopey rounding functions... In any case, any definitive statement about the correct choices here would be a real contribution. I am still unsure about most of it. Finally, if you look deep inside the glibc implementation, you will find an implementation of sin_pi. It is hidden in the gamma function implementation where it is used in the reflection formula. It is not exported, since it is not a standard function. These functions were written at some Sun subsidiary in the early 1990s. They are probably very good, in context. But they are also written in this kind of internal library implementation gibberish style, which is impossible to understand and impossible to export for use elsewhere. You can find one in the file e_lgamma_r.c, for example. There is more than one, because of the support for different platforms. Still, it might be worthwhile to evaluate this. With some work, it might be possible to lift it out of its strange little glibc ghetto and put it to use elsewhere. -- G. Jungman [-- Attachment #2: trig-pi.tar.gz --] [-- Type: application/gzip, Size: 2174 bytes --] ^ permalink raw reply [flat|nested] 3+ messages in thread
end of thread, other threads:[~2017-06-29 7:53 UTC | newest] Thread overview: 3+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- [not found] <PsQhDk0sE8yDv7AdpVknc-L90U5KXCOImz_GJaw_iVm9q9Inq1YdQkE4g0BoWmIAvt_E_NoObbeLb8bbWgkLbrXm-S9J-E-2NnaJbYja63g=@protonmail.com> 2017-06-28 20:51 ` sin_pi and cos_pi Gerard Jungman 2017-06-29 7:53 ` Patrick Alken 2017-06-27 19:22 Gerard Jungman
This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).