#include "config.h"
#include "libgfortran.h"
#include <gthr.h>
extern void random_r4 (GFC_REAL_4 *);
iexport_proto(random_r4);
extern void random_r8 (GFC_REAL_8 *);
iexport_proto(random_r8);
extern void arandom_r4 (gfc_array_r4 *);
export_proto(arandom_r4);
extern void arandom_r8 (gfc_array_r8 *);
export_proto(arandom_r8);
#ifdef HAVE_GFC_REAL_10
extern void random_r10 (GFC_REAL_10 *);
iexport_proto(random_r10);
extern void arandom_r10 (gfc_array_r10 *);
export_proto(arandom_r10);
#endif
#ifdef HAVE_GFC_REAL_16
extern void random_r16 (GFC_REAL_16 *);
iexport_proto(random_r16);
extern void arandom_r16 (gfc_array_r16 *);
export_proto(arandom_r16);
#endif
#ifdef __GTHREAD_MUTEX_INIT
static __gthread_mutex_t random_lock = __GTHREAD_MUTEX_INIT;
#else
static __gthread_mutex_t random_lock;
#endif
static inline void
rnumber_4 (GFC_REAL_4 *f, GFC_UINTEGER_4 v)
{
GFC_UINTEGER_4 mask;
#if GFC_REAL_4_RADIX == 2
mask = ~ (GFC_UINTEGER_4) 0u << (32 - GFC_REAL_4_DIGITS);
#elif GFC_REAL_4_RADIX == 16
mask = ~ (GFC_UINTEGER_4) 0u << ((8 - GFC_REAL_4_DIGITS) * 4);
#else
#error "GFC_REAL_4_RADIX has unknown value"
#endif
v = v & mask;
*f = (GFC_REAL_4) v * (GFC_REAL_4) 0x1.p-32f;
}
static inline void
rnumber_8 (GFC_REAL_8 *f, GFC_UINTEGER_8 v)
{
GFC_UINTEGER_8 mask;
#if GFC_REAL_8_RADIX == 2
mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_8_DIGITS);
#elif GFC_REAL_8_RADIX == 16
mask = ~ (GFC_UINTEGER_8) 0u << (16 - GFC_REAL_8_DIGITS) * 4);
#else
#error "GFC_REAL_8_RADIX has unknown value"
#endif
v = v & mask;
*f = (GFC_REAL_8) v * (GFC_REAL_8) 0x1.p-64;
}
#ifdef HAVE_GFC_REAL_10
static inline void
rnumber_10 (GFC_REAL_10 *f, GFC_UINTEGER_8 v)
{
GFC_UINTEGER_8 mask;
#if GFC_REAL_10_RADIX == 2
mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_10_DIGITS);
#elif GFC_REAL_10_RADIX == 16
mask = ~ (GFC_UINTEGER_10) 0u << ((16 - GFC_REAL_10_DIGITS) * 4);
#else
#error "GFC_REAL_10_RADIX has unknown value"
#endif
v = v & mask;
*f = (GFC_REAL_10) v * (GFC_REAL_10) 0x1.p-64;
}
#endif
#ifdef HAVE_GFC_REAL_16
static inline void
rnumber_16 (GFC_REAL_16 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2)
{
GFC_UINTEGER_8 mask;
#if GFC_REAL_16_RADIX == 2
mask = ~ (GFC_UINTEGER_8) 0u << (128 - GFC_REAL_16_DIGITS);
#elif GFC_REAL_16_RADIX == 16
mask = ~ (GFC_UINTEGER_8) 0u << ((32 - GFC_REAL_16_DIGITS) * 4);
#else
#error "GFC_REAL_16_RADIX has unknown value"
#endif
v2 = v2 & mask;
*f = (GFC_REAL_16) v1 * (GFC_REAL_16) 0x1.p-64
+ (GFC_REAL_16) v2 * (GFC_REAL_16) 0x1.p-128;
}
#endif
#define GFC_SL(k, n) ((k)^((k)<<(n)))
#define GFC_SR(k, n) ((k)^((k)>>(n)))
#define KISS_DEFAULT_SEED_1 123456789, 362436069, 521288629, 316191069
#define KISS_DEFAULT_SEED_2 987654321, 458629013, 582859209, 438195021
#ifdef HAVE_GFC_REAL_16
#define KISS_DEFAULT_SEED_3 573658661, 185639104, 582619469, 296736107
#endif
static GFC_UINTEGER_4 kiss_seed[] = {
KISS_DEFAULT_SEED_1,
KISS_DEFAULT_SEED_2,
#ifdef HAVE_GFC_REAL_16
KISS_DEFAULT_SEED_3
#endif
};
static GFC_UINTEGER_4 kiss_default_seed[] = {
KISS_DEFAULT_SEED_1,
KISS_DEFAULT_SEED_2,
#ifdef HAVE_GFC_REAL_16
KISS_DEFAULT_SEED_3
#endif
};
static const GFC_INTEGER_4 kiss_size = sizeof(kiss_seed)/sizeof(kiss_seed[0]);
static GFC_UINTEGER_4 * const kiss_seed_1 = kiss_seed;
static GFC_UINTEGER_4 * const kiss_seed_2 = kiss_seed + 4;
#ifdef HAVE_GFC_REAL_16
static GFC_UINTEGER_4 * const kiss_seed_3 = kiss_seed + 8;
#endif
static GFC_UINTEGER_4
kiss_random_kernel(GFC_UINTEGER_4 * seed)
{
GFC_UINTEGER_4 kiss;
seed[0] = 69069 * seed[0] + 1327217885;
seed[1] = GFC_SL(GFC_SR(GFC_SL(seed[1],13),17),5);
seed[2] = 18000 * (seed[2] & 65535) + (seed[2] >> 16);
seed[3] = 30903 * (seed[3] & 65535) + (seed[3] >> 16);
kiss = seed[0] + seed[1] + (seed[2] << 16) + seed[3];
return kiss;
}
void
random_r4 (GFC_REAL_4 *x)
{
GFC_UINTEGER_4 kiss;
__gthread_mutex_lock (&random_lock);
kiss = kiss_random_kernel (kiss_seed_1);
rnumber_4 (x, kiss);
__gthread_mutex_unlock (&random_lock);
}
iexport(random_r4);
void
random_r8 (GFC_REAL_8 *x)
{
GFC_UINTEGER_8 kiss;
__gthread_mutex_lock (&random_lock);
kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
kiss += kiss_random_kernel (kiss_seed_2);
rnumber_8 (x, kiss);
__gthread_mutex_unlock (&random_lock);
}
iexport(random_r8);
#ifdef HAVE_GFC_REAL_10
void
random_r10 (GFC_REAL_10 *x)
{
GFC_UINTEGER_8 kiss;
__gthread_mutex_lock (&random_lock);
kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
kiss += kiss_random_kernel (kiss_seed_2);
rnumber_10 (x, kiss);
__gthread_mutex_unlock (&random_lock);
}
iexport(random_r10);
#endif
#ifdef HAVE_GFC_REAL_16
void
random_r16 (GFC_REAL_16 *x)
{
GFC_UINTEGER_8 kiss1, kiss2;
__gthread_mutex_lock (&random_lock);
kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
kiss1 += kiss_random_kernel (kiss_seed_2);
kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
kiss2 += kiss_random_kernel (kiss_seed_3);
rnumber_16 (x, kiss1, kiss2);
__gthread_mutex_unlock (&random_lock);
}
iexport(random_r16);
#endif
void
arandom_r4 (gfc_array_r4 *x)
{
index_type count[GFC_MAX_DIMENSIONS];
index_type extent[GFC_MAX_DIMENSIONS];
index_type stride[GFC_MAX_DIMENSIONS];
index_type stride0;
index_type dim;
GFC_REAL_4 *dest;
GFC_UINTEGER_4 kiss;
int n;
dest = x->data;
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];
__gthread_mutex_lock (&random_lock);
while (dest)
{
kiss = kiss_random_kernel (kiss_seed_1);
rnumber_4 (dest, kiss);
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];
}
}
}
__gthread_mutex_unlock (&random_lock);
}
void
arandom_r8 (gfc_array_r8 *x)
{
index_type count[GFC_MAX_DIMENSIONS];
index_type extent[GFC_MAX_DIMENSIONS];
index_type stride[GFC_MAX_DIMENSIONS];
index_type stride0;
index_type dim;
GFC_REAL_8 *dest;
GFC_UINTEGER_8 kiss;
int n;
dest = x->data;
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];
__gthread_mutex_lock (&random_lock);
while (dest)
{
kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
kiss += kiss_random_kernel (kiss_seed_2);
rnumber_8 (dest, kiss);
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];
}
}
}
__gthread_mutex_unlock (&random_lock);
}
#ifdef HAVE_GFC_REAL_10
void
arandom_r10 (gfc_array_r10 *x)
{
index_type count[GFC_MAX_DIMENSIONS];
index_type extent[GFC_MAX_DIMENSIONS];
index_type stride[GFC_MAX_DIMENSIONS];
index_type stride0;
index_type dim;
GFC_REAL_10 *dest;
GFC_UINTEGER_8 kiss;
int n;
dest = x->data;
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];
__gthread_mutex_lock (&random_lock);
while (dest)
{
kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
kiss += kiss_random_kernel (kiss_seed_2);
rnumber_10 (dest, kiss);
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];
}
}
}
__gthread_mutex_unlock (&random_lock);
}
#endif
#ifdef HAVE_GFC_REAL_16
void
arandom_r16 (gfc_array_r16 *x)
{
index_type count[GFC_MAX_DIMENSIONS];
index_type extent[GFC_MAX_DIMENSIONS];
index_type stride[GFC_MAX_DIMENSIONS];
index_type stride0;
index_type dim;
GFC_REAL_16 *dest;
GFC_UINTEGER_8 kiss1, kiss2;
int n;
dest = x->data;
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];
__gthread_mutex_lock (&random_lock);
while (dest)
{
kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
kiss1 += kiss_random_kernel (kiss_seed_2);
kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
kiss2 += kiss_random_kernel (kiss_seed_3);
rnumber_16 (dest, kiss1, kiss2);
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];
}
}
}
__gthread_mutex_unlock (&random_lock);
}
#endif
void
random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
{
int i;
__gthread_mutex_lock (&random_lock);
if (size == NULL && put == NULL && get == NULL)
{
for (i=0; i<kiss_size; i++)
kiss_seed[i] = kiss_default_seed[i];
}
if (size != NULL)
*size = kiss_size;
if (put != NULL)
{
if (GFC_DESCRIPTOR_RANK (put) != 1)
runtime_error ("Array rank of PUT is not 1.");
if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
runtime_error ("Array size of PUT is too small.");
for (i = 0; i < kiss_size; i++)
kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
}
if (get != NULL)
{
if (GFC_DESCRIPTOR_RANK (get) != 1)
runtime_error ("Array rank of GET is not 1.");
if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
runtime_error ("Array size of GET is too small.");
for (i = 0; i < kiss_size; i++)
get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
}
__gthread_mutex_unlock (&random_lock);
}
iexport(random_seed);
#ifndef __GTHREAD_MUTEX_INIT
static void __attribute__((constructor))
init (void)
{
__GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
}
#endif