public inbox for gsl-discuss@sourceware.org
 help / color / mirror / Atom feed
* 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

* Re: sin_pi and cos_pi
  2017-06-28 20:51 ` 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

* 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

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 --
2017-06-27 19:22 sin_pi and cos_pi Gerard Jungman
     [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

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).