Ticket #1155: patchRandom.txt

Index: tango/math/Random.d
===================================================================
--- tango/math/Random.d (revision 3630)
+++ tango/math/Random.d (working copy)
@@ -11,8 +11,10 @@
 *******************************************************************************/
 
 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;
 
-
 version (Win32)
          private extern(Windows) int QueryPerformanceCounter (ulong *);
 
@@ -21,6 +23,17 @@
         private import tango.stdc.posix.sys.time;
         }
 
+/// 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;
+}
 
 /******************************************************************************
 
@@ -41,9 +54,12 @@
 
         The period of KISS is thus 2^32*(2^32-1)*(2^63+2^32-1) > 2^127
 
+        Conversion to floating point generalizing Dr Jurgen A Doornik
+        "Conversion of high-period random numbers to floating point"
+
 ******************************************************************************/
 
-class Random
+final class Random : IWritable, IReadable
 {
         /**********************************************************************
 
@@ -94,7 +110,7 @@
 
         **********************************************************************/
 
-        final Random seed ()
+        Random seed ()
         {
                 ulong s;
 
@@ -117,7 +133,7 @@
 
         **********************************************************************/
 
-        final Random seed (uint seed)
+        Random seed (uint seed)
         {
                 kiss_x = seed | 1;
                 kiss_y = seed | 2;
@@ -133,7 +149,7 @@
 
         **********************************************************************/
 
-        final uint next ()
+        uint next ()
         {
                 kiss_x = kiss_x * 69069 + 1;
                 kiss_y ^= kiss_y << 13;
@@ -156,7 +172,7 @@
 
         **********************************************************************/
 
-        final uint next (uint max)
+        uint next (uint max)
         {
                 return next() % max;
         }
@@ -170,9 +186,334 @@
 
         **********************************************************************/
 
-        final uint next (uint min, uint max)
+        uint next (uint min, uint max)
         {
                 return next(max-min) + min;
         }
+
+        /**********************************************************************
+                Returns a value between 0 and 1, exclusive
+         maybe here a bit transfer would be better, as it is faster, but then
+         avoiding 0 is more difficult
+        **********************************************************************/
+        double nextFloat ()
+        {
+            static assert(float.mant_dig<=33,"float mantissa expected to be at most 33 (32 normalized) bit");
+            const float factS=powI(0.5f,float.mant_dig-1),shiftS=0.5f+0.5f*factS;
+            const uint maskS=(~0u)>>(33-float.mant_dig);
+            const int shift1=-(1<<(float.mant_dig-2));
+            uint nV=next;
+            int v1=cast(int)(nV & maskS)+shift1;
+            return (v1 * factS + shiftS);
+        }
+        /**********************************************************************
+                Returns a value between 0 and 1, exclusive
+        **********************************************************************/
+        double nextDouble ()
+        {
+            const double fact32=powI(0.5,32),factD=powI(0.5,double.mant_dig-1),
+                shiftD=0.5+0.5*factD;
+            const uint maskD=(~0u)>>(65-double.mant_dig);
+            static assert(32<double.mant_dig && double.mant_dig<=64,
+                "double mantissa expected to be between 32 and 64 bits");
+            return ((cast(int)next) * fact32 + shiftD +(next & maskD) * factD);
+        }
+        
+        /**********************************************************************
+                Returns a value between 0 and 1, exclusive
+        **********************************************************************/
+        real nextReal ()
+        {
+            static assert(33<real.mant_dig && real.mant_dig<=65,
+                "sizes larger than 65 bits or smaller than 33 (64 32 normalized) not suported");
+            const real fact32=powI(0.5L,32),factR=powI(0.5L,real.mant_dig-1),
+                shiftR=0.5L+0.5L*factR;
+            const uint maskR=(~0u)>>(65-real.mant_dig);
+            return ((cast(int)next) * fact32 + shiftR +(next & maskR) * factR);
+        }
+        /// randomizes the given array and returns it (for some types this is potentially
+        /// more efficient, both from the use of random numbers and speedwise)
+        byte[] randomize(byte[] a){
+            uint val=next; /// begin without value?
+            int rest=4;
+            for (size_t i=0;i<a.length;++i) {
+                if (rest!=0){
+                    a[i]=cast(byte)(0xFF&val);
+                    val>>=8;
+                    --rest;
+                } else {
+                    val=next;
+                    a[i]=cast(byte)(0xFF&val);
+                    val>>=8;
+                    rest=3;
+                }
+            }
+            return a;
+        }
+        ubyte[] randomize(ubyte[] a){
+            uint val=next; /// begin without value?
+            int rest=4;
+            for (size_t i=0;i<a.length;++i) {
+                if (rest!=0){
+                    a[i]=cast(ubyte)(0xFF&val);
+                    val>>=8;
+                    --rest;
+                } else {
+                    val=next;
+                    a[i]=cast(ubyte)(0xFF&val);
+                    val>>=8;
+                    rest=3;
+                }
+            }
+            return a;
+        }
+        uint[] randomize(uint[] a){
+            for (size_t i=0;i<a.length;++i)
+                a[i]=next;
+            return a;
+        }
+        /// ditto
+        int[] randomize(int[] a){
+            for (size_t i=0;i<a.length;++i)
+                a[i]=cast(int)next;
+            return a;
+        }
+        /// ditto
+        float[] randomize(float[] a){
+            static assert(float.mant_dig<=33,
+                "float mantissa expected to be at most 33 (32 normalized) bit");
+            const float factS=powI(0.5f,float.mant_dig-1),shiftS=0.5f+0.5f*factS;
+            const uint maskS=(~0u)>>(33-float.mant_dig),
+                maskRest=~maskS;
+            const int shift1=-(1<<(float.mant_dig-2));
+            const int nbitRest=33-float.mant_dig;
+            uint rest=0,nV;
+            int nrest=32-nbitRest;
+            foreach(ref el;a){
+                if (nrest<0){
+                    nrest=32-nbitRest;
+                    nV=rest;
+                } else {
+                    nV=next;
+                    rest|= maskRest & nV;
+                    rest>>=nbitRest;
+                    nrest-=nbitRest;
+                }
+                int v1=cast(int)(nV & maskS)+shift1;
+                el=(v1 * factS + shiftS);
+            }
+            return a;
+        }
+        /// ditto
+        double[] randomize(double[] a){
+            const double fact32=powI(0.5,32),factD=powI(0.5L,double.mant_dig-1),
+                shiftD=0.5L+0.5L*factD;
+            const uint maskD=(~0u)>>(65-double.mant_dig),
+                maskRest= ~maskD;
+            const int nbitRest=65-double.mant_dig;
+            uint rest=0,nV;
+            int nrest=32-nbitRest;
+            foreach (ref el;a){
+                if (nrest<0){
+                    nrest=32-nbitRest;
+                    nV=rest;
+                } else {
+                    nV=next;
+                    rest|= maskRest & nV;
+                    rest>>=nbitRest;
+                    nrest-=nbitRest;
+                }
+                el=((cast(int)next) * fact32 +(nV & maskD) * factD + shiftD);
+            }
+            return a;
+        }
+        /// ditto
+        real[] randomize(real[] a){
+            const real fact32=powI(0.5L,32),factR=powI(0.5L,real.mant_dig-1),
+                shiftR=0.5L+0.5L*factR;
+            const uint maskR=(~0u)>>(65-real.mant_dig),
+                maskRest= ~maskR;
+            static if (real.mant_dig>=63){
+                /// avoid all the overhead for just few bits...
+                foreach (ref el;a){
+                    el=((cast(int)next) * fact32 + shiftR +(maskR & next) * factR);
+                }
+            } else {
+                const int nbitRest=65-real.mant_dig;
+                uint rest=0,nV;
+                int nrest=32-nbitRest;
+                foreach (ref el;a){
+                    if (nrest<0){
+                        nrest=32-nbitRest;
+                        nV=rest;
+                    } else {
+                        nV=next;
+                        rest<<=nbitRest;
+                        rest|= maskRest & nV;
+                        nrest-=nbitRest;
+                    }
+                    el=((cast(int)next) * fact32 +(nV & maskR) * factR + shiftR);
+                }
+            }
+            return a;
+        }
+        
+        ///  IWritable implementation
+        void write (IWriter input){
+            input(kiss_k)(kiss_m)(kiss_x)(kiss_y)(kiss_z)(kiss_w)(kiss_carry);
+        }
+        /// IReadable implementation
+        void read (IReader input){
+            input(kiss_k)(kiss_m)(kiss_x)(kiss_y)(kiss_z)(kiss_w)(kiss_carry);
+        }
+        /// writes the current status in a string
+        char[] toString(){
+            char[] res=new char[12+7*8];
+            int i=0;
+            res[i..5]="KISS1";
+            i+=5;
+            foreach (val;[kiss_k,kiss_m,kiss_x,kiss_y,kiss_z,kiss_w,kiss_carry]){
+                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+6]=="1_","unexpected version, expected 1");
+            i+=5;
+            foreach (val;[&kiss_k,&kiss_m,&kiss_x,&kiss_y,&kiss_z,&kiss_w,&kiss_carry]){
+                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;
+        }
 }
 
+debug(UnitTest){
+    import tango.stdc.stdio:printf;
+    import tango.io.protocol.Writer;
+    import tango.io.protocol.Reader;
+    import tango.io.Buffer;
+
+    void printStat(T)(T[] a, bool checkB=false){
+        T minV,maxV;
+        real meanV=0.0L;
+        if (a.length>0){
+            minV=a[0];
+            maxV=a[0];
+            foreach (el;a){
+                assert(!checkB || (cast(T)0<el && el<cast(T)1),"el out of bounds");
+                if (el<minV) minV=el;
+                if (el>maxV) maxV=el;
+                meanV+=cast(real)el;
+            }
+            meanV/=cast(real)a.length;
+            bool printM=false;
+            if (checkB) {
+                if (minV>0.2) printM=true;
+                if (maxV<0.8) printM=true;
+                if (0.3>meanV || meanV>0.7) printM=true;
+            } else {
+                real offset=0.2L*(cast(real)T.max-cast(real)T.min);
+                real rMean=0.5L*(cast(real)T.max+cast(real)T.min);
+                if (cast(real)minV>cast(real)T.min+offset) printM=true;
+                if (cast(real)maxV<cast(real)T.max-offset) printM=true;
+                if (rMean-offset>meanV || meanV>rMean+offset) printM=true;
+            }
+            if (printM)
+                printf("WARNING tango.math.Random statistic is strange: %.*s[%d] %Le %Le %Le\n\0",cast(int)T.stringof.length,T.stringof.ptr,a.length,cast(real)minV,meanV,cast(real)maxV);
+        }
+    }
+
+    unittest {
+        // test float corner cases
+        uint zeros=0u,ones=~0u,nneg=0x8000_0000,ppos=0x7FFF_FFFF;
+        const float factS=powI(0.5f,float.mant_dig-1),shiftS=0.5f+0.5f*factS;
+        const uint maskS=(~0u)>>(33-float.mant_dig);
+        const int shift1=-1<<(float.mant_dig-2);
+        int v1=cast(int)(zeros & maskS)+shift1;
+        float f1=(v1 * factS + shiftS);
+        assert(f1>0.0f && f1<1.0f,"float conversion out of bounds");
+        v1=cast(int)(ones & maskS)+shift1;
+        float f2=(v1 * factS + shiftS);
+        assert(f2>0.0f && f2<1.0f,"float conversion out of bounds");
+        assert(1.0f-f2==f1,"float conversion of 0 and 0xFFFFFFFF, no symmetry");
+        assert(float.epsilon>f1,"minimum float bigger or equal epsilon");
+        // test double corner cases
+        const double fact32=powI(0.5,32),factD=powI(0.5,double.mant_dig-1),
+        shiftD=0.5+0.5*factD;
+        const uint maskD=(~0u)>>(65-double.mant_dig);
+        double d1=((cast(int)nneg) * fact32 + shiftD +(zeros & maskD) * factD);
+        double d2=((cast(int)ppos) * fact32 + shiftD +(ones & maskD) * factD);
+        assert(d1>0.0 && d2<1.0,"double conversion out of bounds");
+        assert(d2>0.0 && d2<1.0,"double conversion out of bounds");
+        assert(1.0-d2==d1,"double conversion of 0 and 0xFFFFFFFF, no symmetry");
+        assert(double.epsilon>d1,"minimum double bigger or equal epsilon");
+        // test real corner cases
+        const real fact32R=powI(0.5L,32),factR=powI(0.5L,real.mant_dig-1),
+            shiftR=0.5L+0.5L*factR;
+        const uint maskR=(~0u)>>(65-real.mant_dig);
+        real r1=((cast(int)nneg) * fact32R + shiftR +(zeros & maskR) * factR);
+        real r2=((cast(int)ppos) * fact32R + shiftR +(ones & maskR) * factR);
+        assert(r1>0.0L && r2<1.0L,"real conversion out of bounds");
+        assert(r2>0.0L && r2<1.0L,"real conversion out of bounds");
+        assert(1.0L-r2==r1,"real conversion of 0 and 0xFFFFFFFF, no symmetry");
+        assert(real.epsilon>r1,"minimum real bigger or equal epsilon");
+        Random r=new Random();
+        
+        // checkpoint status (str or writer)
+        char[] status=r.toString();
+        Buffer buf=new Buffer(100);
+        auto reader = new Reader (buf);
+        auto writer = new Writer (buf);
+        writer(r);
+        
+        uint tVal=r.next;
+        int count=10_000;
+        for (int i=count;i!=0;--i){
+            float f=r.nextFloat();
+            assert(0<f && f<1,"float out of bounds");
+            double d=r.nextDouble();
+            assert(0<d && d<1,"double out of bounds");
+            real rr=r.nextReal();
+            assert(0<rr && rr<1,"double out of bounds");
+        }
+        byte[1000] barr;
+        ubyte[1000] ubarr;
+        uint[1000] uiarr;
+        int[1000] iarr;
+        float[1000] farr;
+        double[1000] darr;
+        real[1000] rarr;
+        r.randomize(barr);
+        printStat(barr);
+        r.randomize(ubarr);
+        printStat(ubarr);
+        r.randomize(uiarr);
+        printStat(uiarr);
+        r.randomize(iarr);
+        printStat(iarr);
+        r.randomize(farr);
+        printStat(farr,true);
+        r.randomize(darr);
+        printStat(darr,true);
+        r.randomize(rarr);
+        printStat(rarr,true);
+        r.fromString(status);
+        assert(r.next==tVal,"restoring of status from str failed");
+        reader(r);
+        assert(r.next==tVal,"restoring of status from writer failed");
+    }
+
+}