/*
* lcg.c - Implementation of a linear congruential generator with
* log(n) skipping. For documentation, see the header file
* lcg.h.
*
* (C) 2016 Ramses van Zon, SciNet HPC Consortium, Toronto, Canada
*/
#include "lcg.h"
#include <assert.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
const double pi=3.1415926535897932; /* value M_PI not present in c99/c11 */
static void lcg_advance(lcg_t* lcg)
{
/* Perform one step in the generator. Not exported, as
* lcg_skip(lcg,1) does the same */
assert(lcg != NULL);
lcg->s = (lcg->s * lcg->a + lcg->b) % lcg->m;
}
static uint64_t rand2lcgseed(int32_t seed)
{
/* To initialize the internal buffer r(n) of drand48, lrand48, and
* mrand48, the 32 bits of the seed value are copied into the
* upper 32 bits of its state buffer (a "short int r[n]" with
* n=3), with the lower 16 bits of r[n] arbitrarily being set to
* 0x330e. */
uint64_t seedll = 0;
memcpy(&seedll, &seed, sizeof(seed));
return (seedll<<16) + 0x330e;
}
static int32_t u64_2_s32(uint64_t x, int n)
{
/* The algorithm that mrand48 uses to extract a signed 32-bit
* value from the 48-bit unsigned long value is to do a 16 bits
* right-shift. The extra mask is to ensure there is not overflow
* for the 32-bit long value. The top bit becomes the sign, so
* the result may be negative. Here generalized for n bits
* shifted, rand48 uses n=16. */
assert(n >= 0);
return (int32_t) ((x>>n)&0xffffffff);
}
static int32_t u64_2_s32_positive(uint64_t x, int n)
{
/* The algorithm that lrand48 uses to extract a signed long value
* that is positive from the 48-bit unsigned long value is to do a
* 16+1 bits right-shift. The +1 in 16+1 is to clear
* the sign bit, which is further ensured by the 0x7fffffff
* mask. Here generalized for n+1 bits shifted, rand48 uses
* n=16. */
assert(n >= 0);
return (long) ((x>>(n+1)) & 0x7fffffff);
}
static unsigned long u64_2_u32(uint64_t x, int n)
{
/* An algorithm similar to lrand48 to extract an unsigned long
* value from the 48-bit unsigned long value is to do a 16 bits
* right-shift. The extra mask is to ensure there is not overflow
* for the 32 bit long value. Here generalized for n bits
* shifted. */
assert(n >= 0);
return (unsigned long) ((x>>n) & 0xffffffff);
}
static int hibit(uint64_t m)
{
unsigned r = 0;
while (m >>= 1)
r++;
return r;
}
lcg_t lcg_init(uint64_t a, uint64_t b, uint64_t m, uint64_t seed)
{
assert(a > 0);
assert(b >= 0);
assert(m > 0);
assert(0 <= seed && seed < m);
return (lcg_t){a,b,m,seed};
}
lcg_t lcg_init_alt35(int64_t seed)
{
return (lcg_t){
.a = 34359738368LL,
.b = 1LL,
.m = 1LL << 35,
.s = seed
};
}
lcg_t lcg_init_rand48(int32_t seed)
{
uint64_t seed64 = rand2lcgseed(seed);
return (lcg_t){
.a = 25214903917LL,
.b = 11LL,
.m = 1LL << 48,
.s = seed64
};
}
void lcg_seed(lcg_t* lcg, uint64_t seed)
{
assert(lcg->a > 0);
assert(lcg->b >= 0);
assert(lcg->m > 0);
assert(0 <= seed && seed < lcg->m);
lcg->s = seed;
}
int32_t lcg_int32(lcg_t* lcg)
{
assert(lcg != NULL);
lcg_advance(lcg);
int n = hibit(lcg->m)-32;
return u64_2_s32(lcg->s, n);
}
uint32_t lcg_uint32(lcg_t* lcg)
{
assert(lcg != NULL);
lcg_advance(lcg);
int n = hibit(lcg->m)-32;
return u64_2_u32(lcg->s, n);
}
int32_t lcg_pint32(lcg_t* lcg)
{
assert(lcg != NULL);
lcg_advance(lcg);
int n = hibit(lcg->m)-32;
return u64_2_s32_positive(lcg->s, n);
}
uint64_t lcg_int64(lcg_t* lcg)
{
assert(lcg != NULL);
lcg_advance(lcg);
return lcg->s;
}
float lcg_float(lcg_t* lcg)
{
assert(lcg != NULL);
lcg_advance(lcg);
return (float)(lcg->s)/(float)(lcg->m);
}
double lcg_double(lcg_t* lcg)
{
assert(lcg != NULL);
lcg_advance(lcg);
return (double)(lcg->s)/(double)(lcg->m);
}
/* Inline macro for LCG multiplication, i.e. LGC1 *= LGC2.
*
* Based on the formula:
*
* LCG(a1,b1,m) * LCG(a2,b2,m) = LCG((a1*a2)%m, (b1+a1*b2)%m)
*
* Note: A macro rather than a function because we want the result
* in-place. It uses the usual do/while trick for avoiding macro
* side-effects.
*/
#define LCGMULRIGHT(a1,b1,a2,b2,m) \
do { \
const uint64_t tmp = a1; \
a1 = (a1*a2)%m; \
b1 = (b1+tmp*b2)%m; \
} while(0)
void lcg_skip(lcg_t* lcg, int64_t n)
{
/* O(log(n)) skip ahead function for LCG. The algorithm is a
* variant of modular exponentiation.
* (https://en.wikipedia.org/wiki/Modular_nonentiation), but with
* multiplication replaced by multiplication of LCGs (macro
* above). If n<0, the function will abort through an assert, or,
* if compile with NDEBUG defined, will not skip at all. */
assert(lcg != NULL);
assert(n >= 0);
uint64_t a = 1;
uint64_t b = 0;
uint64_t basea = lcg->a % lcg->m;
uint64_t baseb = lcg->b % lcg->m;
while (n > 0) {
if ((n % 2) == 1)
LCGMULRIGHT(a,b,basea,baseb,lcg->m);
n >>= 1;
LCGMULRIGHT(basea,baseb,basea,baseb,lcg->m);
}
lcg->s = (lcg->s * a + b) % lcg->m;
}
/* rand48 replacements */
void lcg_srand48(lcg_t* lcg, int32_t seed)
{
/* srand48 is used to initialize the internal buffer r(n) of
* drand48, lrand48, and mrand48 such that the 32 bits of the seed
* value are copied into the upper 32 bits of r(n), with the lower
* 16 bits of r(n) arbitrarily being set to 0x330e. */
lcg_seed(lcg, rand2lcgseed(seed));
}
long lcg_lrand48(lcg_t* lcg)
{
/* lrand48 and nrand48 return values of type long in the range [0,
* 2**31-1]. The high-order (31) bits of r(n+1) are loaded into
* the lower bits of the returned value, with the topmost (sign)
* bit set to zero. */
return u64_2_s32_positive(lcg_int64(lcg),16);
}
long lcg_mrand48(lcg_t* lcg)
{
/* mrand48 and jrand48 return values of type long in the range
* [-2**31, 2**31-1]. The high-order (32) bits of r(n+1) are
* loaded into the returned value. */
return u64_2_s32(lcg_int64(lcg),16);
}
double lcg_drand48(lcg_t* lcg)
{
/* drand48 and erand48 return values of type double. The full 48
* bits of r(n+1) are loaded into the mantissa of the returned
* value, with the exponent set such that the values produced lie
* in the interval [0.0, 1.0]. */
union {
uint64_t long_long_value;
uint16_t x[4];
} bitfill;
bitfill.long_long_value = lcg_int64(lcg);
const double z = 1.0 / (1L << 16);
return (z * (z * (z * bitfill.x[0] + bitfill.x[1]) + bitfill.x[2]));
}
double lcg_normal(lcg_t* lcg)
{
double fac, y1, y2, x1;
if (lcg->have == 1) { /* already one available ? */
lcg->have = 0;
return lcg->x2;
} else {
/* generate a pair of random variables */
y1 = lcg_drand48(lcg);
y2 = lcg_drand48(lcg);
fac = sqrt(-2*log(y1));
lcg->have = 1;
x1 = fac*sin(2*pi*y2); /* x1 and x2 are now gaussian */
lcg->x2 = fac*cos(2*pi*y2); /* so store one */
return x1; /* and return the other. */
}
}
