[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] Big errors in gsl_cdf_gaussian_P when rounding down
From: |
Abu Yoav |
Subject: |
[Bug-gsl] Big errors in gsl_cdf_gaussian_P when rounding down |
Date: |
Tue, 14 Sep 2010 13:32:20 -0700 |
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
- [Bug-gsl] Big errors in gsl_cdf_gaussian_P when rounding down,
Abu Yoav <=