Re: [Bug-gsl] Big errors in gsl_cdf_gaussian_P when rounding down

[Abu Yoav]

Hi again,

It seems this bug is known (for quite some time):
http://sources.redhat.com/bugzilla/show_bug.cgi?id=3976

I'll code a workaround in my program. Thanks again for your help.

On Wed, Sep 15, 2010 at 1:38 PM, Brian Gough <address@hidden> wrote:

> At Tue, 14 Sep 2010 13:32:20 -0700,
> Abu Yoav wrote:
> > I'm not sure that this is -- strictly speaking -- a bug. However, it
> > does seem like something which can happen in a typical use case, and
> > which gives unexpected results.
> >
> > I have a program in which the rounding mode is set to downward:
> > fesetround(FE_DOWNWARD);
> > This seemingly naive setting gave me an error of about 80% when running
> > gsl_cdf_gaussian_P. More so, for sigma fixed and x1<x2, I got
> > gsl_cdf_gaussian_P(x1,sigma)>gsl_cdf_gaussian_P(x2,sigma).
>
> That's interesting, I investigated it further to see where it
> originates from.  It is purely from the exp() function, which gives
> different results with the exact same argument in the two modes, so I
> don't think there is much we can do about it.  The value of u is the
> one occurring in the computation of the gauss cdf.  If you are able
> to, it would be helpful if you could look in the GLIBC bugzilla and
> see if it is a known problem.
>
>
> #include <stdio.h>
> #include <math.h>
> #include <fenv.h>
> #include <assert.h>
>
> int main(int argc, char **argv) {
>
>     double u = -3.267668770400000144e-01, z;
>     printf("With normal rounding\n");
>     z=exp(u); printf("z=%.18e\n", z);
>     printf("With round downward\n");
>
>    int res;
>    res = fesetround(FE_DOWNWARD); // round downward
>    assert(res == 0);
>     z=exp(u); printf("z=%.18e\n", z);
>
>    return 0;
> }
>
>
> $ gcc -Wall fernd.c -lm
> address@hidden:~/tmp$ ./a.out
> With normal rounding
> z=7.212518637988338810e-01
> With round downward
> z=5.005904787886047425e-01
>
>