[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