[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] gsl_max and nested random variables
From: |
K Okamoto |
Subject: |
[Bug-gsl] gsl_max and nested random variables |
Date: |
Mon, 9 Jun 2008 11:14:37 +0200 |
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:
int IT;
for (IT = 0; IT < 500; IT++)
{
printf("%f %f\n", GSL_MAX(0, gsl_ran_gaussian(r, 2)), mymax(0,
gsl_ran_gaussian(r, 2)));
}
/* Where mymax is: */
double mymax(double x, double y)
{
if (x >= y)
return(x);
else
return(y);
}
The first element is occasionally < 0, while the second element is
never less than zero. I am pretty sure it is a bug, because this
problem goes away if I have:
for (IT = 0; IT < 500; IT++)
{
double use = gsl_ran_gaussian(r,2);
printf("%f %f\n", GSL_MAX(0, use), mymax(0, gsl_ran_gaussian(r, 2)));
}
using GSL_MAX_DBL doesn't change the outcome.
I am compiling this using gcc v.4.1.1 on fedora 7.
Thx,
ken
PS.
The random number generator is initialized in main() as:
int seed; /* will be needed for random number seed generation */
const gsl_rng_type *T;
gsl_rng *r;
gsl_rng_env_setup();
T = gsl_rng_default;
r = gsl_rng_alloc (T);
seedit("s"); /* Ensure that the next run is going to be a new seed*/
seed = seedrd(); /* specify the seed */
gsl_rng_set(r , seed); /* set the random variable environments
for gsl */
drand48();
seedit("end"); /* Ensure that the next run is going to be a new seed*/
And (my) functions are defined as:
/* ------- related functions -------*/
int seedrd(void)
{
FILE *fopen(), *pfseed;
unsigned short op1;
pfseed = fopen("seedms1", "r"); /* open seedms */
fscanf(pfseed, "%hd", &op1); /* take the first value,
store in address of op1 */
fclose(pfseed); /* close seedms */
return op1; /* return op1 */
}
void seedit( char *flag )
{
FILE *fopen(), *pfseed;
unsigned short seedv[3], seedv2[3], *seed48(), *pseed ;
int i;
if( flag[0] == 's' )
{
pfseed = fopen("seedms1","r");
if( pfseed == NULL )
{
seedv[0] = 3579 ; seedv[1] = 27011;
seedv[2] = 59243;
}
else
{
seedv2[0] = 3579; seedv2[1] = 27011;
seedv2[2] = 59243;
for(i=0;i<3;i++)
{
if( fscanf(pfseed," %hd",seedv+i) < 1 )
seedv[i] = seedv2[i] ;
}
fclose( pfseed);
}
seed48( seedv );
/*printf("\n%d %d %d\n", seedv[0],
seedv[1], seedv[2] ); */
}
else
{
pfseed = fopen("seedms1","w");
pseed = seed48(seedv);
fprintf(pfseed,"%d %d %d\n",pseed[0],
pseed[1],pseed[2]);
}
}
- [Bug-gsl] gsl_max and nested random variables,
K Okamoto <=