Re: sin_pi and cos_pi

[Gerard Jungman <jungman@lanl.gov>]

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.