Hi,
I'm attaching a reduced testcase that can be compiled. I'm getting a
negative result in the first iteration.
The reason appears to be that the macro GSL_MAX is defined like this in
gsl/gsl_math.h:
#define GSL_MAX(a,b) ((a) > (b) ? (a) : (b))
The preprocessor therefore translates the line
double max = GSL_MAX (0.0, gsl_ran_gaussian (r, 2));
into
double max = (0.0 > gsl_ran_gaussian (r, 2) ? 0.0 : gsl_ran_gaussian (r,
2));
Because gsl_ran_gaussian () does not always return the same value for the
same input arguments, it is possible that the first function call yields
something positive, so the ( ? : ) expression decides not to return 0.0, but
the result of a second call to the function, which might be negative.
I can't see any obvious way to fix this right now (except for replacing the
macro by a function call). Maybe one could say in the documentation that one
should use the workaround
double random = gsl_ran_gaussian (r, 2);
double max = GSL_MAX (0.0, random);
in cases like this.
Cheers,
Frank
2008/6/9 K Okamoto <[EMAIL PROTECTED]>:
> Hello,
>
> I searched gsl_max in the bug archives to no avail so I am guessing
> this is the first time? In any event I found GSL_MAX to occasionally
> give the wrong answer. Here's the code:
>
/* compile and link with
gcc -o testcase testcase.c -lgsl -lgslcblas -lm
*/
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
int main () {
const gsl_rng_type * T = gsl_rng_default;
gsl_rng * r = gsl_rng_alloc (T);
gsl_rng_set(r , 0);
int IT;
for (IT = 0; IT < 1; IT++)
{
double max = GSL_MAX (0.0, gsl_ran_gaussian (r, 2));
printf ("%g\n", max);
}
}
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
