#ifndef _RAN4_
#define _RAN4_

#define NITER 4
void psdes(unsigned long& lword, unsigned long& irword)
// $-1òüPseudo-DESòý hashing of the 64-bit word (lword,irword). Both 32-bit arguments are returned-A
// hashed on all bits.
{
  unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;
  static unsigned long c1[NITER]={0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};
  static unsigned long c2[NITER]={0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};
  for (i=0;i<NITER;i++) {
    //    Perform niter iterations of DES logic,us ing a simpler (non-cryptographic) nonlinear function
    //    instead of DES$-1òùs.-A
    ia=(iswap=(irword)) ^ c1[i]; //The bit-rich constants c1 and (below) c2 guarantee lots of nonlinear mixing.
    itmpl = ia & 0xffff;
    itmph = ia >> 16;
    ib=itmpl*itmpl+ ~(itmph*itmph);
    irword=(lword) ^ (((ia = (ib >> 16) | ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);
    lword=iswap;
  }
}

float ran4(long& idum)
// Returns a uniform random deviate in the range 0.0 to 1.0,gener ated by pseudo-DES (DESlike)
// hashing of the 64-bit word (idums,idum),where idums was set by a previous call with
// negative idum. Also increments idum. Routine can be used to generate a random sequence
// by successive calls,le aving idum unaltered between calls; or it can randomly access the nth
// deviate in a sequence by calling with idum = n. Different sequences are initialized by calls with
// differing negative values of idum.
{
  unsigned long irword, itemp, lword;
  static long idums = 0;
  //  The hexadecimal constants jflone and jflmsk below are used to produce a floating number
  //    between 1. and 2. by bitwise masking. They are machine-dependent.
  static unsigned long jflone = 0x3f800000;
  static unsigned long jflmsk = 0x007fffff;
  if (idum < 0) { // Reset idums and prepare to return the first deviate in its sequence.
    idums = -idum; 
    idum=1;
  }
  irword=idum;
  lword=idums;
  psdes(lword,irword); //$(B!H(BPseudo-DES$(B!I(B encode the words.
  itemp=jflone | (jflmsk & irword); //Mask to a floating number between 1 and 2
  ++idum; 
  return (*reinterpret_cast<float*>(&itemp))-1.0; //Subtraction moves range to 0. to 1.
}

#endif // _RAN4_