[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] Big errors in gsl_cdf_gaussian_P when rounding down
From: |
Abu Yoav |
Subject: |
Re: [Bug-gsl] Big errors in gsl_cdf_gaussian_P when rounding down |
Date: |
Wed, 15 Sep 2010 14:24:21 -0700 |
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
>
>