#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"
#define MT19937
#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];
void
random_seed (GFC_INTEGER_4 * size, const gfc_array_i4 * put,
const gfc_array_i4 * get)
{
if (get == NULL && put == NULL && size == NULL)
{
int fd;
fd = open ("/dev/urandom", O_RDONLY);
if (fd == 0)
{
GFC_UINTEGER_4 s = (GFC_UINTEGER_4) seed;
for (i = 0; i < N; i++)
{
s = s * 29943829 - 1;
seed[i] = s;
}
}
else
{
read (fd, &seed[0], sizeof (GFC_UINTEGER_4) * N);
close (fd);
}
return;
}
if (size != NULL)
{
*size = N;
return;
}
if (put != NULL)
{
if (GFC_DESCRIPTOR_RANK (put) != 1)
abort ();
if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < N)
abort ();
if (put->dim[0].stride == 0)
return;
for (i = 0; i < N; i++)
seed[i] = put->data[i * put->dim[0].stride];
}
if (get != NULL)
{
if (GFC_DESCRIPTOR_RANK (get) != 1)
abort ();
if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < N)
abort ();
if (get->dim[0].stride == 0)
return;
for (i = 0; i < N; i++)
get->data[i * get->dim[0].stride] = seed[i];
}
}
static void
random_generate (void)
{
GFC_UINTEGER_4 y;
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;
}
#define random_r4 prefix(random_r4)
void
random_r4 (GFC_REAL_4 * harv)
{
if (i >= N)
random_generate ();
*harv = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
(GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
}
#define random_r8 prefix(random_r8)
void
random_r8 (GFC_REAL_8 * harv)
{
if ((i + 1) >= N)
random_generate ();
*harv = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
(GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
i += 2;
}
#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)
{
if (i >= N)
random_generate ();
*dest = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
(GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
dest += stride0;
count[0]++;
n = 0;
while (count[n] == extent[n])
{
count[n] = 0;
dest -= stride[n] * extent[n];
n++;
if (n == dim)
{
dest = NULL;
break;
}
else
{
count[n]++;
dest += stride[n];
}
}
}
}
#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)
{
if ((i + 1) >= N)
random_generate ();
*dest = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
(GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
i += 2;
dest += stride0;
count[0]++;
n = 0;
while (count[n] == extent[n])
{
count[n] = 0;
dest -= stride[n] * extent[n];
n++;
if (n == dim)
{
dest = NULL;
break;
}
else
{
count[n]++;
dest += stride[n];
}
}
}
}