/* -*- Mode: c; c-basic-offset: 2 -*-
 *
 * mt.c - Mersenne Twister functions
 *
 * This is free and unencumbered software released into the public domain.
 *
 * Anyone is free to copy, modify, publish, use, compile, sell, or
 * distribute this software, either in source code form or as a compiled
 * binary, for any purpose, commercial or non-commercial, and by any
 * means.
 *
 * In jurisdictions that recognize copyright laws, the author or authors
 * of this software dedicate any and all copyright interest in the
 * software to the public domain. We make this dedication for the benefit
 * of the public at large and to the detriment of our heirs and
 * successors. We intend this dedication to be an overt act of
 * relinquishment in perpetuity of all present and future rights to this
 * software under copyright law.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 * IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
 * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 * OTHER DEALINGS IN THE SOFTWARE.
 *
 * For more information, please refer to <http://unlicense.org/>
 * 
 */


#ifdef MTWIST_CONFIG
#include <mtwist_config.h>
#endif

#include <stdio.h>

#ifdef HAVE_STDLIB_H
#include <stdlib.h>
#endif

#ifdef HAVE_STDINT_H
/* To get UINT32_C */
#define __STDC_CONSTANT_MACROS 1
#include <stdint.h>
#else
/* pessimistic - use an unsigned long */
typedef unsigned long uint32_t;
#define UINT32_C(v) (v ## UL)
#endif


#include <mtwist.h>

#include <mtwist_internal.h>

/*
 * Mersenne Twister (MT19937) algorithm
 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
 *
 * http://en.wikipedia.org/wiki/Mersenne_twister
 *
 */

#define MTWIST_UPPER_MASK    UINT32_C(0x80000000)
#define MTWIST_LOWER_MASK    UINT32_C(0x7FFFFFFF)
#define MTWIST_FULL_MASK     UINT32_C(0xFFFFFFFF)

#define MTWIST_MATRIX_A      UINT32_C(0x9908B0DF)

#define MTWIST_MIXBITS(u, v) ( ( (u) & MTWIST_UPPER_MASK) | ( (v) & MTWIST_LOWER_MASK) )
#define MTWIST_TWIST(u, v)  ( (MTWIST_MIXBITS(u, v) >> 1) ^ ( (v) & UINT32_C(1) ? MTWIST_MATRIX_A : UINT32_C(0)) )


/**
 * mtwist_new:
 *
 * Construct a Mersenne Twister object
 *
 * Return value: new MT object or NULL on failure
 */
mtwist*
mtwist_new(void)
{
  mtwist* mt;
  
  mt = (mtwist*)calloc(1, sizeof(*mt));
  if(!mt)
    return NULL;
  
  mt->remaining = 0;
  mt->next = NULL;
  mt->seeded = 0;

  return mt;
}


/**
 * mtwist_free:
 * @mt: mt object
 *
 * Destroy a Mersenne Twister object
 */
void
mtwist_free(mtwist* mt) 
{
  if(mt)
    free(mt);
}


/**
 * mtwist_init:
 * @mt: mt object
 * @seed: seed (lower 32 bits used)
 *
 * Initialise a Mersenne Twister with an unsigned 32 bit int seed
 */
void
mtwist_init(mtwist* mt, unsigned long seed)
{
  int i;

  if(!mt)
    return;
  
  mt->state[0] = (uint32_t)(seed & MTWIST_FULL_MASK);
  for(i = 1; i < MTWIST_N; i++) {
    mt->state[i] = (UINT32_C(1812433253) * (mt->state[i - 1] ^ (mt->state[i - 1] >> 30)) + i); 
    mt->state[i] &= MTWIST_FULL_MASK;
  }

  mt->remaining = 0;
  mt->next = NULL;

  mt->seeded = 1;
}


static void
mtwist_update_state(mtwist* mt)
{
  int count;
  uint32_t *p = mt->state;

  for(count = (MTWIST_N - MTWIST_M + 1); --count; p++)
    *p = p[MTWIST_M] ^ MTWIST_TWIST(p[0], p[1]);

  for(count = MTWIST_M; --count; p++)
    *p = p[MTWIST_M - MTWIST_N] ^ MTWIST_TWIST(p[0], p[1]);

  *p = p[MTWIST_M - MTWIST_N] ^ MTWIST_TWIST(p[0], mt->state[0]);

  mt->remaining = MTWIST_N;
  mt->next = mt->state;
}


/**
 * mtwist_u32rand:
 * @mt: mt object
 *
 * Get a random unsigned 32 bit integer from the random number generator
 *
 * Return value: unsigned long with 32 valid bits
 */
unsigned long
mtwist_u32rand(mtwist* mt)
{
  uint32_t r;

  if(!mt)
    return 0UL;

  if(!mt->seeded)
    mtwist_init(mt, mtwist_seed_from_system(mt));

  if(!mt->remaining)
    mtwist_update_state(mt);
  
  r = *mt->next++;
  mt->remaining--;

  /* Tempering */
  r ^= (r >> 11);
  r ^= (r << 7) & UINT32_C(0x9D2C5680);
  r ^= (r << 15) & UINT32_C(0xEFC60000);
  r ^= (r >> 18);

  r &= MTWIST_FULL_MASK;

  return (unsigned long)r;
}


/**
 * mtwist_drand:
 * @mt: mt object
 *
 * Get a random double from the random number generator
 *
 * Return value: random double in the range 0.0 inclusive to 1.0 exclusive; [0.0, 1.0) */
double
mtwist_drand(mtwist* mt)
{
  unsigned long r;
  double d;

  if(!mt)
    return 0.0;

  r = mtwist_u32rand(mt);

  d = r / 4294967296.0; /* 2^32 */

  return d;
}