summaryrefslogtreecommitdiff
path: root/libgfortran/intrinsics/random.c
diff options
context:
space:
mode:
Diffstat (limited to 'libgfortran/intrinsics/random.c')
-rw-r--r--libgfortran/intrinsics/random.c362
1 files changed, 362 insertions, 0 deletions
diff --git a/libgfortran/intrinsics/random.c b/libgfortran/intrinsics/random.c
new file mode 100644
index 00000000000..b578148f469
--- /dev/null
+++ b/libgfortran/intrinsics/random.c
@@ -0,0 +1,362 @@
+/* Implementation of the RANDOM intrinsics
+ Copyright 2002 Free Software Foundation, Inc.
+ Contributed by Lars Segerlund <seger@linuxmail.org>
+
+ 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 ).
+
+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. */
+
+#include "config.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <assert.h>
+#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. */
+
+#define random_seed prefix(random_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 (get->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];
+ }
+ }
+ }
+}
+