/* Hoang-Trong Minh Tuan (Oct, 14, 2010):
   + move the code from .h to .cpp
   + modified the code to add three initialize() functions; making it more easier to separate the allocation
   from the initialization
   
*/
#include <cuda.h>
#include "cutil_inline.h"

// includes, project
#include "saruprngCUDA.h"

/*
 * Copyright (c) 2008 Steve Worley < m a t h g e e k@(my last name).com >
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */

/*
  C++ Saru PRNG by Steve Worley.   This is version 0.91, November 7 2008.

  Saru is a 32 bit PRNG with period of 3666320093*2^32. It passes even
  the most stringent randomness test batteries, including DIEHARD and
  the much more comprehensive BigCrush battery. It is designed for
  evaluation efficiency, a (relatively) compact state to allow fast
  copying, and most notably the ability to very efficiently advance
  (or rewind) its state. This advancement feature is useful for game
  and graphics tools, and is especially attractive for multithreaded
  and GPU applications.
  
  The Saru generator is an original work and (deliberately)
  unpatented. See the algorithm description at
  worley.com/mathgeek/saru.html. Updates/fixes to this code (and a plain C
  implementation) can be found there too.

  -----------------

  Usage:

  Constructors for 0, 1, 2, or 3 integer seeds.
  Saru z, s(12345), t(123, 456), u(123, 456, 789);

  Advance state by 1, and output a 32 bit integer pseudo-random value.
  cout << s.u32() << endl; //  Passes BigCrush and DIEHARD

  Advance state by 1, and output a double precision [0..1) floating point 
  cout << s.d() << endl;

  Advance state by 1, and output a single precision [0..1) floating point 
  cout << s.f() << endl;

  Move the generator state forwards a variable number of steps
  s.advance(steps);

  Efficient state advancement or rewind when delta is known at compiletime
  s.advance<123>(); 
  s.rewind<123>(); 

  Advance state by a different step (other than 1) and output a pseudorandom value.
  cout << s.u32<4>() << endl; // skips forward 4 values and outputs the prand.

  Small structure size (64 bits of plain old data) means it's easy and fast
  to copy, store, pass, or revert the state.
  z=s;

  Fork the PRNG, creating a new independent stream, seeded using
  current generator's state. Template seeding allows multiple
  independent children forks. 
  Saru n=s.fork<123>();


  In practice you will likely extend or wrap this class to provide
  more specific functionality such as different output distributions.

  ---------------------

*/

// partial specialization to make a special case for step of 1
template <> __device__ inline void  Saru::advanceWeyl<1>() /* especially efficient single step  */
{ state.y=state.y+oWeylOffset+((((signed int)state.y)>>31)&oWeylPeriod); }


/* This seeding was carefully tested for good churning with 1, 2, and
   3 bit flips.  All 32 incrementing counters (each of the circular
   shifts) pass the TestU01 Crush tests. */
Saru::Saru(unsigned int seed)
{
    this->initialize(seed);
}
/* unsafe implementation (initialize should be private) -but it's good for this case */
__device__ void Saru::initialize(unsigned int seed){
    state.x  = 0x79dedea3*(seed^(((signed int)seed)>>14));
    state.y = seed ^ (((signed int)state.x)>>8);
    state.x  = state.x + (state.y*(state.y^0xdddf97f5));
    state.y = 0xABCB96F7 + (state.y>>1);   
}
/* seeding from 2 samples. We lose one bit of entropy since our input
   seeds have 64 bits but at the end, after mixing, we have just 63. */
Saru::Saru(unsigned int seed1, unsigned int seed2)
{
    this->initialize(seed1, seed2);
}
__device__ void Saru::initialize(unsigned int seed1, unsigned int seed2) {
    seed2+=seed1<<16;
    seed1+=seed2<<11;
    seed2+=((signed int)seed1)>>7;
    seed1^=((signed int)seed2)>>3;
    seed2*=0xA5366B4D;
    seed2^=seed2>>10;
    seed2^=((signed int)seed2)>>19;
    seed1+=seed2^0x6d2d4e11;
  
    state.x  = 0x79dedea3*(seed1^(((signed int)seed1)>>14));
    state.y = (state.x+seed2) ^ (((signed int)state.x)>>8);
    state.x  = state.x + (state.y*(state.y^0xdddf97f5));
    state.y = 0xABCB96F7 + (state.y>>1);    
}

/* 3 seeds. We have to premix the seeds before dropping to 64 bits.
   TODO: this may be better optimized in a future version */
Saru::Saru(unsigned int seed1, unsigned int seed2, unsigned int seed3) {
    this->initialize(seed1, seed2, seed3);
}
__device__ void Saru::initialize(unsigned int seed1, unsigned int seed2, unsigned int seed3) 
{
    seed3^=(seed1<<7)^(seed2>>6);
    seed2+=(seed1>>4)^(seed3>>15);
    seed1^=(seed2<<9)+(seed3<<8);
    seed3^=0xA5366B4D*((seed2>>11) ^ (seed1<<1));
    seed2+=0x72BE1579*((seed1<<4)  ^ (seed3>>16));
    seed1^=0X3F38A6ED*((seed3>>5)  ^ (((signed int)seed2)>>22));
    seed2+=seed1*seed3;
    seed1+=seed3 ^ (seed2>>2);
    seed2^=((signed int)seed2)>>17;
  
    state.x  = 0x79dedea3*(seed1^(((signed int)seed1)>>14));
    state.y = (state.x+seed2) ^ (((signed int)state.x)>>8);
    state.x  = state.x + (state.y*(state.y^0xdddf97f5));
    state.y = 0xABCB96F7 + (state.y>>1);
}

template <unsigned int offset, unsigned int delta, 
	  unsigned int modulus, unsigned int steps> 
__device__ inline unsigned int Saru::advanceAnyWeyl(unsigned int x) 
{
    const unsigned int fullDelta=CTmultmod<delta, steps%modulus, modulus>::value;
    /* the runtime code boils down to this single constant-filled line. */
    return x+((x-offset>modulus-fullDelta) ? fullDelta-modulus : fullDelta);
}

__device__ inline void Saru::advance(unsigned int steps)
{
    // Computes the LCG advancement AND the Weyl D*E mod m simultaneously

    unsigned int currentA=LCGA;
    unsigned int currentC=LCGC;

    unsigned int currentDelta=oWeylDelta;
    unsigned int netDelta=0;
  
    while (steps) {    
        if (steps&1) {
            state.x=currentA*state.x+currentC; // LCG step
            if (netDelta<oWeylPeriod-currentDelta) netDelta+=currentDelta; 
            else netDelta+=currentDelta-oWeylPeriod;
        }
    
        // Change the LCG to step at twice the rate as before
        currentC+=currentA*currentC;
        currentA*=currentA;          
    
        // Change the Weyl delta to step at 2X rate
        if (currentDelta<oWeylPeriod-currentDelta) currentDelta+=currentDelta;
        else currentDelta+=currentDelta-oWeylPeriod;
    
        steps/=2;
    }
  
    // Apply the net delta to the Weyl state.
    if (state.y-oWeylOffset<oWeylPeriod-netDelta) state.y+=netDelta; 
    else state.y+=netDelta-oWeylPeriod;
}


template <unsigned int seed> 
Saru Saru::fork() const
{
    const unsigned int churned1=0xDEADBEEF ^ (0x1fc4ce47*(seed^(seed>>13)));
    const unsigned int churned2=0x1234567+(0x82948463*(churned1^(churned1>>20)));
    const unsigned int churned3=0x87654321^(0x87655677*(churned2^(churned2>>16)));
    /* The above lines are FREE since they'll be precomputed at compile time. They
       take user values like 1 2 3 and hash them to become roughly uncorrelated. */
    Saru z;

    z.state.x=churned2+state.x+(churned3^state.y);
    unsigned int add=(z.state.y+churned1)>>1;
    if (z.state.y-oWeylOffset<oWeylPeriod-add) z.state.y+=add;
    else z.state.y+=add-oWeylPeriod;
    return z;
}

/* Core PRNG evaluation. Very simple! */
template <unsigned int steps>
__device__ inline unsigned int Saru::u32()
{
    advanceLCG<steps>();
    advanceWeyl<steps>();
    unsigned int v=(state.x ^ (state.x>>26))+state.y;
    //  return v;  // not fully random (won't pass crush tests), but faster!
    return (v^(v>>20))*0x6957f5a7;
}

__device__ inline unsigned int Saru::u32()
{
    return u32<1>();
}


/* Floats have 23 bits of mantissa. We take 31 p-rand bits, cast to
   signed int and simply multiply to get the (0,1] range. We shift and cast
   to signedint to boost x86 conversion speed,
   see worley.com/mathgeek/floatconvert.html. */

template <unsigned int steps>
__device__ inline float Saru::f()
{
    //  return ((signed int)(u32<steps>()>>1))*(1.0f/0x80000000); 
    return u32<steps>()*(1.0f/4294967296.0f);

}

/* for a range that doesn't start at 0, we use the full 32 bits since
   we need to add an offset anyway. We still use the int cast method. */
template <unsigned int steps>
__device__ inline float Saru::f(float low, float high)
{
    const float TWO_N32 = 0.232830643653869628906250e-9f; /* 2^-32 */
    return ((signed int)(u32<steps>()))*(TWO_N32*(high-low))+0.5f*(high+low);
}


__device__ inline float Saru::f()
{
    return f<1>();
}

__device__ inline float Saru::f(float low, float high)
{
    return f<1>(low, high);
}


/* Doubles have 52 bits of mantissa. Casting to an int allows faster
   conversion, even with the extra needed fp addition. We use less-random
   "state" bits for the lowest order bits simply for speed. Output is
   in the (0,1] range. See worley.com/mathgeek/floatconvert.html. */

template <unsigned int steps>
__device__ inline double Saru::d()
{
    const double TWO_N32 = 0.232830643653869628906250e-9; /* 2^-32 */
    signed int v=(signed int)u32<steps>(); // deliberate cast to signed int for conversion speed

    return (v*TWO_N32+(0.5+0.5*TWO_N32))+((signed int)state.x)*(TWO_N32*TWO_N32);
}

template <unsigned int steps>
__device__ inline double Saru::d(double low, double high)
{
    const double TWO_N32 = 0.232830643653869628906250e-9; /* 2^-32 */
    signed int v=(signed int)u32<steps>(); // deliberate cast to signed int for conversion speed
    return (v*TWO_N32*(high-low)+(high+low)*(0.5+0.5*TWO_N32))+
        ((signed int)state.x)*(TWO_N32*TWO_N32*(high-low));
}


__device__ inline double Saru::d()
{
    return d<1>();
}

__device__ inline double Saru::d(double low, double high)
{
    return d<1>(low, high);
}

/////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////
__global__ void
ker_gen_rand(Saru * saruObj, double* d_idata, int size, int times)
{
    int threadGID = threadIdx.x+blockIdx.x*blockDim.x;
    
    if (threadGID < size) {
        for (int i = 0; i < times; i++)
            d_idata[threadGID + i * size] = saruObj[threadGID].d();
    }
}
/* Given the array of PRNG objects (created by a call to init_saru())
   Each call to this function will generate an array of "size" number of pseudo-random values (in double precision)
    + "size" is also the number of elements in the array saruObj and d_idata
    + d_idata is already allocated
*/
extern "C"  void gen_rand_saru_c1060(Saru * saruObj, double * d_idata, int size, int times){
    int threads = 192;
    int blocks = (size+threads-1)/threads;
    dim3 dimBlock(threads, 1, 1);
    dim3 dimGrid(blocks, 1, 1);

    ker_gen_rand<<< dimGrid, dimBlock>>>(saruObj, d_idata, size, times);
}
extern "C"  void gen_rand_saru_c2050(Saru * saruObj, double * d_idata, int size, int times){
    int threads = 256;
    int blocks = (size+threads-1)/threads;
    dim3 dimBlock(threads, 1, 1);
    dim3 dimGrid(blocks, 1, 1);

    ker_gen_rand<<< dimGrid, dimBlock>>>(saruObj, d_idata, size, times);
}


/* This function must be called by init_saru()
    + saruObj is an array of Saru objects; allocated on GPU
   This function provide the seed value to each object so that each one provide
   a unique pseudo-random sequences
    + global_seed is a global seed value
    + the local seed should be the global thread ID
*/   
__global__ void ker_init_saru(Saru *saruObj, int size, int global_seed) {
    int threadGID = threadIdx.x+blockIdx.x*blockDim.x;

    if (threadGID < size) {
        saruObj[threadGID].initialize(threadGID, global_seed);
    }
}

/*
  This function initialize the array of PRNG objects
  The array of PRNG objects will be reused until free by a call to finalize_saru()
 */
extern "C" Saru* init_saru(int size, int global_seed){
    Saru *saruObj;

    cudaMalloc((void**) &saruObj, size*sizeof(Saru));
    
    int threads = 192;
    int blocks = (size+threads-1)/threads;
    dim3 dimBlock(threads, 1, 1);
    dim3 dimGrid(blocks, 1, 1);
    
    ker_init_saru<<<dimGrid, dimBlock>>>(saruObj, size, global_seed);
    return saruObj;
}


/* This function free all saruObjs
   (must be called explicitly after using the PRNG)
 */
extern "C" void finalize_saru(Saru* saruObj) {
    cudaFree(saruObj);
}
