sin_pi and cos_pi

[Gerard Jungman <jungman at lanl dot gov>]

> 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: trig-pi.tar.gz
Description: application/gzip