gsl_max and nested random variables

["K Okamoto" <okamotok@...........>]

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]);
                        }
        }