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 <[email protected]> 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 > b...@fencepost:~/tmp$ ./a.out > With normal rounding > z=7.212518637988338810e-01 > With round downward > z=5.005904787886047425e-01 > > _______________________________________________ Bug-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/bug-gsl
