Dear Gerard and Patrick, thanks for your insightful answers, it is indeed quite tempting to take the standard trigonometric functions for granted. If I find the time maybe as a first step separate functions sin_pi and cos_pi with proper argument reduction could be added and once this is implemented we can come back to this issue.
Cheers, Konrad On 06/16/2017 07:30 AM, Patrick Alken wrote: > Thanks Gerard, its good to know you've thought about this already. I've > added your comments to the bug tracker for this issue. > > On 06/15/2017 09:45 PM, Gerard Jungman wrote: >> On 06/15/2017 03:59 AM, Patrick Alken wrote: >>> >>> recently I needed some of the Bessel functions Y_nu, but noticed that >>> the GSL-implementation raises a domain error for nu < 0. I'm not sure >>> why this is the case, maybe there are good reasons I have overlooked, >>> but to me it seems we could use the rotation between Bessel functions of >>> first and second kind to cover all values of nu. >> >> >> There are three issues. One is mathematical, one is a related >> observation involving the mirror implementation for J, and the >> last is a smallish but important design issue related to the >> form of the proposed change. >> >> First the math issue. Suppose nu is exactly integral. The naive >> evaluation of sin(M_PI*nu) and cos(M_PI*nu) could be very bad, >> since you really want these to be exactly 0.0 and 1.0, to >> avoid strange admixtures of possibly singular values. >> The behavior near integral nu is also suspect, since >> you must validate the behavior of the trig functions >> near zeros. >> >> Therefore, you really need a set of functions which are sometimes >> called "sin_pi(x)" and "cos_pi(x)" which perform the argument >> reduction of x exactly. The built-in argument reduction of M_PI*x >> is not acceptable. >> >> Now for the observation involving J, which is directly related >> to this point. It may be that admixtures of the near-singular >> solution are ok when evaluating Y, since Y dominates any >> solution, but that the problem is much more difficult in >> the case of J. >> >> The difference may be large enough that you can imagine >> easily implementing Y but punting on J. Then you would >> be left with two functions which are mathematically >> similar but for which the exported interfaces have >> different preconditions. This is confusing enough >> that it becomes an argument for leaving them both out. >> >> Now for the smallish design issue. >> >> If these problems can be solved in a coherent way, then >> it would be fine to provide Y and J for negative nu by >> rotation. But the code needs to be refactored. The proposed >> change involves the function calling itself, with flipped argument, >> which is bad levelization. Rather, we need separate functions >> which handle positive nu, which are called by wrappers which >> handle both positive and negative nu. >> >> The main task of the wrappers would be to manage the >> delegation to properly argument-reduced sin_pi and >> cos_pi functions when nu < 0, so their job is not >> entirely trivial after all. >> >> Whether or not we continue to export the nu>0 interfaces >> is a separate question. I don't think I have any special >> feelings about that. >> >> >> I do not remember much about what I was thinking 20 years ago. >> But I can guess that the main problem was this business of >> proper argument reduction. Maybe it is a red herring after >> all, but this would need to be demonstrated. >> >> As a general comment, it is always necessary to be very careful >> about the quality of the underlying trig functions in this kind >> of situation. Sadly, quality will vary from platform to platform. >> >> -- >> G. Jungman > >
