gcc/libgfortran/intrinsics/random.c

/* Implementation of the RANDOM intrinsics
   Copyright 2002, 2004 Free Software Foundation, Inc.
   Contributed by Lars Segerlund <seger@linuxmail.org>
   and Steve Kargl.

This file is part of the GNU Fortran 95 runtime library (libgfortran).

Libgfortran is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.

Ligbfortran is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with libgfor; see the file COPYING.LIB.  If not,
write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA.  */

#if 0

/*  The Mersenne Twister code is currently commented out due to

    (1) Simple user specified seeds lead to really bad sequences for
        nearly 100000 random numbers.
    (2) open(), read(), and close() are not properly declared via header
        files.
    (3) The global index i is abused and causes unexpected behavior with
        GET and PUT.
    (4) See PR 15619.

  The algorithm was taken from the paper :

	Mersenne Twister:	623-dimensionally equidistributed
				uniform pseudorandom generator.

	by:	Makoto Matsumoto
		Takuji Nishimura

	Which appeared in the: ACM Transactions on Modelling and Computer
	Simulations: Special Issue on Uniform Random Number
	Generation. ( Early in 1998 ).  */


#include "config.h"
#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>

#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif

#include "libgfortran.h"

/*Use the 'big' generator by default ( period -> 2**19937 ).  */

#define MT19937

/* Define the necessary constants for the algorithm.  */

#ifdef  MT19937
enum constants
{
  N = 624, M = 397, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
};
#define M_A	0x9908B0DF
#define T_B	0x9D2C5680
#define T_C	0xEFC60000
#else
enum constants
{
  N = 351, M = 175, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
};
#define M_A	0xE4BD75F5
#define T_B	0x655E5280
#define T_C	0xFFD58000
#endif

static int i = N;
static unsigned int seed[N];

/* This is the routine which handles the seeding of the generator,
   and also reading and writing of the seed.  */

void
random_seed (GFC_INTEGER_4 * size, const gfc_array_i4 * put,
	     const gfc_array_i4 * get)
{
  /* Initialize the seed in system dependent manner.  */
  if (get == NULL && put == NULL && size == NULL)
    {
      int fd;
      fd = open ("/dev/urandom", O_RDONLY);
      if (fd == 0)
	{
	  /* We dont have urandom.  */
	  GFC_UINTEGER_4 s = (GFC_UINTEGER_4) seed;
	  for (i = 0; i < N; i++)
	    {
	      s = s * 29943829 - 1;
	      seed[i] = s;
	    }
	}
      else
	{
	  /* Using urandom, might have a length issue.  */
	  read (fd, &seed[0], sizeof (GFC_UINTEGER_4) * N);
	  close (fd);
	}
      return;
    }

  /* Return the size of the seed */
  if (size != NULL)
    {
      *size = N;
      return;
    }

  /* if we have gotten to this pount we have a get or put
   * now we check it the array fulfills the demands in the standard .
   */

  /* Set the seed to PUT data */
  if (put != NULL)
    {
      /* if the rank of the array is not 1 abort */
      if (GFC_DESCRIPTOR_RANK (put) != 1)
	abort ();

      /* if the array is too small abort */
      if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < N)
	abort ();

      /* If this is the case the array is a temporary */
      if (put->dim[0].stride == 0)
	return;

      /*  This code now should do correct strides. */
      for (i = 0; i < N; i++)
	seed[i] = put->data[i * put->dim[0].stride];
    }

  /* Return the seed to GET data */
  if (get != NULL)
    {
      /* if the rank of the array is not 1 abort */
      if (GFC_DESCRIPTOR_RANK (get) != 1)
	abort ();

      /* if the array is too small abort */
      if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < N)
	abort ();

      /* If this is the case the array is a temporary */
      if (get->dim[0].stride == 0)
	return;

      /*  This code now should do correct strides. */
      for (i = 0; i < N; i++)
	get->data[i * get->dim[0].stride] = seed[i];
    }
}

/* Here is the internal routine which generates the random numbers
   in 'batches' based upon the need for a new batch.
   It's an integer based routine known as 'Mersenne Twister'.
   This implementation still lacks 'tempering' and a good verification,
   but gives very good metrics.  */

static void
random_generate (void)
{
  /* 32 bits.  */
  GFC_UINTEGER_4 y;

  /* Generate batch of N.  */
  int k, m;
  for (k = 0, m = M; k < N - 1; k++)
    {
      y = (seed[k] & (-1 << R)) | (seed[k + 1] & ((1u << R) - 1));
      seed[k] = seed[m] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
      if (++m >= N)
	m = 0;
    }

  y = (seed[N - 1] & (-1 << R)) | (seed[0] & ((1u << R) - 1));
  seed[N - 1] = seed[M - 1] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
  i = 0;
}

/* A routine to return a REAL(KIND=4).  */

#define random_r4 prefix(random_r4)
void
random_r4 (GFC_REAL_4 * harv)
{
  /* Regenerate if we need to.  */
  if (i >= N)
    random_generate ();

  /* Convert uint32 to REAL(KIND=4).  */
  *harv = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
			(GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
}

/* A routine to return a REAL(KIND=8).  */

#define random_r8 prefix(random_r8)
void
random_r8 (GFC_REAL_8 * harv)
{
  /* Regenerate if we need to, may waste one 32-bit value.  */
  if ((i + 1) >= N)
    random_generate ();

  /* Convert two uint32 to a REAL(KIND=8).  */
  *harv = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
	  (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
  i += 2;
}

/* Code to handle arrays will follow here.  */

/* REAL(KIND=4) REAL array.  */

#define arandom_r4 prefix(arandom_r4)
void
arandom_r4 (gfc_array_r4 * harv)
{
  index_type count[GFC_MAX_DIMENSIONS - 1];
  index_type extent[GFC_MAX_DIMENSIONS - 1];
  index_type stride[GFC_MAX_DIMENSIONS - 1];
  index_type stride0;
  index_type dim;
  GFC_REAL_4 *dest;
  int n;

  dest = harv->data;

  if (harv->dim[0].stride == 0)
    harv->dim[0].stride = 1;

  dim = GFC_DESCRIPTOR_RANK (harv);

  for (n = 0; n < dim; n++)
    {
      count[n] = 0;
      stride[n] = harv->dim[n].stride;
      extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
      if (extent[n] <= 0)
	return;
    }

  stride0 = stride[0];

  while (dest)
    {
      /* Set the elements.  */

      /* regenerate if we need to */
      if (i >= N)
	random_generate ();

      /* Convert uint32 to float in a hopefully g95 compiant manner */
      *dest = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
			    (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));

      /* Advance to the next element.  */
      dest += stride0;
      count[0]++;
      /* Advance to the next source element.  */
      n = 0;
      while (count[n] == extent[n])
	{
	  /* When we get to the end of a dimension,
	     reset it and increment
	     the next dimension.  */
	  count[n] = 0;
	  /* We could precalculate these products,
	     but this is a less
	     frequently used path so proabably not worth it.  */
	  dest -= stride[n] * extent[n];
	  n++;
	  if (n == dim)
	    {
	      dest = NULL;
	      break;
	    }
	  else
	    {
	      count[n]++;
	      dest += stride[n];
	    }
	}
    }
}

/* REAL(KIND=8) array.  */

#define arandom_r8 prefix(arandom_r8)
void
arandom_r8 (gfc_array_r8 * harv)
{
  index_type count[GFC_MAX_DIMENSIONS - 1];
  index_type extent[GFC_MAX_DIMENSIONS - 1];
  index_type stride[GFC_MAX_DIMENSIONS - 1];
  index_type stride0;
  index_type dim;
  GFC_REAL_8 *dest;
  int n;

  dest = harv->data;

  if (harv->dim[0].stride == 0)
    harv->dim[0].stride = 1;

  dim = GFC_DESCRIPTOR_RANK (harv);

  for (n = 0; n < dim; n++)
    {
      count[n] = 0;
      stride[n] = harv->dim[n].stride;
      extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
      if (extent[n] <= 0)
	return;
    }

  stride0 = stride[0];

  while (dest)
    {
      /* Set the elements.  */

      /* regenerate if we need to, may waste one 32-bit value */
      if ((i + 1) >= N)
	random_generate ();

      /* Convert two uint32 to a REAL(KIND=8).  */
      *dest = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
	      (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
      i += 2;

      /* Advance to the next element.  */
      dest += stride0;
      count[0]++;
      /* Advance to the next source element.  */
      n = 0;
      while (count[n] == extent[n])
	{
	  /* When we get to the end of a dimension,
	     reset it and increment
	     the next dimension.  */
	  count[n] = 0;
	  /* We could precalculate these products,
	     but this is a less
	     frequently used path so proabably not worth it.  */
	  dest -= stride[n] * extent[n];
	  n++;
	  if (n == dim)
	    {
	      dest = NULL;
	      break;
	    }
	  else
	    {
	      count[n]++;
	      dest += stride[n];
	    }
	}
    }
}
#endif /* Mersenne Twister code */


/* George Marsaglia's KISS (Keep It Simple Stupid) random number generator.

   This PRNG combines:

   (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
       of 2^32,
   (2) A 3-shift shift-register generator with a period of 2^32-1,
   (3) Two 16-bit multiply-with-carry generators with a period of
       597273182964842497 > 2^59.

   The overall period exceeds 2^123.

   http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu

   The above web site has an archive of a newsgroup posting from George
   Marsaglia with the statement:

   Subject: Random numbers for C: Improvements.
   Date: Fri, 15 Jan 1999 11:41:47 -0500
   From: George Marsaglia <geo@stat.fsu.edu>
   Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
   References: <369B5E30.65A55FD1@stat.fsu.edu>
   Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
   Lines: 93

   As I hoped, several suggestions have led to
   improvements in the code for RNG's I proposed for
   use in C. (See the thread "Random numbers for C: Some
   suggestions" in previous postings.) The improved code
   is listed below.

   A question of copyright has also been raised.  Unlike
   DIEHARD, there is no copyright on the code below. You
   are free to use it in any way you want, but you may
   wish to acknowledge the source, as a courtesy.

"There is no copyright on the code below." included the original
KISS algorithm.  */

#include "config.h"
#include "libgfortran.h"

#define GFC_SL(k, n)	((k)^((k)<<(n)))
#define GFC_SR(k, n)	((k)^((k)>>(n)))

static const GFC_INTEGER_4 kiss_size = 4;
#define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069};
static const GFC_UINTEGER_4 kiss_default_seed[4] = KISS_DEFAULT_SEED;
static GFC_UINTEGER_4 kiss_seed[4] = KISS_DEFAULT_SEED;

/* kiss_random_kernel() returns an integer value in the range of
   (0, GFC_UINTEGER_4_HUGE].  The distribution of pseudorandom numbers
   should be uniform.  */

static GFC_UINTEGER_4
kiss_random_kernel(void)
{

  GFC_UINTEGER_4 kiss;

  kiss_seed[0] = 69069 * kiss_seed[0] + 1327217885;
  kiss_seed[1] = GFC_SL(GFC_SR(GFC_SL(kiss_seed[1],13),17),5);
  kiss_seed[2] = 18000 * (kiss_seed[2] & 65535) + (kiss_seed[2] >> 16);
  kiss_seed[3] = 30903 * (kiss_seed[3] & 65535) + (kiss_seed[3] >> 16);
  kiss = kiss_seed[0] + kiss_seed[1] + (kiss_seed[2] << 16) + kiss_seed[3];

  return kiss;

}

/*  This function produces a REAL(4) value from the uniform distribution
    with range [0,1).  */

void
prefix(random_r4) (GFC_REAL_4 *x)
{

  GFC_UINTEGER_4 kiss;

  kiss = kiss_random_kernel ();
  /* Burn a random number, so the REAL*4 and REAL*8 functions
     produce similar sequences of random numbers.  */
  kiss_random_kernel ();
  *x = normalize_r4_i4 (kiss, ~(GFC_UINTEGER_4) 0);
}

/*  This function produces a REAL(8) value from the uniform distribution
    with range [0,1).  */

void
prefix(random_r8) (GFC_REAL_8 *x)
{

  GFC_UINTEGER_8 kiss;

  kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32;
  kiss += kiss_random_kernel ();
  *x = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0);
}

/*  This function fills a REAL(4) array with values from the uniform
    distribution with range [0,1).  */

void
prefix(arandom_r4) (gfc_array_r4 *x)
{

  index_type count[GFC_MAX_DIMENSIONS - 1];
  index_type extent[GFC_MAX_DIMENSIONS - 1];
  index_type stride[GFC_MAX_DIMENSIONS - 1];
  index_type stride0;
  index_type dim;
  GFC_REAL_4 *dest;
  int n;

  dest = x->data;

  if (x->dim[0].stride == 0)
    x->dim[0].stride = 1;

  dim = GFC_DESCRIPTOR_RANK (x);

  for (n = 0; n < dim; n++)
    {
      count[n] = 0;
      stride[n] = x->dim[n].stride;
      extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
      if (extent[n] <= 0)
        return;
    }

  stride0 = stride[0];

  while (dest)
    {
      prefix(random_r4) (dest);

      /* Advance to the next element.  */
      dest += stride0;
      count[0]++;
      /* Advance to the next source element.  */
      n = 0;
      while (count[n] == extent[n])
        {
          /* When we get to the end of a dimension, reset it and increment
             the next dimension.  */
          count[n] = 0;
          /* We could precalculate these products, but this is a less
             frequently used path so probably not worth it.  */
          dest -= stride[n] * extent[n];
          n++;
          if (n == dim)
            {
              dest = NULL;
              break;
            }
          else
            {
              count[n]++;
              dest += stride[n];
            }
        }
    }
}

/*  This function fills a REAL(8) array with values from the uniform
    distribution with range [0,1).  */

void
prefix(arandom_r8) (gfc_array_r8 *x)
{

  index_type count[GFC_MAX_DIMENSIONS - 1];
  index_type extent[GFC_MAX_DIMENSIONS - 1];
  index_type stride[GFC_MAX_DIMENSIONS - 1];
  index_type stride0;
  index_type dim;
  GFC_REAL_8 *dest;
  int n;

  dest = x->data;

  if (x->dim[0].stride == 0)
    x->dim[0].stride = 1;

  dim = GFC_DESCRIPTOR_RANK (x);

  for (n = 0; n < dim; n++)
    {
      count[n] = 0;
      stride[n] = x->dim[n].stride;
      extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
      if (extent[n] <= 0)
        return;
    }

  stride0 = stride[0];

  while (dest)
    {
      prefix(random_r8) (dest);

      /* Advance to the next element.  */
      dest += stride0;
      count[0]++;
      /* Advance to the next source element.  */
      n = 0;
      while (count[n] == extent[n])
        {
          /* When we get to the end of a dimension, reset it and increment
             the next dimension.  */
          count[n] = 0;
          /* We could precalculate these products, but this is a less
             frequently used path so probably not worth it.  */
          dest -= stride[n] * extent[n];
          n++;
          if (n == dim)
            {
              dest = NULL;
              break;
            }
          else
            {
              count[n]++;
              dest += stride[n];
            }
        }
    }
}

/* prefix(random_seed) is used to seed the PRNG with either a default
   set of seeds or user specified set of seeds.  prefix(random_seed) 
   must be called with no argument or exactly one argument.  */

void
random_seed (GFC_INTEGER_4 *size, gfc_array_i4 * put, 
		     gfc_array_i4 * get)
{

  int i;

  if (size == NULL && put == NULL && get == NULL)
    {
      /* From the standard: "If no argument is present, the processor assigns
         a processor-dependent value to the seed."  */
      kiss_seed[0] = kiss_default_seed[0];
      kiss_seed[1] = kiss_default_seed[1];
      kiss_seed[2] = kiss_default_seed[2];
      kiss_seed[3] = kiss_default_seed[3];
    }

  if (size != NULL)
    *size = kiss_size;

  if (put != NULL)
    {
      /* If the rank of the array is not 1, abort.  */
      if (GFC_DESCRIPTOR_RANK (put) != 1)
        runtime_error ("Array rank of PUT is not 1.");

      /* If the array is too small, abort.  */
      if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
        runtime_error ("Array size of PUT is too small.");

      if (put->dim[0].stride == 0)
	put->dim[0].stride = 1;

      /*  This code now should do correct strides.  */
      for (i = 0; i < kiss_size; i++)
        kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
    }

  /* Return the seed to GET data.  */
  if (get != NULL)
    {
      /* If the rank of the array is not 1, abort.  */
      if (GFC_DESCRIPTOR_RANK (get) != 1)
	runtime_error ("Array rank of GET is not 1.");

      /* If the array is too small, abort.  */
      if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
	runtime_error ("Array size of GET is too small.");

      if (get->dim[0].stride == 0)
	get->dim[0].stride = 1;

      /*  This code now should do correct strides.  */
      for (i = 0; i < kiss_size; i++)
        get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
    }
}