[m-rev.] diff: monte

Zoltan Somogyi zs at csse.unimelb.edu.au
Mon Aug 9 12:27:36 AEST 2010


extras/monte/*:
	Add the files of this example program.

cvs diff: Diffing .
Index: ORIGINAL
===================================================================
RCS file: ORIGINAL
diff -N ORIGINAL
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ORIGINAL	9 Aug 2010 02:23:18 -0000
@@ -0,0 +1,302 @@
+Article: 16024 of sci.math.num-analysis
+Xref: taurus.cs.nps.navy.mil sci.stat.consult:7790 sci.math.num-analysis:16024
+Path: taurus.cs.nps.navy.mil!lll-winken.llnl.gov!uwm.edu!news.alpha.net!news.mathworks.com!udel!ssnet.com!usenet
+From: Bob Wheeler <bwheeler at ssnet.com>
+Newsgroups: sci.stat.consult,sci.math.num-analysis
+Subject: Marsaglia's Mother of all RNG's (Long?)
+Date: Fri, 28 Oct 94 19:32:08 EDT
+Organization: SSNet -- Public Internet Access in Delaware!
+Lines: 285
+Distribution: inet
+Message-ID: <38s2p1$qaf at marlin.ssnet.com>
+NNTP-Posting-Host: echip.com
+Mime-Version: 1.0
+Content-Type: TEXT/PLAIN; charset=US-ASCII
+X-Newsreader: NEWTNews & Chameleon -- TCP/IP for MS Windows from NetManage
+
+
+Several people have asked me to post this:
+First the C program, then George Marsaliga's post with details
+about the RNG.  He claims a period of about 2^250 for this and
+that it passes all of the usual tests.  I've tried it enough
+to be sure his claim is reasonable.
+
+The program:
+
+
+#include <string.h>
+
+static short mother1[10];
+static short mother2[10];
+static short mStart=1;
+
+#define m16Long 65536L             /* 2^16 */
+#define m16Mask 0xFFFF          /* mask for lower 16 bits */
+#define m15Mask 0x7FFF         /* mask for lower 15 bits */
+#define m31Mask 0x7FFFFFFF     /* mask for 31 bits */
+#define m32Double  4294967295.0  /* 2^32-1 */
+
+/* Mother **************************************************************
+|  George Marsaglia's The mother of all random number generators
+|      producing uniformly distributed pseudo random 32 bit values
+with
+|      period about 2^250.
+|  The text of Marsaglia's posting is appended at the end of the function.
+|
+|  The arrays mother1 and mother2 store carry values in their
+|      first element, and random 16 bit numbers in elements 1 to 8.
+|      These random numbers are moved to elements 2 to 9 and a new
+|      carry and number are generated and placed in elements 0 and 1.
+|  The arrays mother1 and mother2 are filled with random 16 bit values
+|      on first call of Mother by another generator.  mStart is the
+switch.
+|
+|  Returns:
+|  A 32 bit random number is obtained by combining the output of the
+|      two generators and returned in *pSeed.  It is also scaled by
+|      2^32-1 and returned as a double between 0 and 1
+|
+|  SEED:
+|  The inital value of *pSeed may be any long value
+|
+|  Bob Wheeler 8/8/94
+*/
+
+
+double Mother(unsigned long *pSeed)
+{
+   unsigned long  number,
+                       number1,
+                       number2;
+   short n,
+           *p;
+   unsigned short sNumber;
+
+       /* Initialize motheri with 9 random values the first time */
+   if (mStart) {
+       sNumber= *pSeed&m16Mask;   /* The low 16 bits */
+       number= *pSeed&m31Mask;   /* Only want 31 bits */
+
+       p=mother1;
+       for (n=18;n--;) {
+           number=30903*sNumber+(number>>16);   /* One line
+multiply-with-cary */
+           *p++=sNumber=number&m16Mask;
+           if (n==9)
+               p=mother2;
+       }
+       /* make cary 15 bits */
+       mother1[0]&=m15Mask;
+       mother2[0]&=m15Mask;
+       mStart=0;
+   }
+
+       /* Move elements 1 to 8 to 2 to 9 */
+   memmove(mother1+2,mother1+1,8*sizeof(short));
+   memmove(mother2+2,mother2+1,8*sizeof(short));
+
+       /* Put the carry values in numberi */
+   number1=mother1[0];
+   number2=mother2[0];
+
+       /* Form the linear combinations */
+
+number1+=1941*mother1[2]+1860*mother1[3]+1812*mother1[4]+1776*mother1[5]+
+
+1492*mother1[6]+1215*mother1[7]+1066*mother1[8]+12013*mother1[9];
+
+number2+=1111*mother2[2]+2222*mother2[3]+3333*mother2[4]+4444*mother2[5]+
+
+5555*mother2[6]+6666*mother2[7]+7777*mother2[8]+9272*mother2[9];
+
+       /* Save the high bits of numberi as the new carry */
+   mother1[0]=number1/m16Long;
+   mother2[0]=number2/m16Long;
+       /* Put the low bits of numberi into motheri[1] */
+   mother1[1]=m16Mask&number1;
+   mother2[1]=m16Mask&number2;
+
+       /* Combine the two 16 bit random numbers into one 32 bit */
+   *pSeed=(((long)mother1[1])<<16)+(long)mother2[1];
+
+       /* Return a double value between 0 and 1 */
+   return ((double)*pSeed)/m32Double;
+}
+
+
+
+/*
+
+*********************
+Marsaglia's comments
+
+        Yet another RNG
+Random number generators are frequently posted on
+the network; my colleagues and I posted ULTRA in
+1992 and, from the number of requests for releases
+to use it in software packages, it seems to be
+widely used.
+
+I have long been interested in RNG's and several
+of my early ones are used as system generators or
+in statistical packages.
+
+So why another one?  And why here?
+
+Because I want to describe a generator, or
+rather, a class of generators, so promising
+I am inclined to call it
+
+   The Mother of All Random Number Generators
+
+and because the generator seems promising enough
+to justify shortcutting the many months, even
+years, before new developments are widely
+known through publication in a journal.
+
+This new class leads to simple, fast programs that
+produce sequences with very long periods.  They
+use multiplication, which experience has shown
+does a better job of mixing bits than do +,- or
+exclusive-or, and they do it with easily-
+implemented arithmetic modulo a power of 2, unlike
+arithmetic modulo a prime.  The latter, while
+satisfactory, is difficult to implement.  But the
+arithmetic here modulo 2^16 or 2^32 does not suffer
+the flaws of ordinary congruential generators for
+those moduli: trailing bits too regular.  On the
+contrary, all bits of the integers produced by
+this new method, whether leading or trailing, have
+passed extensive tests of randomness.
+
+Here is an idea of how it works, using, say, integers
+of six decimal digits from which we return random 3-
+digit integers.  Start with n=123456, the seed.
+
+Then form a new n=672*456+123=306555 and return 555.
+Then form a new n=672*555+306=373266 and return 266.
+Then form a new n=672*266+373=179125 and return 125,
+
+and so on.  Got it?  This is a multiply-with-carry
+sequence x(n)=672*x(n-1)+ carry mod b=1000, where
+the carry is the number of b's dropped in the
+modular reduction. The resulting sequence of 3-
+digit x's has period 335,999.  Try it.
+
+No big deal, but that's just an example to give
+the idea. Now consider the sequence of 16-bit
+integers produced by the two C statements:
+
+k=30903*(k&65535)+(k>>16); return(k&65535);
+
+Notice that it is doing just what we did in the
+example: multiply the bottom half (by 30903,
+carefully chosen), add the top half and return the
+new bottom.
+
+That will produce a sequence of 16-bit integers
+with period > 2^29, and if we concatenate two
+such:
+     k=30903*(k&65535)+(k>>16);
+     j=18000*(j&65535)+(j>>16);
+     return((k<<16)+j);
+we get a sequence of more than 2^59 32-bit integers
+before cycling.
+
+The following segment in a (properly initialized)
+C procedure will generate more than 2^118
+32-bit random integers from six random seed values
+i,j,k,l,m,n:
+         k=30903*(k&65535)+(k>>16);
+         j=18000*(j&65535)+(j>>16);
+         i=29013*(i&65535)+(i>>16);
+         l=30345*(l&65535)+(l>>16);
+         m=30903*(m&65535)+(m>>16);
+         n=31083*(n&65535)+(n>>16);
+         return((k+i+m)>>16)+j+l+n);
+
+And it will do it much faster than any of several
+widely used generators designed to use 16-bit
+integer arithmetic, such as that of Wichman-Hill
+that combines congruential sequences for three
+15-bit primes (Applied Statistics, v31, p188-190,
+1982), period about 2^42.
+
+I call these multiply-with-carry generators. Here
+is an extravagant 16-bit example that is easily
+implemented in C or Fortran. It does such a
+thorough job of mixing the bits of the previous
+eight values that it is difficult to imagine a
+test of randomness it could not pass:
+
+x[n]=12013x[n-8]+1066x[n-7]+1215x[n-6]+1492x[n-5]+1776x[n-4]
+ +1812x[n-3]+1860x[n-2]+1941x[n-1]+carry mod 2^16.
+
+The linear combination occupies at most 31 bits of
+a 32-bit integer. The bottom 16 is the output, the
+top 15 the next carry. It is probably best to
+implement with 8 case segments. It takes 8
+microseconds on my PC. Of course it just provides
+16-bit random integers, but awfully good ones. For
+32 bits you would have to combine it with another,
+such as
+
+x[n]=9272x[n-8]+7777x[n-7]+6666x[n-6]+5555x[n-5]+4444x[n-4]
+    +3333x[n-3]+2222x[n-2]+1111x[n-1]+carry mod 2^16.
+
+Concatenating those two gives a sequence of 32-bit
+random integers (from 16 random 16-bit seeds),
+period about 2^250. It is so awesome it may merit
+the Mother of All RNG's title.
+
+The coefficients in those two linear combinations
+suggest that it is easy to get long-period
+sequences, and that is true.  The result is due to
+Cemal Kac, who extended the theory we gave for
+add-with-carry sequences: Choose a base b and give
+r seed values x[1],...,x[r] and an initial 'carry'
+c. Then the multiply-with-carry sequence
+
+ x[n]=a1*x[n-1]+a2*x[n-2]+...+ar*x[n-r]+carry mod b,
+
+where the new carry is the number of b's dropped
+in the modular reduction, will have period the
+order of b in the group of residues relatively
+prime to m=ar*b^r+...+a1b^1-1.  Furthermore, the
+x's are, in reverse order, the digits in the
+expansion of k/m to the base b, for some 0<k<m.
+
+In practice b=2^16 or b=2^32 allows the new
+integer and the new carry to be the bottom and top
+half of a 32- or 64-bit linear combination of  16-
+or 32-bit integers.  And it is easy to find
+suitable m's if you have a primality test:  just
+search through candidate coefficients until you
+get an m that is a safeprime---both m and (m-1)/2
+are prime.  Then the period of the multiply-with-
+carry sequence will be the prime (m-1)/2. (It
+can't be m-1 because b=2^16 or 2^32 is a square.)
+
+Here is an interesting simple MWC generator with
+period> 2^92, for 32-bit arithmetic:
+
+x[n]=1111111464*(x[n-1]+x[n-2]) + carry mod 2^32.
+
+Suppose you have functions, say top() and bot(),
+that give the top and bottom halves of a 64-bit
+result.  Then, with initial 32-bit x, y and carry
+c,  simple statements such as
+     y=bot(1111111464*(x+y)+c)
+     x=y
+     c=top(y)
+will, repeated, give over 2^92 random 32-bit y's.
+
+Not many machines have 64 bit integers yet.  But
+most assemblers for modern CPU's permit access to
+the top and bottom halves of a 64-bit product.
+
+I don't know how to readily access the top half of
+a 64-bit product in C.  Can anyone suggest how it
+might be done? (in integer arithmetic)
+
+George Marsaglia geo at stat.fsu.edu
Index: doit.m
===================================================================
RCS file: doit.m
diff -N doit.m
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ doit.m	6 Aug 2010 04:02:01 -0000
@@ -0,0 +1,29 @@
+%-----------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sw=4 et
+%-----------------------------------------------------------------------------%
+
+:- module doit.
+
+:- interface.
+
+:- import_module io.
+
+:- pred main(io::di, io::uo) is det.
+
+:- implementation.
+
+:- import_module hg.
+:- import_module geom.
+:- import_module monte.
+:- import_module rnd.
+
+:- import_module float.
+:- import_module int.
+:- import_module list.
+:- import_module string.
+
+main(!IO) :-
+    NumPoints = 1000000,
+    monte(hgBox, hg, NumPoints, Volume),
+    io.format("The computed volume after testing %d points is %f\n",
+        [i(NumPoints), f(Volume)], !IO).
Index: dots.m
===================================================================
RCS file: dots.m
diff -N dots.m
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ dots.m	6 Aug 2010 04:04:55 -0000
@@ -0,0 +1,55 @@
+%-----------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sw=4 et
+%-----------------------------------------------------------------------------%
+
+:- module dots.
+
+:- interface.
+
+:- import_module io.
+
+:- pred main(io::di, io::uo) is det.
+
+:- implementation.
+
+:- import_module hg.
+:- import_module geom.
+:- import_module monte.
+:- import_module rnd.
+
+:- import_module float.
+:- import_module int.
+:- import_module list.
+:- import_module string.
+
+main(!IO) :-
+    io.format("newgraph\n", [], !IO),
+    io.format("xaxis min -2.5 max 2.5 size %2.2f\n", [f(5.0/3.0)], !IO),
+    io.format("yaxis min -2.5 max 4.5 size %2.2f\n", [f(7.0/3.0)], !IO),
+    io.format("newcurve marksize 0.01 0.01 pts ", [], !IO),
+    dots(hgBox, hg, 100000, !IO).
+
+:- pred dots(box::in, monte.shape::in(monte.shape), precision::in,
+    io::di, io::uo) is det.
+
+dots(Box, Shape, Precision, !IO) :-
+    rnd__init(17, Rnd0),
+    dots2(Box, Shape, Precision, Rnd0, !IO).
+
+:- pred dots2(box::in, monte.shape::in(monte.shape), int::in, rnd::in,
+    io::di, io::uo) is det.
+
+dots2(Box, Shape, N, !.Rnd, !IO) :-
+    ( N > 0 ->
+        frange(Box ^ xmin, Box ^ xmax, X, !Rnd),
+        frange(Box ^ ymin, Box ^ ymax, Y, !Rnd),
+        frange(Box ^ zmin, Box ^ zmax, Z, !Rnd),
+        ( call(Shape, X, Y, Z) ->
+            io.format("%2.2f %2.2f ", [f(X), f(Z)], !IO)
+        ;
+            true
+        ),
+        dots2(Box, Shape, N - 1, !.Rnd, !IO)
+    ;
+        true
+    ).
Index: geom.m
===================================================================
RCS file: geom.m
diff -N geom.m
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ geom.m	6 Aug 2010 07:56:26 -0000
@@ -0,0 +1,124 @@
+%-----------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sw=4 et
+%-----------------------------------------------------------------------------%
+
+:- module geom.
+
+:- interface.
+
+:- import_module float.
+
+:- type shape == pred(float, float, float).
+:- inst shape == (pred(in, in, in) is semidet).
+
+:- pred sphere(float::in, float::in, float::in) is semidet.
+
+:- pred cylinder(float::in, float::in, float::in) is semidet.
+
+:- pred cube(float::in, float::in, float::in) is semidet.
+
+:- pred cone(float::in, float::in, float::in) is semidet.
+
+:- pred torus(float::in, float::in, float::in, float::in) is semidet.
+
+:- pred planeX(float::in, float::in, float::in) is semidet.
+
+:- pred planeY(float::in, float::in, float::in) is semidet.
+
+:- pred planeZ(float::in, float::in, float::in) is semidet.
+
+:- pred union(shape::in(shape), shape::in(shape),
+    float::in, float::in, float::in) is semidet.
+
+:- pred intersection(shape::in(shape), shape::in(shape),
+    float::in, float::in, float::in) is semidet.
+
+:- pred difference(shape::in(shape), shape::in(shape),
+    float::in, float::in, float::in) is semidet.
+
+:- pred translate(shape::in(shape), float::in, float::in, float::in,
+    float::in, float::in, float::in) is semidet.
+
+:- pred scale(shape::in(shape), float::in, float::in, float::in,
+    float::in, float::in, float::in) is semidet.
+
+:- pred rotateX(shape::in(shape), float::in,
+    float::in, float::in, float::in) is semidet.
+
+:- pred rotateY(shape::in(shape), float::in,
+    float::in, float::in, float::in) is semidet.
+
+:- pred rotateZ(shape::in(shape), float::in,
+    float::in, float::in, float::in) is semidet.
+
+:- implementation.
+
+:- import_module math.
+
+sphere(X, Y, Z) :-
+    sqr(X) + sqr(Y) + sqr(Z) =< 1.0.
+
+cylinder(X, Y, Z) :-
+    sqr(X) + sqr(Y) =< 1.0,
+    0.0 =< Z, Z =< 1.0.
+
+cube(X, Y, Z) :-
+    0.0 =< X, X =< 1.0,
+    0.0 =< Y, Y =< 1.0,
+    0.0 =< Z, Z =< 1.0.
+
+cone(X, Y, Z) :-
+    0.0 =< Z, Z =< 1.0,
+    sqr(X) + sqr(Y) =< 1.0 - Z.
+
+torus(R, X, Y, Z) :-
+    sqr(sqrt(sqr(X) + sqr(Z)) - 1.0) + sqr(Y) =< sqr(R).
+
+planeX(X, _, _) :-
+    X >= 0.0.
+
+planeY(_, Y, _) :-
+    Y >= 0.0.
+
+planeZ(_, _, Z) :-
+    Z >= 0.0.
+
+union(A, B, X, Y, Z) :-
+    (
+        call(A, X, Y, Z)
+    ;
+        call(B, X, Y, Z)
+    ).
+
+intersection(A, B, X, Y, Z) :-
+    call(A, X, Y, Z),
+    call(B, X, Y, Z).
+
+difference(A, B, X, Y, Z) :-
+    call(A, X, Y, Z),
+    not call(B, X, Y, Z).
+
+translate(Shape, Dx, Dy, Dz, X, Y, Z) :-
+    call(Shape, X - Dx, Y - Dy, Z - Dz).
+
+scale(Shape, Sx, Sy, Sz, X, Y, Z) :-
+    call(Shape, X / Sx, Y / Sy, Z / Sz).
+
+rotateX(Shape, Th, X, Y, Z) :-
+    Sth = sin(-Th),
+    Cth = cos(-Th),
+    call(Shape, X, Cth * Y - Sth * Z, Sth * Y + Cth * Z).
+
+rotateY(Shape, Th, X, Y, Z) :-
+    Sth = sin(-Th),
+    Cth = cos(-Th),
+    call(Shape, Cth * X + Sth * Z, Y, -Sth * X + Cth * Z).
+
+rotateZ(Shape, Th, X, Y, Z) :-
+    Sth = sin(-Th),
+    Cth = cos(-Th),
+    call(Shape, Cth * X - Sth * Y, Sth * X + Sth * Y, Z).
+
+:- func sqr(float) = float.
+
+sqr(X) = X * X.
Index: hg.m
===================================================================
RCS file: hg.m
diff -N hg.m
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ hg.m	6 Aug 2010 07:56:43 -0000
@@ -0,0 +1,106 @@
+%-----------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sw=4 et
+%-----------------------------------------------------------------------------%
+
+:- module hg.
+
+:- interface.
+
+:- import_module monte.
+:- import_module float.
+
+:- func hgBox = box.
+
+:- pred hg(float::in, float::in, float::in) is semidet.
+
+:- pred head(float::in, float::in, float::in) is semidet.
+
+:- pred body(float::in, float::in, float::in) is semidet.
+
+:- implementation.
+
+:- import_module geom.
+
+hgBox = box(-2.5, -2.0, -2.5, 2.5, 2.0, 4.5).
+
+hg(X, Y, Z) :-
+    union(
+        translate(head, 0.0, 0.0, 2.0),
+    body,
+    X, Y, Z).
+
+head(X, Y, Z) :-
+    union(
+        torus(0.25),
+        translate(
+            union(
+                difference(torus(0.25), planeZ),
+                union(
+                    translate(
+                        scale(sphere, 0.25, 0.25, 0.25),
+                        -1.0, 0.0, 0.0
+                    ),
+                    translate(
+                        scale(sphere, 0.25, 0.25, 0.25),
+                        1.0, 0.0, 0.0
+                    )
+                )
+            ),
+            0.0,
+            0.0,
+            2.0
+        ),
+        X,
+        Y,
+        Z
+    ).
+
+body(X, Y, Z) :-
+    union(
+        translate(
+            scale(cylinder, 0.25, 0.25, 3.0),
+            0.0,
+            0.0,
+            -2.0
+        ),
+        union(
+            rotateY(
+                translate(
+                    scale(cylinder, 0.25, 0.25, 2.0),
+                    0.0,
+                    0.0,
+                    -1.0
+                ),
+                pi / 2.0
+            ),
+            union(
+                translate(
+                    scale(sphere, 0.25, 0.25, 0.25),
+                    0.0,
+                    0.0,
+                    -2.0
+                ),
+                union(
+                    translate(
+                        scale(sphere, 0.25, 0.25, 0.25),
+                        -1.0,
+                        0.0,
+                        0.0
+                    ),
+                    translate(
+                        scale(sphere, 0.25, 0.25, 0.25),
+                        1.0,
+                        0.0,
+                        0.0
+                    )
+                )
+            )
+        ),
+        X,
+        Y,
+        Z
+    ).
+
+:- func pi = float.
+
+pi = 3.14159265359.
Index: monte.m
===================================================================
RCS file: monte.m
diff -N monte.m
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ monte.m	6 Aug 2010 07:57:27 -0000
@@ -0,0 +1,60 @@
+%-----------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sw=4 et
+%-----------------------------------------------------------------------------%
+
+:- module monte.
+
+:- interface.
+
+:- import_module float, int.
+
+:- type box
+    --->    box(
+            xmin :: float,
+            ymin :: float,
+            zmin :: float,
+            xmax :: float,
+            ymax :: float,
+            zmax :: float
+        ).
+
+:- type shape == pred(float, float, float).
+:- inst shape == (pred(in, in, in) is semidet).
+
+:- type precision == int.
+
+:- type volume == float.
+
+:- pred monte(box::in, shape::in(shape), precision::in, volume::out) is det.
+
+:- implementation.
+
+:- import_module rnd.
+
+monte(Box, Shape, Precision, Volume) :-
+    rnd.init(17, Rnd0),
+    monte_loop(Box, Shape, Precision, Rnd0, 0, Hits),
+    Prop = float(Hits) / float(Precision),
+    BoxVolume =
+        (Box ^ xmax - Box ^ xmin) *
+        (Box ^ ymax - Box ^ ymin) *
+        (Box ^ zmax - Box ^ zmin),
+    Volume = BoxVolume * Prop.
+
+:- pred monte_loop(box::in, shape::in(shape), int::in, rnd::in,
+    int::in, int::out) is det.
+
+monte_loop(Box, Shape, N, !.Rnd, !Hits) :-
+    ( N > 0 ->
+        frange(Box ^ xmin, Box ^ xmax, X, !Rnd),
+        frange(Box ^ ymin, Box ^ ymax, Y, !Rnd),
+        frange(Box ^ zmin, Box ^ zmax, Z, !Rnd),
+        ( call(Shape, X, Y, Z) ->
+            !:Hits = !.Hits + 1
+        ;
+            true
+        ),
+        monte_loop(Box, Shape, N - 1, !.Rnd, !Hits)
+    ;
+        true
+    ).
Index: rnd.m
===================================================================
RCS file: rnd.m
diff -N rnd.m
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ rnd.m	6 Aug 2010 07:55:38 -0000
@@ -0,0 +1,313 @@
+%-----------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sw=4 et
+%-----------------------------------------------------------------------------%
+% file: rnd.m
+% main author: conway.
+%
+% This pseudo-random number generator is derived from C source due to
+% George Marsaglia, geo at stat.fsu.edu.
+%
+% The original source is included in the file ORIGINAL as a reference.
+%
+%-----------------------------------------------------------------------------%
+
+:- module rnd.
+
+:- interface.
+
+:- import_module float.
+:- import_module int.
+:- import_module list.
+:- import_module pair.
+
+:- type rnd.
+
+    % Initialize a random number supply.
+    %
+:- pred rnd.init(int::in, rnd::out) is det.
+
+    % rnd(Random, !Rnd):
+    %
+    % Get a random float on the range [0, 1).
+    %
+:- pred rnd(float::out, rnd::in, rnd::out) is det.
+
+    % irange(Lower, Upper, Random, !Rnd):
+    %
+    % Get a random int on the range [Lower, Upper].
+    %
+:- pred irange(int::in, int::in, int::out, rnd::in, rnd::out) is det.
+
+    % frange(Lower, Upper, Random, !Rnd):
+    %
+    % Get a random float on the range [Lower, Upper).
+    %
+:- pred frange(float::in, float::in, float::out, rnd::in, rnd::out) is det.
+
+    % Generate a random permutation of a list.
+    %
+:- pred shuffle(list(T)::in, list(T)::out, rnd::in, rnd::out) is det.
+
+    % Get a random element of a list.
+    %
+:- pred oneof(list(T)::in, T::out, rnd::in, rnd::out) is det.
+
+    % get a random element of a weighted list.
+    %
+:- pred wghtd_oneof(list(pair(int, T))::in, T::out, rnd::in, rnd::out) is det.
+
+    % gaussian(X, Y, !Rnd):
+    %
+    % Generate a pair of gaussian deviates `X' and `Y'.
+    %
+:- pred gaussian(float::out, float::out, rnd::in, rnd::out) is det.
+
+%-----------------------------------------------------------------------------%
+
+:- implementation.
+
+:- import_module math.
+:- import_module require.
+
+irange(Min, Max, Val, !R) :-
+    frange(rfloat(Min), rfloat(Max+1), FVal, !R),
+    Val = rint(FVal).
+
+frange(Min, Max, Val, !R) :-
+    rnd(J, !R),
+    Val = J*(Max - Min)+Min.
+
+shuffle(Ins, Outs, !R) :-
+    list.length(Ins, N),
+    shuffle2(N, Ins, [], T0, !R),
+    shuffle2(N, T0, [], T1, !R),
+    shuffle2(N, T1, [], T2, !R),
+    shuffle2(N, T2, [], T3, !R),
+    shuffle2(N, T3, [], U, !R),
+    shuffle2(N, U, [], Outs, !R).
+
+:- pred shuffle2(int::in, list(T)::in, list(T)::in, list(T)::out,
+    rnd::in, rnd::out) is det.
+
+shuffle2(N, Ins, !Acc, !R) :-
+    ( N > 0 ->
+        irange(0, N - 1, J, !R),
+        delnth(Ins, J, Rest, T),
+        shuffle2(N - 1, Rest, [T | !.Acc], !:Acc, !R)
+    ;
+        true
+    ).
+
+:- pred delnth(list(T)::in, int::in, list(T)::out, T::out) is det.
+
+delnth([], _, _, _) :-
+    error("delnth: no enough elems!").
+delnth([X | Xs], N, Zs, Z) :-
+    ( N =< 0 ->
+        Z = X,
+        Zs = Xs
+    ;
+        Zs = [X|Ys],
+        delnth(Xs, N-1, Ys, Z)
+    ).
+
+oneof(Things, Thing, !R) :-
+    list.length(Things, Num),
+    irange(0, Num-1, X, !R),
+    list.index0_det(Things, X, Thing).
+
+wghtd_oneof(WghtdThings, Thing, !R) :-
+    cumu(WghtdThings, 0, Sum),
+    irange(0, Sum, X, !R),
+    pick(WghtdThings, X, Thing).
+
+:- pred cumu(list(pair(int, T))::in, int::in, int::out) is det.
+
+cumu([], Sum, Sum).
+cumu([X - _T | Rest0], Sum, Sum1) :-
+    cumu(Rest0, X+Sum, Sum1).
+
+:- pred pick(list(pair(int, T))::in, int::in, T::out) is det.
+
+pick([], _, _) :-
+    error("pick: no things to pick from!").
+pick([N - T | Rest], P, Thing) :-
+    ( N >= P ->
+        Thing = T
+    ;
+        pick(Rest, P - N, Thing)
+    ).
+
+gaussian(X, Y, !Rnd) :-
+    frange(-1.0, 1.0, V1, !Rnd),
+    frange(-1.0, 1.0, V2, !Rnd),
+    R = V1*V1 + V2*V2,
+    ( R >= 1.0, R \= 0.0  ->
+        gaussian(X, Y, !Rnd)
+    ;
+        Fac = sqrt(-2.0 * ln(R) / R),
+        X = V1 * Fac,
+        Y = V2 * Fac
+    ).
+
+%-----------------------------------------------------------------------------%
+
+:- type vec
+    --->    vec(int, int, int, int, int, int, int, int, int, int).
+
+:- type rnd
+    --->    rnd(
+                vec,
+                vec,
+                int
+            ).
+
+rnd.init(Seed, rnd(M1, M2, Seed)) :-
+    SN = Seed /\ ((1 `unchecked_left_shift` 15) - 1),
+    N  = Seed /\ ((1 `unchecked_left_shift` 30) - 1),
+    M1a = vec(0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
+    M2a = vec(0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
+    seed1(17, SN, N, M1a, M2a, M1b, M2b),
+    set(M1b, 0, (M1b ** 0) /\ ((1 `unchecked_left_shift` 15) - 1), M1),
+    set(M2b, 0, (M2b ** 0) /\ ((1 `unchecked_left_shift` 15) - 1), M2).
+
+:- pred seed1(int, int, int, vec, vec, vec, vec).
+:- mode seed1(in, in, in, in, in, out, out) is det.
+
+seed1(N, SNum0, Num0, M1a, M2a, M1, M2) :-
+    (N > 0 ->
+        Num1 = 30903 * SNum0 + (Num0 `unchecked_right_shift` 15),
+        SNum1 = Num1 /\ ((1 `unchecked_left_shift` 15) - 1),
+        ( N >= 9 ->
+            M2b = M2a,
+            set(M1a, 17 - N, SNum1, M1b)
+        ;
+            M1b = M1a,
+            set(M2a, 8 - N, SNum1, M2b)
+        ),
+        seed1(N-1, SNum1, Num1, M1b, M2b, M1, M2)
+    ;
+        M1 = M1a,
+        M2 = M2a
+    ).
+
+rnd(Res, rnd(M1a, M2a, _Seed0), rnd(M1d, M2d, Seed1)) :-
+    shift(M1a, M1b),
+    shift(M2a, M2b),
+    N1a = (M1b ** 0),
+    N2a = (M2b ** 0),
+
+    N1b = N1a + 1941 * (M1b ** 2) + 1860 * (M1b ** 3) + 1812 * (M1b ** 4)
+        + 1776 * (M1b ** 5) + 1492 * (M1b ** 6) + 1215 * (M1b ** 7)
+        + 1066 * (M1b ** 8) + 12013 * (M1b ** 9),
+
+    N2b = N2a + 1111 * (M2b ** 2) + 2222 * (M2b ** 3) + 3333 * (M2b ** 4)
+        + 4444 * (M2b ** 5) + 5555 * (M2b ** 6) + 6666 * (M2b ** 7)
+        + 7777 * (M2b ** 8) + 9272 * (M2b ** 9),
+
+    set(M1b, 0, (N1b `unchecked_right_shift` 15) /\
+        ((1 `unchecked_left_shift` 15) - 1), M1c),
+    set(M2b, 0, (N2b `unchecked_right_shift` 15) /\
+        ((1 `unchecked_left_shift` 15) - 1), M2c),
+
+    set(M1c, 1, N1b /\ ((1 `unchecked_left_shift` 15) - 1), M1d),
+    set(M2c, 1, N2b /\ ((1 `unchecked_left_shift` 15) - 1), M2d),
+
+    Seed1 = ((M1d ** 1) `unchecked_left_shift` 15) + (M2d ** 1),
+
+    Res = rfloat(Seed1)/rfloat((1 `unchecked_left_shift` 30) - 1).
+
+:- pred shift(vec::in, vec::out) is det.
+:- pragma inline(shift/2).
+
+shift(Vec0, Vec1) :-
+    Vec0 = vec(A, B, C, D, E, F, G, H, I, _),
+    Vec1 = vec(A, B, B, C, D, E, F, G, H, I).
+
+%-----------------------------------------------------------------------------%
+
+:- func (vec ** int) = int.
+:- mode ((in ** in) = out) is det.
+:- mode ((in ** in(bound(0))) = out) is det.
+:- mode ((in ** in(bound(1))) = out) is det.
+:- mode ((in ** in(bound(2))) = out) is det.
+:- mode ((in ** in(bound(3))) = out) is det.
+:- mode ((in ** in(bound(4))) = out) is det.
+:- mode ((in ** in(bound(5))) = out) is det.
+:- mode ((in ** in(bound(6))) = out) is det.
+:- mode ((in ** in(bound(7))) = out) is det.
+:- mode ((in ** in(bound(8))) = out) is det.
+:- mode ((in ** in(bound(9))) = out) is det.
+:- pragma inline((**)/2).
+
+( Vec ** Ind ) = Res :-
+    Vec = vec(A, B, C, D, E, F, G, H, I, J),
+    (
+        ( Ind = 0, Res0 = A
+        ; Ind = 1, Res0 = B
+        ; Ind = 2, Res0 = C
+        ; Ind = 3, Res0 = D
+        ; Ind = 4, Res0 = E
+        ; Ind = 5, Res0 = F
+        ; Ind = 6, Res0 = G
+        ; Ind = 7, Res0 = H
+        ; Ind = 8, Res0 = I
+        ; Ind = 9, Res0 = J
+        )
+    ->
+        Res = Res0
+    ;
+        error("**: out of range")
+    ).
+
+:- pred set(vec, int, int, vec).
+:- mode set(in, in, in, out) is det.
+:- mode set(in, in(bound(0)), in, out) is det.
+:- mode set(in, in(bound(1)), in, out) is det.
+:- mode set(in, in(bound(2)), in, out) is det.
+:- mode set(in, in(bound(3)), in, out) is det.
+:- mode set(in, in(bound(4)), in, out) is det.
+:- mode set(in, in(bound(5)), in, out) is det.
+:- mode set(in, in(bound(6)), in, out) is det.
+:- mode set(in, in(bound(7)), in, out) is det.
+:- mode set(in, in(bound(8)), in, out) is det.
+:- mode set(in, in(bound(9)), in, out) is det.
+:- pragma inline(rnd.set/4).
+
+set(Vec0, Ind, V, Vec) :-
+    Vec0 = vec(A, B, C, D, E, F, G, H, I, J),
+    (
+        ( Ind = 0, Vec1 = vec(V, B, C, D, E, F, G, H, I, J)
+        ; Ind = 1, Vec1 = vec(A, V, C, D, E, F, G, H, I, J)
+        ; Ind = 2, Vec1 = vec(A, B, V, D, E, F, G, H, I, J)
+        ; Ind = 3, Vec1 = vec(A, B, C, V, E, F, G, H, I, J)
+        ; Ind = 4, Vec1 = vec(A, B, C, D, V, F, G, H, I, J)
+        ; Ind = 5, Vec1 = vec(A, B, C, D, E, V, G, H, I, J)
+        ; Ind = 6, Vec1 = vec(A, B, C, D, E, F, V, H, I, J)
+        ; Ind = 7, Vec1 = vec(A, B, C, D, E, F, G, V, I, J)
+        ; Ind = 8, Vec1 = vec(A, B, C, D, E, F, G, H, V, J)
+        ; Ind = 9, Vec1 = vec(A, B, C, D, E, F, G, H, I, V)
+        )
+    ->
+        Vec = Vec1
+    ;
+        error("set: out of range")
+    ).
+
+%-----------------------------------------------------------------------------%
+
+:- func rfloat(int) = float.
+:- pragma foreign_proc("C",
+    rfloat(I::in) = (F::out),
+    [will_not_call_mercury, thread_safe, promise_pure],
+"
+    F = I;
+").
+
+:- func rint(float) = int.
+:- pragma foreign_proc("C",
+    rint(F::in) = (I::out),
+    [will_not_call_mercury, thread_safe, promise_pure],
+"
+    I = F;
+").
--------------------------------------------------------------------------
mercury-reviews mailing list
Post messages to:       mercury-reviews at csse.unimelb.edu.au
Administrative Queries: owner-mercury-reviews at csse.unimelb.edu.au
Subscriptions:          mercury-reviews-request at csse.unimelb.edu.au
--------------------------------------------------------------------------



More information about the reviews mailing list