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

[Abu Yoav]

Hello,

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

A test program with corresponding output follows. Hope this helps...

--- program ---
// compilation: gcc gslbug.c -lgsl -lgslcblas

#include <stdio.h>
#include <gsl/gsl_cdf.h>
#include <math.h>
#include <fenv.h>
#include <assert.h>

int main(int argc, char **argv) {

    double sigma = sqrt(0.1);
    double x1 = -1.146892;
    double x2 = -1.111964;

    double x1cdf, x2cdf;

    printf("With normal rounding\n");

    x1cdf = gsl_cdf_gaussian_P(x1,sigma);
    x2cdf = gsl_cdf_gaussian_P(x2,sigma);

    printf( "x1 = %e, x2 = %e\n", x1, x2);
    printf( "x1cdf = %e, x2cdf = %e\n", x1cdf, x2cdf);

    printf("With round downward\n");

    int res;
    res = fesetround(FE_DOWNWARD); // round downward
    assert(res == 0);

    x1cdf = gsl_cdf_gaussian_P(x1,sigma);
    x2cdf = gsl_cdf_gaussian_P(x2,sigma);

    printf( "x1 = %e, x2 = %e\n", x1, x2);
    printf( "x1cdf = %e, x2cdf = %e\n", x1cdf, x2cdf);


    return 0;
}

--- output ---
With normal rounding
x1 = -1.146892e+00, x2 = -1.111964e+00
x1cdf = 1.434827e-04, x2cdf = 2.187710e-04
With round downward
x1 = -1.146892e+00, x2 = -1.111964e+00
x1cdf = 1.434827e-04, x2cdf = 1.220588e-04