#ifndef __ISAAC_HPP
#define __ISAAC_HPP


/*

    C++ TEMPLATE VERSION OF Robert J. Jenkins Jr.'s
    ISAAC Random Number Generator.

    Ported from vanilla C to to template C++ class
    by Quinn Tyler Jackson on 16-23 July 1998.

        quinn@qtj.net

    The function for the expected period of this
    random number generator, according to Jenkins is:

        f(a,b) = 2**((a+b*(3+2^^a)-1)

        (where a is ALPHA and b is bitwidth)
        
    So, for a bitwidth of 32 and an ALPHA of 8,
    the expected period of ISAAC is:

        2^^(8+32*(3+2^^8)-1) = 2^^8295

    Jackson has been able to run implementations
    with an ALPHA as high as 16, or

        2^^2097263

*/

/*
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE 
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
*/

typedef unsigned int UINT32;
const UINT32 GOLDEN_RATIO = UINT32(0x9e3779b9);


template <UINT32 ALPHA = (8)>
class QTIsaac
{
  public:
   
    typedef unsigned char byte;
   
    struct randctx
    {
        randctx(void)
        {
            randrsl = new UINT32[N];
            randmem = new UINT32[N];
        }
      
        ~randctx(void)
        {
            delete [] randrsl;
            delete [] randmem;
        }
      
        UINT32 randcnt;
        UINT32* randrsl;
        UINT32* randmem;
        UINT32 randa;
        UINT32 randb;
        UINT32 randc;
	};
   
    QTIsaac(UINT32 a = 0, UINT32 b = 0, UINT32 c = 0);
    virtual ~QTIsaac(void);
   
    UINT32 rand(void);
    virtual void randinit(randctx* ctx, bool bUseSeed);
    virtual void srand(
		UINT32 a = 0, UINT32 b = 0, UINT32 c = 0, UINT32* s = NULL);
   
    enum {N = (1<<ALPHA)};
   
  protected:
   
     virtual void isaac(randctx* ctx);
   
     UINT32 ind(UINT32* mm, UINT32 x);
     void rngstep(
		UINT32 mix, UINT32& a, UINT32& b, UINT32*& mm, UINT32*& m, 
		UINT32*& m2, UINT32*& r, UINT32& x, UINT32& y);
     virtual void shuffle(
		UINT32& a, UINT32& b, UINT32& c, UINT32& d, UINT32& e, UINT32& f, 
		UINT32& g, UINT32& h);
   
  private:
    randctx m_rc;   
};


template<UINT32 ALPHA>
QTIsaac<ALPHA>::QTIsaac(UINT32 a, UINT32 b, UINT32 c) : m_rc()
{
    srand(a, b, c);
}


template<UINT32 ALPHA>
QTIsaac<ALPHA>::~QTIsaac(void)
{
    // DO NOTHING
}


template<UINT32 ALPHA>
void QTIsaac<ALPHA>::srand(UINT32 a, UINT32 b, UINT32 c, UINT32* s)
{
	for(int i = 0; i < N; i++)
	{
		m_rc.randrsl[i] = s != NULL ? s[i] : 0;
	}
   
	m_rc.randa = a;
	m_rc.randb = b;
	m_rc.randc = c;
   
	randinit(&m_rc, true);
}


template<UINT32 ALPHA>
inline UINT32 QTIsaac<ALPHA>::rand(void)
{
	return 0x7fffffff & (!m_rc.randcnt-- ? 
		(isaac(&m_rc), m_rc.randcnt=(N-1), m_rc.randrsl[m_rc.randcnt]) : 
		m_rc.randrsl[m_rc.randcnt]);
}


template<UINT32 ALPHA>
inline void QTIsaac<ALPHA>::randinit(randctx* ctx, bool bUseSeed)
{
	UINT32 a,b,c,d,e,f,g,h;
	int i;

	a = b = c = d = e = f = g = h = GOLDEN_RATIO;

	UINT32* m = (ctx->randmem);
	UINT32* r = (ctx->randrsl);

	if(!bUseSeed)
	{
		ctx->randa = 0;
		ctx->randb = 0;
		ctx->randc = 0;
	}

	// scramble it
	for(i=0; i < 4; ++i)         
	{
		shuffle(a,b,c,d,e,f,g,h);
	}

	if(bUseSeed) 
	{
		// initialize using the contents of r[] as the seed

		for(i=0; i < N; i+=8)
		{
			a+=r[i  ]; b+=r[i+1]; c+=r[i+2]; d+=r[i+3];
			e+=r[i+4]; f+=r[i+5]; g+=r[i+6]; h+=r[i+7];

			shuffle(a,b,c,d,e,f,g,h);

			m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
			m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
		}           

		//do a second pass to make all of the seed affect all of m

		for(i=0; i < N; i += 8)
		{
			a+=m[i  ]; b+=m[i+1]; c+=m[i+2]; d+=m[i+3];
			e+=m[i+4]; f+=m[i+5]; g+=m[i+6]; h+=m[i+7];

			shuffle(a,b,c,d,e,f,g,h);

			m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
			m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
		}
	}
	else
	{
		// fill in mm[] with messy stuff

		shuffle(a,b,c,d,e,f,g,h);

		m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
		m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;

	}

	isaac(ctx);         // fill in the first set of results 
	ctx->randcnt = N;   // prepare to use the first set of results 
}


template<UINT32 ALPHA>
inline UINT32 QTIsaac<ALPHA>::ind(UINT32* mm, UINT32 x)
{
	return (*(UINT32*)((byte*)(mm) + ((x) & ((N-1)<<2))));
}


template<UINT32 ALPHA>
inline void QTIsaac<ALPHA>::rngstep(UINT32 mix, UINT32& a, UINT32& b, UINT32*& mm, UINT32*& m, UINT32*& m2, UINT32*& r, UINT32& x, UINT32& y)
{
	x = *m;  
	a = (a^(mix)) + *(m2++); 
	*(m++) = y = ind(mm,x) + a + b; 
	*(r++) = b = ind(mm,y>>ALPHA) + x; 
}


template<UINT32 ALPHA>
inline void QTIsaac<ALPHA>::shuffle(UINT32& a, UINT32& b, UINT32& c, UINT32& d, UINT32& e, UINT32& f, UINT32& g, UINT32& h)
{ 
	a^=b<<11; d+=a; b+=c; 
	b^=c>>2;  e+=b; c+=d; 
	c^=d<<8;  f+=c; d+=e; 
	d^=e>>16; g+=d; e+=f; 
	e^=f<<10; h+=e; f+=g; 
	f^=g>>4;  a+=f; g+=h; 
	g^=h<<8;  b+=g; h+=a; 
	h^=a>>9;  c+=h; a+=b; 
}


template<UINT32 ALPHA>
inline void QTIsaac<ALPHA>::isaac(randctx* ctx)
{
	UINT32 x,y;

	UINT32* mm = ctx->randmem;
	UINT32* r  = ctx->randrsl;

	UINT32 a = (ctx->randa);
	UINT32 b = (ctx->randb + (++ctx->randc));

	UINT32* m    = mm; 
	UINT32* m2   = (m+(N/2));
	UINT32* mend = m2;

	for(; m<mend; )
	{
		rngstep((a<<13), a, b, mm, m, m2, r, x, y);
		rngstep((a>>6) , a, b, mm, m, m2, r, x, y);
		rngstep((a<<2) , a, b, mm, m, m2, r, x, y);
		rngstep((a>>16), a, b, mm, m, m2, r, x, y);
	}

	m2 = mm;

	for(; m2<mend; )
	{
		rngstep((a<<13), a, b, mm, m, m2, r, x, y);
		rngstep((a>>6) , a, b, mm, m, m2, r, x, y);
		rngstep((a<<2) , a, b, mm, m, m2, r, x, y);
		rngstep((a>>16), a, b, mm, m, m2, r, x, y);
	}

	ctx->randb = b;
	ctx->randa = a;
}


#endif // __ISAAC_HPP

