Ticket #1155: Random.d

/*******************************************************************************
    Random number generators
    
    This is an attempt at having a good flexible and easy to use random number
    generator.
    ease of use:
    - shared generator for quick usage available through the "rand" object
      int i=rand.uniformR(10); // a random number from [0;10)
    - simple Random (non threadsafe) and RandomSync (threadsafe) types to 
      create new generators (for heavy use a good idea is one Random object per thread)
    - several distributions can be requested like this
        rand.distributionD!(type)(paramForDistribution)
      the type can often be avoided if the parameters make it clear.
      From it single numbers can be generated with .getRandom(), and variables
      initialized either with call style (var) or with .randomize(var).
      Utility functions to generate numbers directly are also available.
      The choice to put all the distribution in a single object that caches them
      has made (for example) the gamma distribution very easy to implement.
    - sample usage:
      auto r=new Random();
      int i; float f; real r; real[100] ar0; real[] ar=ar0[];
      // initialize with uniform distribution
      r(i)(f)(r)(ar);
      // unfortunetely one cannot use directly ar0...
      // unifrom distribution 0..10
      r.uniformR(10)(i)(f)(r)(ar);
      // uniform distribution -10..10
      r.uniformRSymm(10)(i)(f)(r)(ar);
      // the distribution can be stored
      auto r2=r.uniformRSymm(10);
      // and used later
      r2(ar);
      // complex distributions (normal,exp,gamma) are produced for the requested type
      r.normalSource!(float)()(f);
      // with sigma=2
      r.normalD(2.0f)(f);
      // and can be used also to initialize other types
      r.normalSource!(float)()(r)(ar);
      r.normalD(2.0f)(r)(ar);
      // but this is different from
      r.normalSource!(real)()(i)(r)(ar);
      r.normalD(2.0L)(i)(r)(ar);
      // as the source generates numbers of its type that then are simply casted to
      // the type needed.
      // Uniform distribution (as its creation for different types has no overhead)
      // never cast, so that (for example) bounds exclusion for floats is really
      // guaranteed.
      // For the other distribution using a distribution of different type than
      // the variable should be done with care, as underflow/overflow might ensue.
      //
      // some utility functions are also available
      int i2=r.uniform!(int)();
      int i2=r.randomize(i); // both i and i2 are initialized to the same value
      float f2=r.normalSigma(3.0f);
    flexibility:
    - easily swappable basic source
      // a random generator that uses the system provided random generator:
      auto r=RandomG!(Urandom)();
    - ziggurat generator can be easily adapted to any decreasing derivable
      distribution, the hard parametrization (to find xLast) can be done
      automatically
    - several distributions available "out of the box"
    quality:
    - the default Source (KISS99) passes all statistical tests 
      (P. L'Ecuyer and R. Simard, ACM Transactions on Mathematical Software (2007),
      33, 4, Article 22)
    - floating point uniform generator always initializes the full mantissa, the
      only flaw is a (*very* small) predilection of 0 as least important bit 
      (IEEE rounds to 0 in case of tie).
      Using a method that initializes the full mantissa was shown to improve the
      quality of subsequntly derived normal distribued numbers
      (Thomas et al. Gaussian random number generators. Acm Comput Surv (2007)
      vol. 39 (4) pp. 11))
    - Ziggurat method, a very fast and accurate method was used for both Normal and
      exp distribued numbers.
    - gamma distribued numbers uses a method recently proposed by Marsaglia and
      Tsang. The method is very fast, and should be good.
      My (Fawzi's) feeling is that the transformation h(x)=(1+d*x)^3 might loose
      a couple of bits of precision in some cases, but it is unclear if this
      might become visible in (*very* extensive) tests or not.
    - the basic source can be easily be changed with something else
    efficiency:
    - very fast methods have been used, and some effort has been put into
      optimizing some of them, but not all, but the interface has been choosen
      so that close to optimal implementation can be provided through the same
      interface.
    - Normal and Exp sources allocated only upon request: no memory waste, but
      a (*very* small) speed hit, that can be avoided by storing the source in
      a variable and using it (not going through the RandomG)
    annoyances:
    - I have put quite some functionality in a single module, maybe some might
      prefer to split it up
    - I have added two "next" methods to RandomG for backward compatibility
      reasons, and replaced the .hared instance from Random has been
      replaced by the "rand" object. The idea behind this is that RandomG is
      a template and rand it should be shared across all templates.
      If the name rand is considered bad one could change it. 
      I kept .shared static method that returns rand, so this remain a dropin
      replacement of the old random.
    - you cannot initialize a static array directly, this because randomize is
      declared like this:
        U randomize(U)(ref U a) { }
      and a static array cannot be passed by reference. Removing the ref would
      make arrays initialized, and scalar not, which is much worse.

        copyright:      Copyright (c) 2008. Fawzi Mohamed
        license:        BSD style: $(LICENSE)
        version:        Initial release: July 2008
        author:         Fawzi Mohamed

*******************************************************************************/
module tango.math.Random;
import tango.io.protocol.model.IWriter:IWritable,IWriter;
import tango.io.protocol.model.IReader:IReadable,IReader;
import Integer = tango.text.convert.Integer;
import tango.math.Bracket:findRoot;
import tango.core.sync.Mutex: Mutex;
import tango.io.FileConduit; // use stdc read/write?
import tango.math.Math:exp,sqrt,log,PI,abs;
//import tango.math.Probability: normalDistributionCompl;
import tango.math.ErrorFunction:erfc;
//import tango.stdc.stdio:printf;

version (Win32)
         private extern(Windows) int QueryPerformanceCounter (ulong *);

version (Posix)
        {
        private import tango.stdc.posix.sys.time;
        }

version(darwin) { version=has_urandom; }
version(linux)  { version=has_urandom; }

/// compile time integer power
private T powI(T)(T x,int p){
    T xx=cast(T)1;
    if (p<0){
        p=-p;
        x=1/x;
    }
    for (int i=0;i<p;++i)
        xx*=x;
    return xx;
}

/// if T is a float
template isFloat(T){
    static if(is(T==float)||is(T==double)||is(T==real)){
        const bool isFloat=true;
    } else {
        const bool isFloat=false;
    }
}

/// Strips the []'s off of a type.
template arrayBaseT(T)
{
    static if( is( T S : S[]) ) {
        alias arrayBaseT!(S)  arrayBaseT;
    }
    else {
        alias T arrayBaseT;
    }
}


/// very simple array based source (use with care, some methods in non uniform distributions
/// expect a random source with correct statistics, and could loop forever with such a source)
struct ArraySource{
    uint[] a;
    size_t i;
    static ArraySource opCall(uint[] a,size_t i=0)
    in { assert(a.length>0,"array needs at least one element"); }
    body {
        ArraySource res;
        res.a=a;
        res.i=i;
        return res;
    }
    uint next(){
        assert(a.length>i,"error, array out of bounds");
        uint el=a[i];
        i=(i+1)%a.length;
        return el;
    }
    ubyte nextB(){
        return cast(ubyte)(0xFF&next);
    }
    ulong nextL(){
        return ((cast(ulong)next)<<32)+cast(ulong)next;
    }
}

version(has_urandom) {
    /// basic source that takes data from system random device
    /// This is an engine, do not use directly, use RandomG!(Urandom)
    /// should use stdc rad/write?
    struct URandom{
        static FileConduit.Style readStyle;
        static Mutex lock;
        static this(){
            readStyle.access=FileConduit.Access.Read;
            readStyle.open  =FileConduit.Open.Exists;
            readStyle.share =FileConduit.Share.Read;
            readStyle.cache =FileConduit.Cache.None;

            lock=new Mutex();
        }
        const int canCheckpoint=false;
        const int canSeed=false;
    
        void skip(uint n){ }
        ubyte nextB(){
            union ToVoidA{
                ubyte i;
                void[1] a;
            }
            ToVoidA el;
            synchronized(lock){
                auto fn = new FileConduit("/dev/urandom", readStyle); 
                if(fn.read(el.a)!=el.a.length){
                    throw new Exception("could not write the requested bytes from urandom");
                }
                fn.close();
            }
            return el.i;
        }
        uint next(){
            union ToVoidA{
                uint i;
                void[4] a;
            }
            ToVoidA el;
            synchronized(lock){
                auto fn = new FileConduit("/dev/urandom", readStyle); 
                if(fn.read(el.a)!=el.a.length){
                    throw new Exception("could not write the requested bytes from urandom");
                }
                fn.close();
            }
            return el.i;
        }
        ulong nextL(){
            union ToVoidA{
                ulong l;
                void[8] a;
            }
            ToVoidA el;
            synchronized(lock){
                auto fn = new FileConduit("/dev/urandom", readStyle); 
                if(fn.read(el.a)!=el.a.length){
                    throw new Exception("could not write the requested bytes from urandom");
                }
                fn.close();
            }
            return el.l;
        }

        void seed(uint delegate() r) { }
        ///  IWritable implementation
        void write (IWriter input){
            char[] r="URandom";
            input(r);
        }
        /// IReadable implementation
        void read (IReader input){
            char[] r;
            input(r);
            assert(r=="URandom","unexpected read");
        }
        /// writes the current status in a string
        char[] toString(){
            return "URandom";
        }
        /// reads the current status from a string (that should have been trimmed)
        /// returns the number of chars read
        uint fromString(char[] s){
            char[] r="URandom";
            assert(s[0.. r.length]==r,"unxepected string instad of URandom:"~s);
            return r.length;
        }
    }
}

/+ Kiss99 random number generator, by Marisaglia
+ a simple RNG that passes all statistical tests
+ This is the engine, *never* use it directly, always use it though a RandomG class
+/
struct Kiss99{
    private uint kiss_x = 123456789;
    private uint kiss_y = 362436000;
    private uint kiss_z = 521288629;
    private uint kiss_c = 7654321;
    private uint nBytes = 0;
    private uint restB  = 0;
    
    const int canCheckpoint=true;
    const int canSeed=true;
    
    void skip(uint n){
        for (int i=n;i!=n;--i){
            next;
        }
    }
    ubyte nextB(){
        if (nBytes>0) {
            ubyte res=cast(ubyte)(restB & 0xFF);
            restB >>= 8;
            --nBytes;
            return res;
        } else {
            restB=next;
            ubyte res=cast(ubyte)(restB & 0xFF);
            restB >>= 8;
            nBytes=3;
            return res;
        }
    }
    uint next(){
        const ulong a = 698769069UL;
        ulong t;
        kiss_x = 69069*kiss_x+12345;
        kiss_y ^= (kiss_y<<13); kiss_y ^= (kiss_y>>17); kiss_y ^= (kiss_y<<5);
        t = a*kiss_z+kiss_c; kiss_c = (t>>32);
        kiss_z=cast(uint)t;
        return kiss_x+kiss_y+kiss_z;
    }
    ulong nextL(){
        return ((cast(ulong)next)<<32)+cast(ulong)next;
    }
    
    void seed(uint delegate() r){
        kiss_x = r();
        for (int i=0;i<100;++i){
            kiss_y=r();
            if (kiss_y!=0) break;
        }
        if (kiss_y==0) kiss_y=362436000;
        kiss_z=r();
        /* Don’t really need to seed c as well (is reset after a next),
           but doing it allows to completely restore a given internal state */
        kiss_c = r() % 698769069; /* Should be less than 698769069 */
        nBytes = 0;
        restB=0;
    }
    ///  IWritable implementation
    void write (IWriter input){
        input(kiss_x)(kiss_y)(kiss_z)(kiss_c)(nBytes)(restB);
    }
    /// IReadable implementation
    void read (IReader input){
        input(kiss_x)(kiss_y)(kiss_z)(kiss_c)(nBytes)(restB);
    }
    /// writes the current status in a string
    char[] toString(){
        char[] res=new char[6+6*9];
        int i=0;
        res[i..i+6]="KISS99";
        i+=6;
        foreach (val;[kiss_x,kiss_y,kiss_z,kiss_c,nBytes,restB]){
            res[i]='_';
            ++i;
            Integer.format(res[i..i+8],val,"x8");
            i+=8;
        }
        assert(i==res.length,"unexpected size");
        return res;
    }
    /// reads the current status from a string (that should have been trimmed)
    /// returns the number of chars read
    uint fromString(char[] s){
        uint i=0;
        assert(s[i..i+4]=="KISS","unexpected kind, expected KISS");
        assert(s[i+4..i+7]=="99_","unexpected version, expected 99");
        i+=6;
        foreach (val;[&kiss_x,&kiss_y,&kiss_z,&kiss_c,&nBytes,&restB]){
            assert(s[i]=='_',"no separator _ found");
            ++i;
            uint ate;
            *val=cast(uint)Integer.convert(s[i..i+8],16,&ate);
            assert(ate==8,"unexpected read size");
            i+=8;
        }
        return i;
    }
}

/+ like Kiss99, but synchronized, so multiple thread access is ok
+ (but if you need multiple thread access think about having a random number generator per thread)
+ This is the engine, *never* use it directly, always use it though a RandomG class
+/
struct Kiss99Sync{
    uint kiss_x = 123456789;
    uint kiss_y = 362436000;
    uint kiss_z = 521288629;
    uint kiss_c = 7654321;
    uint nBytes = 0;
    uint restB  = 0;
    Mutex lock;
    
    const int canCheckpoint=true;
    const int canSeed=true;
    
    void skip(uint n){
        for (int i=n;i!=n;--i){
            next;
        }
    }
    ubyte nextB(){
        if (nBytes>0) {
            ubyte res=cast(ubyte)(restB & 0xFF);
            restB >>= 8;
            --nBytes;
            return res;
        } else {
            restB=next;
            ubyte res=cast(ubyte)(restB & 0xFF);
            restB >>= 8;
            nBytes=3;
            return res;
        }
    }
    uint next(){
        uint res;
        const ulong a = 698769069UL;
        ulong t;
        synchronized(lock){
            kiss_x = 69069*kiss_x+12345;
            kiss_y ^= (kiss_y<<13); kiss_y ^= (kiss_y>>17); kiss_y ^= (kiss_y<<5);
            t = a*kiss_z+kiss_c; kiss_c = (t>>32);
            kiss_z=cast(uint)t;
            res=kiss_x+kiss_y+kiss_z;
        }
        return res;
    }
    ulong nextL(){
        return ((cast(ulong)next)<<32)+cast(ulong)next;
    }
    
    void seed(uint delegate() r){
        if (!lock) lock=new Mutex();
        synchronized(lock){
            kiss_x = r();
            for (int i=0;i<100;++i){
                kiss_y=r();
                if (kiss_y!=0) break;
            }
            if (kiss_y==0) kiss_y=362436000;
            kiss_z=r();
            /* Don’t really need to seed c as well (is reset after a next),
               but doing it allows to completely restore a given internal state */
            kiss_c = r() % 698769069; /* Should be less than 698769069 */
            nBytes=0;
            restB=0;
        }
    }
    ///  IWritable implementation
    void write (IWriter input){
        synchronized(lock){
            input(kiss_x)(kiss_y)(kiss_z)(kiss_c)(nBytes)(restB);
        }
    }
    /// IReadable implementation
    void read (IReader input){
        synchronized(lock){
            input(kiss_x)(kiss_y)(kiss_z)(kiss_c)(nBytes)(restB);
        }
    }
    /// writes the current status in a string
    char[] toString(){
        char[] res=new char[6+6*9];
        int i=0;
        res[i..i+6]="KISS99";
        i+=6;
        synchronized(lock){
            foreach (val;[kiss_x,kiss_y,kiss_z,kiss_c,nBytes,restB]){
                res[i]='_';
                ++i;
                Integer.format(res[i..i+8],val,"x8");
                i+=8;
            }
        }
        assert(i==res.length,"unexpected size");
        return res;
    }
    /// reads the current status from a string (that should have been trimmed)
    /// returns the number of chars read
    uint fromString(char[] s){
        uint i=0;
        assert(s[i..i+4]=="KISS","unexpected kind, expected KISS");
        assert(s[i+4..i+7]=="99_","unexpected version, expected 99");
        i+=6;
        synchronized(lock){
            foreach (val;[&kiss_x,&kiss_y,&kiss_z,&kiss_c,&nBytes,&restB]){
                assert(s[i]=='_',"no separator _ found");
                ++i;
                uint ate;
                *val=cast(uint)Integer.convert(s[i..i+8],16,&ate);
                assert(ate==8,"unexpected read size");
                i+=8;
            }
        }
        return i;
    }
}

/// ziggurat method for decreasing distributions.
/// Marsaglia, Tsang, Journal of Statistical Software, 2000
/// If has negative is true the distribution is assumed to be symmetric with respect to 0, 
/// otherwise it is assumed to be from 0 to infinity.
/// Struct based to avoid extra indirection when wrapped in a class (and it should be wrapped
/// in a class and not used directly).
/// Call style initialization avoided on purpose (this is a big structure, you don't want to return it)
struct Ziggurat(RandG,T,alias probDensityF,alias tailGenerator,bool hasNegative=true){
    static assert(isFloat!(T),T.stringof~" not acceptable, only floating point variables supported");
    const int nBlocks=256;
    T[nBlocks+1] posBlock;
    T[nBlocks+1] fVal;
    RandG r;
    alias Ziggurat SourceT;
    /// initializes the ziggurat
    static Ziggurat create(alias invProbDensityF, alias cumProbDensityFCompl)(RandG rGenerator,real xLast=-1.0L,bool check_error=true){
        /// function to find xLast
        real findXLast(real xLast){
            real v=xLast*probDensityF(xLast)+cumProbDensityFCompl(xLast);
            real fMax=probDensityF(0.0L);
            real pAtt=xLast;
            real fAtt=probDensityF(xLast);
            for (int i=nBlocks-2;i!=0;--i){
                fAtt+=v/pAtt;
                if (fAtt>fMax) return fAtt+(i-1)*fMax;
                pAtt=invProbDensityF(fAtt);
                assert(pAtt>=0,"invProbDensityF is supposed to return positive values");
            }
            return fAtt+v/pAtt-fMax;
        }
        void findBracket(ref real xMin,ref real xMax){
            real vMin=cumProbDensityFCompl(0.0L)/nBlocks;
            real pAtt=0.0L;
            for (int i=1;i<nBlocks;++i){
                pAtt+=vMin/probDensityF(pAtt);
            }
            real df=findXLast(pAtt);
            if (df>0) {
                // (most likely)
                xMin=pAtt;
                real vMax=cumProbDensityFCompl(0.0L);
                xMax=pAtt+vMax/probDensityF(pAtt);
            } else {
                xMax=pAtt;
                xMin=vMin/probDensityF(0.0L);
            }
        }
        if (xLast<=0){
            real xMin,xMax;
            findBracket(xMin,xMax);
            xLast=findRoot(&findXLast,xMin,xMax);
            // printf("xLast:%La => %La\n",xLast,findXLast(xLast));
        }
        Ziggurat res;
        with (res){
            r=rGenerator;
            real v=probDensityF(xLast)*xLast+cumProbDensityFCompl(xLast);
            real pAtt=xLast;
            real fMax=probDensityF(0.0L);
            posBlock[1]=cast(T)xLast;
            real fAtt=probDensityF(xLast);
            fVal[1]=cast(T)fAtt;
            for (int i=2;i<nBlocks;++i){
                fAtt+=v/pAtt;
                assert(fAtt<=fMax,"Ziggurat contruction shoot out");
                pAtt=invProbDensityF(fAtt);
                assert<