#ifndef SARUPRNG_H
#define SARUPRNG_H

class Saru {
 public:

    Saru() {state.x=0x12345678; state.y=12345678;}
    __device__ inline Saru(unsigned int seed);
    __device__ Saru(unsigned int seed1, unsigned int seed2);
    __device__ Saru(unsigned int seed1, unsigned int seed2, unsigned int seed3);

    __device__ inline void initialize(unsigned int seed1, unsigned int seed2);
    __device__ inline void initialize(unsigned int seed);
    __device__ inline void initialize(unsigned int seed1, unsigned int seed2, unsigned int seed3);
    /* Efficient compile-time computed advancements */
    template <unsigned int steps> __device__ inline  void advance() 
        { advanceWeyl<steps>(); advanceLCG<steps>(); }
   
    // OK to advance negative steps in LCG, it's done mod 2^32 so it's
    // the same as advancing 2^32-1 steps, which is correct!
    template <unsigned int steps> __device__ inline void rewind() 
        { rewindWeyl<steps>(); advanceLCG<-steps>(); }
   
    /* Slower (but still reasonable) run-time computed advancements */
    __device__ inline void advance(unsigned int steps);

    template <unsigned int seed> __device__ Saru fork() const; 
   
    template <unsigned int steps> __device__ inline unsigned int u32();
    template <unsigned int steps> __device__ inline float f();
    template <unsigned int steps> __device__ inline double d();
    template <unsigned int steps> __device__ inline float f(float low, float high);
    template <unsigned int steps> __device__ inline double d(double low, double high);
    __device__ inline unsigned int u32();
    __device__ inline float f();
    __device__ inline double d();
    __device__ inline float f(float low, float high);
    __device__ inline double d(double low, double high);
   
 private:
   
    /* compile-time metaprograms to compute LCG and Weyl advancement */
   
   
    /* Computes A^N mod 2^32 */
    template<unsigned int A, unsigned int N> struct CTpow 
                                             { static const unsigned int value=(N&1?A:1)*CTpow<A*A, N/2>::value; };
    /* Specialization to terminate recursion: A^0 = 1 */
    template<unsigned int A> struct CTpow<A, 0> { static const unsigned int value=1; };
   
    /* CTpowseries<A,N> computes 1+A+A^2+A^3+A^4+A^5..+A^(N-1) mod 2^32.
       We do NOT use the more elegant formula (a^N-1)/(a-1) (see Knuth
       3.2.1), because it's more awkward to compute with the CPU's
       inherent mod 2^32.
      
       Based on recursion:
       g(A,n)= (1+A)*g(A*A, n/2);      if n is even
       g(A,n)= 1+A*(1+A)*g(A*A, n/2);  if n is ODD (since n/2 truncates)   */
    template<unsigned int A, unsigned int N>
        struct CTpowseries {
            static const unsigned int recurse=(1+A)*CTpowseries<A*A, N/2>::value;
            static const unsigned int value=  (N&1) ? 1+A*recurse : recurse;
        };
    template<unsigned int A>
        struct CTpowseries<A, 0> { static const unsigned int value=0; };  
    template<unsigned int A>
        struct CTpowseries<A, 1> { static const unsigned int value=1; };
   
   
    /* Compute A*B mod m.  Tricky only because of implicit 2^32 modulus.
       Uses recursion. 
       if A is even, then A*B mod m =  (A/2)*(B+B mod m) mod m.
       if A is odd,  then A*B mod m =  (B+((A/2)*(B+B mod m) mod m)) mod m.  */
    template <unsigned int A, unsigned int B, unsigned int m>
        struct CTmultmod {
            static const unsigned int temp= // (A/2)*(B*2) mod m
                CTmultmod< A/2, (B>=m-B ? B+B-m : B+B), m>::value;
            static const unsigned int value= A&1 ? ((B>=m-temp) ? B+temp-m: B+temp) : temp;
        };
    template <unsigned int B, unsigned int m>  /* terminate the recursion */
        struct CTmultmod<0, B, m> { static const unsigned int value=0; };
   
    template <unsigned int offset, unsigned int delta, 
        unsigned int modulus, unsigned int steps> 
        __device__ inline unsigned int advanceAnyWeyl(unsigned int);

   
    static const unsigned int LCGA=0x4beb5d59; // Full period 32 bit LCG
    static const unsigned int LCGC=0x2600e1f7; 
    static const unsigned int oWeylPeriod=0xda879add; // Prime period 3666320093
    static const unsigned int oWeylOffset=0x8009d14b;
    static const unsigned int oWeylDelta=(oWeylPeriod-0x80000000)+(oWeylOffset-0x80000000); // wraps mod 2^32
   
    /* Compile-time template function to efficently advance a state x with
       a LCG (mod 2^32) N steps.  Runtime, this all becomes a super-simple
       single multiply and add. */ 
    template <unsigned int steps> __device__ inline  void advanceLCG()
        { state.x=CTpow<LCGA, steps>::value*state.x+LCGC*CTpowseries<LCGA, steps>::value; }
   
    template <unsigned int steps> __device__ inline  void advanceWeyl()
        { state.y=advanceAnyWeyl<oWeylOffset, oWeylDelta, oWeylPeriod, steps>(state.y); }
   
    template <unsigned int steps> __device__ inline  void rewindWeyl()
        { state.y=advanceAnyWeyl<oWeylOffset, oWeylPeriod-oWeylDelta, 
                oWeylPeriod, steps>(state.y); }
      
    uint2 state;
};

extern "C" Saru* init_saru(int size, int global_seed);
extern "C" void gen_rand_saru_c1060(Saru * saruObj, double * d_idata, int size, int times);
extern "C" void gen_rand_saru_c2050(Saru * saruObj, double * d_idata, int size, int times);
extern "C" void finalize_saru(Saru* saruObj);

#endif /* SARUPRNG_H */

  
  
