[pure-lang-svn] SF.net SVN: pure-lang: [423] pure/trunk
Status: Beta
Brought to you by:
agraef
Menu
▾
▴
[pure-lang-svn] SF.net SVN: pure-lang: [423] pure/trunk
|
From: <ag...@us...> - 2008-07-08 13:12:32
|
Revision: 423
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=423&view=rev
Author: agraef
Date: 2008-07-08 06:12:34 -0700 (Tue, 08 Jul 2008)
Log Message:
-----------
Add Mersenne twister to math.pure.
Modified Paths:
--------------
pure/trunk/ChangeLog
pure/trunk/lib/math.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
pure/trunk/test/test014.log
pure/trunk/test/test014.pure
Modified: pure/trunk/ChangeLog
===================================================================
--- pure/trunk/ChangeLog 2008-07-08 10:55:04 UTC (rev 422)
+++ pure/trunk/ChangeLog 2008-07-08 13:12:34 UTC (rev 423)
@@ -1,5 +1,8 @@
2008-07-08 Albert Graef <Dr....@t-...>
+ * runtime.cc/h, lib/math.pure: Add random number generator
+ (Mersenne twister). Suggested by Jiri Spitz.
+
* examples/avltree.pure: Added to examples.
* lib/math.pure: Moved abs, sgn, min, max, pred and succ from
Modified: pure/trunk/lib/math.pure
===================================================================
--- pure/trunk/lib/math.pure 2008-07-08 10:55:04 UTC (rev 422)
+++ pure/trunk/lib/math.pure 2008-07-08 13:12:34 UTC (rev 423)
@@ -22,6 +22,12 @@
def inf = 1.0e307 * 1.0e307; def nan = inf-inf;
+/* Random number generator. This uses the Mersenne twister, in order to avoid
+ bad generators present in some C libraries. Returns pseudo random ints in
+ the range -0x80000000..0x7fffffff. */
+
+extern int pure_random() = random, void pure_srandom(int) = srandom;
+
/* Rounding functions. */
extern double floor(double), double ceil(double);
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-07-08 10:55:04 UTC (rev 422)
+++ pure/trunk/runtime.cc 2008-07-08 13:12:34 UTC (rev 423)
@@ -1802,6 +1802,112 @@
return pure_tuplel(2, u, v);
}
+// This is the ``Mersenne Twister'' random number generator MT19937, which
+// generates pseudorandom integers uniformly distributed in 0..(2^32 - 1)
+// starting from any odd seed in 0..(2^32 - 1). This version is a recode
+// by Shawn Cokus (Co...@ma...) on March 8, 1998 of a version by
+// Takuji Nishimura (who had suggestions from Topher Cooper and Marc Rieffel in
+// July-August 1997).
+//
+// Effectiveness of the recoding (on Goedel2.math.washington.edu, a DEC Alpha
+// running OSF/1) using GCC -O3 as a compiler: before recoding: 51.6 sec. to
+// generate 300 million random numbers; after recoding: 24.0 sec. for the same
+// (i.e., 46.5% of original time), so speed is now about 12.5 million random
+// number generations per second on this machine.
+//
+// According to the URL <http://www.math.keio.ac.jp/~matumoto/emt.html>
+// (and paraphrasing a bit in places), the Mersenne Twister is ``designed
+// with consideration of the flaws of various existing generators,'' has
+// a period of 2^19937 - 1, gives a sequence that is 623-dimensionally
+// equidistributed, and ``has passed many stringent tests, including the
+// die-hard test of G. Marsaglia and the load test of P. Hellekalek and
+// S. Wegenkittl.'' It is efficient in memory usage (typically using 2506
+// to 5012 bytes of static data, depending on data type sizes, and the code
+// is quite short as well). It generates random numbers in batches of 624
+// at a time, so the caching and pipelining of modern systems is exploited.
+// It is also divide- and mod-free.
+//
+// This library is free software; you can redistribute it and/or modify it
+// under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation (either version 2 of the License or, at your
+// option, any later version). This library is distributed in the hope that
+// it will be useful, but WITHOUT ANY WARRANTY, without even the implied
+// warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See
+// the GNU Library General Public License for more details. You should have
+// received a copy of the GNU Library General Public License along with this
+// library; if not, write to the Free Software Foundation, Inc., 59 Temple
+// Place, Suite 330, Boston, MA 02111-1307, USA.
+//
+// The code as Shawn received it included the following notice:
+//
+// Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. When
+// you use this, send an e-mail to <mat...@ma...> with
+// an appropriate reference to your work.
+//
+// It would be nice to CC: <Co...@ma...> when you write.
+//
+
+// See http://www.math.keio.ac.jp/~matumoto/emt.html for the original sources.
+
+#define N (624)
+#define M (397)
+#define K (0x9908B0DFU)
+#define hiBit(u) ((u) & 0x80000000U)
+#define loBit(u) ((u) & 0x00000001U)
+#define loBits(u) ((u) & 0x7FFFFFFFU)
+#define mixBits(u, v) (hiBit(u)|loBits(v))
+
+// TLD?
+static uint32_t stateMT[N+1];
+static uint32_t *nextMT;
+static int leftMT = -1;
+
+void pure_srandom(uint32_t seed)
+{
+ // MT works best with odd seeds, so we enforce that here.
+ register uint32_t x = (seed | 1U) & 0xFFFFFFFFU, *s = stateMT;
+ register int j;
+
+ for (leftMT=0, *s++=x, j=N; --j; *s++ = (x*=69069U) & 0xFFFFFFFFU);
+}
+
+static uint32_t reloadMT(void)
+{
+ register uint32_t *p0=stateMT, *p2=stateMT+2, *pM=stateMT+M, s0, s1;
+ register int j;
+
+ if (leftMT < -1)
+ pure_srandom(4357U);
+
+ leftMT=N-1, nextMT=stateMT+1;
+
+ for (s0=stateMT[0], s1=stateMT[1], j=N-M+1; --j; s0=s1, s1=*p2++)
+ *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
+
+ for (pM=stateMT, j=M; --j; s0=s1, s1=*p2++)
+ *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
+
+ s1=stateMT[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
+ s1 ^= (s1 >> 11);
+ s1 ^= (s1 << 7) & 0x9D2C5680U;
+ s1 ^= (s1 << 15) & 0xEFC60000U;
+ return(s1 ^ (s1 >> 18));
+}
+
+uint32_t pure_random(void)
+{
+ uint32_t y;
+
+ if(--leftMT < 0)
+ return reloadMT();
+
+ y = *nextMT++;
+ y ^= (y >> 11);
+ y ^= (y << 7) & 0x9D2C5680U;
+ y ^= (y << 15) & 0xEFC60000U;
+ return (y ^ (y >> 18));
+}
+
extern "C"
pure_expr *bigint_neg(mpz_t x)
{
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-07-08 10:55:04 UTC (rev 422)
+++ pure/trunk/runtime.h 2008-07-08 13:12:34 UTC (rev 423)
@@ -420,6 +420,14 @@
pure_expr *pure_rational(double d);
+/* Random number generator. This uses the Mersenne twister, in order to avoid
+ bad generators present in some C libraries. pure_random returns a
+ pseudorandom 32 bit integer, pure_srandom sets the seed of the
+ generator. */
+
+uint32_t pure_random(void);
+void pure_srandom(uint32_t seed);
+
/* Construct a "byte string" from a string. The result is a raw pointer object
pointing to the converted string. The original string is copied (and, in
the case of pure_byte_cstring, converted to the system encoding). The
Modified: pure/trunk/test/test014.log
===================================================================
--- pure/trunk/test/test014.log 2008-07-08 10:55:04 UTC (rev 422)
+++ pure/trunk/test/test014.log 2008-07-08 13:12:34 UTC (rev 423)
@@ -1,3 +1,15 @@
+drop 97 (catmap (\i/*0:*/ -> [random] {
+ rule #0: i = [random]
+ state 0: #0
+ <var> state 1
+ state 1: #0
+}) (1..100)) when () = srandom 0 {
+ rule #0: () = srandom 0
+ state 0: #0
+ () state 1
+ state 1: #0
+} end;
+[1863734801,-639116898,52532575]
{
rule #0: q = 44%(-14)
state 0: #0
Modified: pure/trunk/test/test014.pure
===================================================================
--- pure/trunk/test/test014.pure 2008-07-08 10:55:04 UTC (rev 422)
+++ pure/trunk/test/test014.pure 2008-07-08 13:12:34 UTC (rev 423)
@@ -1,8 +1,12 @@
+using math;
+
+// Random number generator.
+
+drop 97 [random; i=1..100] when () = srandom 0 end;
+
// Some rational arithmetic tests, pilfered from Rob Hubbard's Q+Q manual.
-using math;
-
// Basic arithmetic and pow function.
let q = 44%(-14); q;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
