gb_rand.hxx

Go to the documentation of this file.

#ifndef CH_TOOLS__GB_RAND_HXX
#define CH_TOOLS__GB_RAND_HXX

#include <iostream>

namespace CH_Tools {


class GB_rand
{
private:
  long A[56];   
  int ind;      

  long mod_diff(long x,long y){return ((x-y)&0x7fffffff);}

  long flip_cycle ()
  {
    long *ii, *jj;
    for (ii = &A[1], jj = &A[32]; jj <= &A[55]; ii++, jj++)
      *ii = ((*ii-*jj) & 0x7fffffff);
    for (jj = &A[1]; ii <= &A[55]; ii++, jj++)
      *ii = ((*ii-*jj) & 0x7fffffff);
    ind = 54;
    return A[55];
  }
    
public:
  void init(long seed=1)
  {
    long i;
    long prev = seed, next = 1;
    seed = prev = (prev & 0x7fffffff);
    A[55] = prev;
    for (i = 21; i; i = (i + 21) % 55) {
      A[i] = next;
      
      next = ((prev-next) & 0x7fffffff);
      if (seed & 1)
        seed = 0x40000000 + (seed >> 1);
      else
        seed >>= 1;
      next = ((next- seed) & 0x7fffffff);
      
      prev = A[i];
    }
    
    (void) flip_cycle ();
    (void) flip_cycle ();
    (void) flip_cycle ();
    (void) flip_cycle ();
    (void) flip_cycle ();
    
  }
  
  GB_rand(long seed=1)
  {ind=0;A[0]=-1;init(seed);}

  ~GB_rand(){}
    
  long unif_long(long m)
  {
    const unsigned long two_to_the_31 = (unsigned long)0x80000000;
    unsigned long t = two_to_the_31 - (two_to_the_31 % (unsigned long)(m));
    long r;
    do {
      r = (A[ind]>=0?(A[ind--]):flip_cycle());
    } while (t <= (unsigned long) r);
    return r % m;
  }

  double next()
  {
    long r=unif_long(0x40000000);
    return (double(r)+.5)/double(0x40000000);
  }
  
  std::ostream& save(std::ostream& out) const
  {
    for(int i=0;i<56;i++) out<<A[i]<<"\n";
    out<<ind<<"\n"; return out;
  }

  std::istream& restore(std::istream& in)
  {
    for(int i=0;i<56;i++)
      in>>A[i];
    in>>ind; return in;
  }

};


}

#endif