#include <gsl/gsl_cdf.h>
#include <stdio.h>

#include <gsl/gsl_math.h>

/* tolerance taken from cdf/test.c */
#define TEST_TOL6  (1048576.0*GSL_DBL_EPSILON)

int main () {
  double P [4];  /* P values calculated using gsl_cdf_fdist_P () */
  double eP [4]; /* expected values taken from cdf/test.c */
  
  P [0] = gsl_cdf_fdist_P (3.479082213465832574, 1, 4040712);
  eP [0] = 0.93785072763723411967;

  P [1] = gsl_cdf_fdist_P (3.002774644786533109, 1, 4040712);
  eP [1] = 0.91687787379476055771;

  P [2] = gsl_cdf_fdist_P (3.000854441173130827, 1, 4040712);
  eP [2] = 0.91677930719813578619;

  P [3] = gsl_cdf_fdist_P (3.000064021622133037, 1, 4040712);
  eP [3] = 0.91678021757407962678;

  int i;
  
  for (i = 0; i < 4; i++)
    printf ("(P - eP) [%i] = %.16e, tolerance = %.16e\n",
	    i, eP [i] - P [i], TEST_TOL6);

  return 0;
}