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

[Brian Gough]

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