From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 33565 invoked by alias); 28 Jun 2017 20:51:32 -0000 Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org Received: (qmail 33545 invoked by uid 89); 28 Jun 2017 20:51:31 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-1.8 required=5.0 tests=AWL,BAYES_00,RP_MATCHES_RCVD,SPF_PASS autolearn=ham version=3.3.2 spammy=feelings, regime, decades, HTo:D*protonmail.com X-HELO: proofpoint8.lanl.gov Received: from proofpoint8.lanl.gov (HELO proofpoint8.lanl.gov) (204.121.3.47) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with ESMTP; Wed, 28 Jun 2017 20:51:29 +0000 Received: from pps.filterd (proofpoint8.lanl.gov [127.0.0.1]) by proofpoint8.lanl.gov (8.16.0.20/8.16.0.20) with SMTP id v5SKlWTm145158; Wed, 28 Jun 2017 14:51:26 -0600 Received: from mailrelay2.lanl.gov (mailrelay2.lanl.gov [128.165.4.103]) by proofpoint8.lanl.gov with ESMTP id 2bbcvm92du-1; Wed, 28 Jun 2017 14:51:26 -0600 Received: from localhost (localhost.localdomain [127.0.0.1]) by mailrelay2.lanl.gov (Postfix) with ESMTP id 388FCE8C03E; Wed, 28 Jun 2017 14:51:26 -0600 (MDT) X-NIE-2-Virus-Scanner: amavisd-new at mailrelay2.lanl.gov Received: from nostromo.lanl.gov (nostromo.lanl.gov [130.55.124.72]) by mailrelay2.lanl.gov (Postfix) with ESMTP id 0365EE8C031; Wed, 28 Jun 2017 14:51:26 -0600 (MDT) Subject: Re: sin_pi and cos_pi To: Konrad , "gsl-discuss@sourceware.org" , "alken@colorado.edu" References: From: Gerard Jungman Message-ID: Date: Wed, 28 Jun 2017 20:51:00 -0000 User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:52.0) Gecko/20100101 Thunderbird/52.0 MIME-Version: 1.0 In-Reply-To: Content-Type: text/plain; charset=utf-8; format=flowed Content-Transfer-Encoding: 7bit X-Proofpoint-Virus-Version: vendor=fsecure engine=2.50.10432:,, definitions=2017-06-28_12:,, signatures=0 X-Proofpoint-Spam-Details: rule=notspam policy=default score=0 spamscore=0 suspectscore=0 malwarescore=0 phishscore=0 adultscore=0 bulkscore=0 classifier=spam adjust=0 reason=mlx scancount=1 engine=8.0.1-1703280000 definitions=main-1706280332 X-SW-Source: 2017-q2/txt/msg00001.txt.bz2 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.