1 for (I = problems ; I->f != 0; I++)
2 {
3   size_t i;
4   double res, err;
5 
6   gsl_rng * r;
7 
8   if (I->dim > 3)
9     {
10       continue ;
11     }
12 
13   r = gsl_rng_alloc (gsl_rng_default);
14 
15   for (i = 0; i < TRIALS ; i++)
16     {
17       MONTE_STATE *s = MONTE_ALLOC (I->dim);
18 #ifdef MONTE_PARAMS
19       MONTE_PARAMS params;
20 #endif
21 
22       I->f->dim = I->dim;
23 
24       MONTE_INTEGRATE (I->f, I->xl, I->xu,
25                        I->dim, I->calls / MONTE_SPEEDUP, r, s,
26                        &res, &err);
27 
28       gsl_test_abs (res, I->expected_result,
29                     5 * GSL_MAX(err, 1024*GSL_DBL_EPSILON),
30                     NAME ", %s, result[%d]", I->description, i);
31 
32       MONTE_ERROR_TEST (err, I->expected_error);
33 
34       result[i] = res;
35       error[i] = err;
36 
37       MONTE_FREE (s);
38     }
39 
40 
41   /* Check the results for consistency as an ensemble */
42 
43   {
44     double mean = 0, sumd2 = 0, sd;
45 
46     /* We need to compute the mean exactly when all terms are equal,
47        to get an exact zero for the standard deviation (this is a
48        common case when integrating a constant). */
49 
50     for (i = 0; i < TRIALS; i++)
51       {
52         mean += (result[i] - mean) / (i + 1.0);
53       }
54 
55     for (i = 0; i < TRIALS; i++)
56       {
57         sumd2 += pow(result[i] - mean, 2.0);
58       }
59 
60     sd = sqrt(sumd2 / (TRIALS-1.0)) ;
61 
62     for (i = 0; i < TRIALS; i++)
63       {
64         gsl_test_factor (error[i], sd, 5.0,
65                          NAME ", %s, abserr[%d] vs sd", I->description, i);
66       }
67   }
68 
69   gsl_rng_free (r);
70 }
71 
72