Added SRFI 25, 26, 27, 28, 30.
This commit is contained in:
mainzelm
parent
87846eef58
commit
024d938e39
|
|
@ -0,0 +1,227 @@
|
|||
/* 54-BIT (double) IMPLEMENTATION IN C OF THE "MRG32K3A" GENERATOR
|
||||
===============================================================
|
||||
|
||||
Sebastian.Egner@philips.com, Mar-2002, in ANSI-C and Scheme 48 0.57
|
||||
|
||||
This code is a C-implementation of Pierre L'Ecuyer's MRG32k3a generator.
|
||||
The code uses (double)-arithmetics, assuming that it covers the range
|
||||
{-2^53..2^53-1} exactly (!). The code of the generator is based on the
|
||||
L'Ecuyer's own implementation of the generator. Please refer to the
|
||||
file 'mrg32k3a.scm' for more information about the method.
|
||||
|
||||
The method provides the following functions via the C/Scheme
|
||||
interface of Scheme 48 0.57 to 'mrg32k3a-b.scm':
|
||||
|
||||
s48_value mrg32k3a_pack_state1(s48_value state);
|
||||
s48_value mrg32k3a_unpack_state1(s48_value state);
|
||||
s48_value mrg32k3a_random_range();
|
||||
s48_value mrg32k3a_random_integer(s48_value state, s48_value range);
|
||||
s48_value mrg32k3a_random_real(s48_value state);
|
||||
|
||||
As Scheme48 FIXNUMs cannot cover the range {0..m1-1}, we break up
|
||||
all values x in the state into x0+x1*w, where w = 2^16 = 65536.
|
||||
The procedures in Scheme correct for that.
|
||||
|
||||
compile this file with:
|
||||
gcc -c -I $SCHEME48 mrg32k3a-b.c
|
||||
|
||||
history of this file:
|
||||
SE, 18-Mar-2002: initial version
|
||||
SE, 22-Mar-2002: interface changed
|
||||
SE, 25-Mar-2002: tested with Scheme 48 0.57 in c/srfi-27
|
||||
SE, 27-Mar-2002: cleaned
|
||||
SE, 13-May-2002: bug found by Shiro Kawai removed
|
||||
*/
|
||||
|
||||
#include "scheme48.h" /* $SCHEME48/c/scheme48.h */
|
||||
#include <sys/time.h>
|
||||
|
||||
/* maximum value for random_integer: min(S48_MAX_FIXNUM_VALUE, m1) */
|
||||
#define m_max (((long)1 << 29) - 1)
|
||||
|
||||
/* The Generator
|
||||
=============
|
||||
*/
|
||||
|
||||
/* moduli of the components */
|
||||
#define m1 4294967087.0
|
||||
#define m2 4294944443.0
|
||||
|
||||
/* representation of the state in C */
|
||||
typedef struct {
|
||||
double
|
||||
x10, x11, x12,
|
||||
x20, x21, x22;
|
||||
} state_t;
|
||||
|
||||
/* recursion coefficients of the components */
|
||||
#define a12 1403580.0
|
||||
#define a13n 810728.0
|
||||
#define a21 527612.0
|
||||
#define a23n 1370589.0
|
||||
|
||||
/* normalization factor 1/(m1 + 1) */
|
||||
#define norm 2.328306549295728e-10
|
||||
|
||||
|
||||
/* the actual generator */
|
||||
|
||||
static double mrg32k3a(state_t *s) { /* (double), in {0..m1-1} */
|
||||
double x10, x20, y;
|
||||
long k10, k20;
|
||||
|
||||
/* #define debug 1 */
|
||||
|
||||
#if defined(debug)
|
||||
printf(
|
||||
"state = {%g %g %g %g %g %g};\n",
|
||||
s->x10, s->x11, s->x12,
|
||||
s->x20, s->x21, s->x22
|
||||
);
|
||||
#endif
|
||||
|
||||
/* component 1 */
|
||||
x10 = a12*(s->x11) - a13n*(s->x12);
|
||||
k10 = x10 / m1;
|
||||
x10 -= k10 * m1;
|
||||
if (x10 < 0.0)
|
||||
x10 += m1;
|
||||
s->x12 = s->x11;
|
||||
s->x11 = s->x10;
|
||||
s->x10 = x10;
|
||||
|
||||
/* component 2 */
|
||||
x20 = a21*(s->x20) - a23n*(s->x22);
|
||||
k20 = x20 / m2;
|
||||
x20 -= k20 * m2;
|
||||
if (x20 < 0.0)
|
||||
x20 += m2;
|
||||
s->x22 = s->x21;
|
||||
s->x21 = s->x20;
|
||||
s->x20 = x20;
|
||||
|
||||
/* combination of component */
|
||||
y = x10 - x20;
|
||||
if (y < 0.0)
|
||||
y += m1;
|
||||
return y;
|
||||
}
|
||||
|
||||
/* Exported Interface
|
||||
==================
|
||||
*/
|
||||
|
||||
s48_value mrg32k3a_pack_state1(s48_value state) {
|
||||
s48_value result;
|
||||
state_t s;
|
||||
|
||||
#define REF(i) (double)s48_extract_integer(S48_VECTOR_REF(state, (long)(i)))
|
||||
|
||||
/* copy the numbers from state into s */
|
||||
s.x10 = REF( 0) + 65536.0 * REF( 1);
|
||||
s.x11 = REF( 2) + 65536.0 * REF( 3);
|
||||
s.x12 = REF( 4) + 65536.0 * REF( 5);
|
||||
s.x20 = REF( 6) + 65536.0 * REF( 7);
|
||||
s.x21 = REF( 8) + 65536.0 * REF( 9);
|
||||
s.x22 = REF(10) + 65536.0 * REF(11);
|
||||
|
||||
#undef REF
|
||||
|
||||
/* box s into a Scheme object */
|
||||
result = S48_MAKE_VALUE(state_t);
|
||||
S48_SET_VALUE(result, state_t, s);
|
||||
return result;
|
||||
}
|
||||
|
||||
s48_value mrg32k3a_unpack_state1(s48_value state) {
|
||||
s48_value result;
|
||||
state_t s;
|
||||
|
||||
/* unbox s from the Scheme object */
|
||||
s = S48_EXTRACT_VALUE(state, state_t);
|
||||
|
||||
/* make and fill a Scheme vector with the numbers */
|
||||
result = s48_make_vector((long)12, S48_FALSE);
|
||||
|
||||
#define SET(i, x) { \
|
||||
long x1 = (long)((x) / 65536.0); \
|
||||
long x0 = (long)((x) - 65536.0 * (double)x1); \
|
||||
S48_VECTOR_SET(result, (long)(i+0), s48_enter_integer(x0)); \
|
||||
S48_VECTOR_SET(result, (long)(i+1), s48_enter_integer(x1)); }
|
||||
|
||||
SET( 0, s.x10);
|
||||
SET( 2, s.x11);
|
||||
SET( 4, s.x12);
|
||||
SET( 6, s.x20);
|
||||
SET( 8, s.x21);
|
||||
SET(10, s.x22);
|
||||
|
||||
#undef SET
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
s48_value mrg32k3a_random_range(void) {
|
||||
return s48_enter_fixnum(m_max);
|
||||
}
|
||||
|
||||
s48_value mrg32k3a_random_integer(s48_value state, s48_value range) {
|
||||
long result;
|
||||
state_t s;
|
||||
long n;
|
||||
double x, q, qn, xq;
|
||||
|
||||
s = S48_EXTRACT_VALUE(state, state_t);
|
||||
n = s48_extract_integer(range);
|
||||
if (!( ((long)1 <= n) && (n <= m_max) ))
|
||||
s48_raise_range_error(n, (long)1, m_max);
|
||||
|
||||
/* generate result in {0..n-1} using the rejection method */
|
||||
q = (double)( (unsigned long)(m1 / (double)n) );
|
||||
qn = q * n;
|
||||
do {
|
||||
x = mrg32k3a(&s);
|
||||
} while (x >= qn);
|
||||
xq = x / q;
|
||||
|
||||
/* check the range */
|
||||
if (!( (0.0 <= xq) && (xq < (double)m_max) ))
|
||||
s48_raise_range_error((long)xq, (long)0, m_max);
|
||||
|
||||
/* return result */
|
||||
result = (long)xq;
|
||||
S48_SET_VALUE(state, state_t, s);
|
||||
return s48_enter_fixnum(result);
|
||||
}
|
||||
|
||||
s48_value mrg32k3a_random_real(s48_value state) {
|
||||
state_t s;
|
||||
double x;
|
||||
|
||||
s = S48_EXTRACT_VALUE(state, state_t);
|
||||
x = (mrg32k3a(&s) + 1.0) * norm;
|
||||
S48_SET_VALUE(state, state_t, s);
|
||||
return s48_enter_double(x);
|
||||
}
|
||||
|
||||
/* Kludge for scsh */
|
||||
static s48_value current_time(void){
|
||||
struct timeval tv;
|
||||
gettimeofday(&tv, NULL);
|
||||
return s48_enter_integer(tv.tv_sec);
|
||||
}
|
||||
|
||||
|
||||
/* Exporting the C values to Scheme
|
||||
================================
|
||||
*/
|
||||
|
||||
void s48_init_srfi_27(void) {
|
||||
S48_EXPORT_FUNCTION(mrg32k3a_pack_state1);
|
||||
S48_EXPORT_FUNCTION(mrg32k3a_unpack_state1);
|
||||
S48_EXPORT_FUNCTION(mrg32k3a_random_range);
|
||||
S48_EXPORT_FUNCTION(mrg32k3a_random_integer);
|
||||
S48_EXPORT_FUNCTION(mrg32k3a_random_real);
|
||||
S48_EXPORT_FUNCTION(current_time);
|
||||
}
|
||||
|
||||
|
|
@ -532,3 +532,15 @@
|
|||
|
||||
))
|
||||
|
||||
(define-interface srfi-27-interface
|
||||
(export random-integer
|
||||
random-real
|
||||
default-random-source
|
||||
make-random-source
|
||||
random-source?
|
||||
random-source-state-ref
|
||||
random-source-state-set!
|
||||
random-source-randomize!
|
||||
random-source-pseudo-randomize!
|
||||
random-source-make-integers
|
||||
random-source-make-reals))
|
||||
|
|
@ -728,7 +728,8 @@
|
|||
(begin
|
||||
(define available-srfis
|
||||
'(srfi-1 srfi-2 srfi-5 srfi-6 srfi-7 srfi-8 srfi-9
|
||||
srfi-11 srfi-13 srfi-14 srfi-16 srfi-17 srfi-23))
|
||||
srfi-11 srfi-13 srfi-14 srfi-16 srfi-17 srfi-23
|
||||
srfi-25 srfi-26 srfi-27 srfi-28 srfi-30))
|
||||
|
||||
; Some SRFI's redefine Scheme variables.
|
||||
(define shadowed
|
||||
|
|
@ -807,11 +808,50 @@
|
|||
; SRFI-19 - implementation is specific to MzScheme
|
||||
; SRFI-20 - withdrawn
|
||||
; SRFI-21 - no implementation given
|
||||
; SRFI-22 - not final yet
|
||||
; SRFI-22 - no implementation given
|
||||
|
||||
(define-structure srfi-23 (export error)
|
||||
(open (subset signals (error))))
|
||||
|
||||
; SRFI-24 - withdrawn
|
||||
|
||||
(define-structure srfi-25 (export
|
||||
array? make-array shape array
|
||||
array-rank
|
||||
array-start array-end
|
||||
array-ref array-set! share-array)
|
||||
(open scheme
|
||||
srfi-23
|
||||
srfi-9)
|
||||
(files (srfi srfi-25)))
|
||||
|
||||
|
||||
(define-structure srfi-26 (export ((cut cute) :syntax))
|
||||
(open scheme)
|
||||
(files (srfi srfi-26)))
|
||||
|
||||
(define-structure srfi-27 srfi-27-interface
|
||||
(open
|
||||
scheme-level-1
|
||||
floatnums
|
||||
external-calls
|
||||
(subset srfi-9 (define-record-type))
|
||||
(subset srfi-23 (error)))
|
||||
;; scsh doesn't have S48's posix subsystem yet:
|
||||
; (subset posix-time (current-time))
|
||||
; (subset posix (time-seconds)))
|
||||
(files (srfi srfi-27)))
|
||||
|
||||
(define-structure srfi-28 (export format)
|
||||
(open scheme
|
||||
srfi-23
|
||||
srfi-6)
|
||||
(files (srfi srfi-28)))
|
||||
|
||||
; SRFI-29 - requires access to the current locale
|
||||
|
||||
; SRFI-30 - scheme/rts/read.scm contains the reader for #|...|# comments
|
||||
|
||||
; ... end of package definitions.
|
||||
|
||||
; Temporary compatibility stuff
|
||||
|
|
@ -897,7 +937,7 @@
|
|||
; SRFI packages
|
||||
srfi-1 srfi-2 srfi-5 srfi-6 srfi-7 srfi-8 srfi-9
|
||||
srfi-11 srfi-13 srfi-14 srfi-16 srfi-17
|
||||
srfi-23
|
||||
srfi-23 srfi-25 srfi-26 srfi-27 srfi-28
|
||||
)
|
||||
:structure)
|
||||
((define-signature define-package) :syntax)))
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,96 @@
|
|||
; REFERENCE IMPLEMENTATION FOR SRFI-26 "CUT"
|
||||
; ==========================================
|
||||
;
|
||||
; Sebastian.Egner@philips.com, 5-Jun-2002.
|
||||
; adapted from the posting by Al Petrofsky <al@petrofsky.org>
|
||||
; placed in the public domain
|
||||
;
|
||||
; The code to handle the variable argument case was originally
|
||||
; proposed by Michael Sperber and has been adapted to the new
|
||||
; syntax of the macro using an explicit rest-slot symbol. The
|
||||
; code to evaluate the non-slots for cute has been proposed by
|
||||
; Dale Jordan. The code to allow a slot for the procedure position
|
||||
; and to process the macro using an internal macro is based on
|
||||
; a suggestion by Al Petrofsky. The code found below is, with
|
||||
; exception of this header and some changes in variable names,
|
||||
; entirely written by Al Petrofsky.
|
||||
;
|
||||
; compliance:
|
||||
; Scheme R5RS (including macros).
|
||||
;
|
||||
; loading this file into Scheme 48 0.57:
|
||||
; ,load cut.scm
|
||||
;
|
||||
; history of this file:
|
||||
; SE, 6-Feb-2002: initial version as 'curry' with ". <>" notation
|
||||
; SE, 14-Feb-2002: revised for <...>
|
||||
; SE, 27-Feb-2002: revised for 'cut'
|
||||
; SE, 03-Jun-2002: revised for proc-slot, cute
|
||||
; SE, 04-Jun-2002: rewritten with internal transformer (no "loop" pattern)
|
||||
; SE, 05-Jun-2002: replace my code by Al's; substituted "constant" etc.
|
||||
; to match the convention in the SRFI-document
|
||||
|
||||
; (srfi-26-internal-cut slot-names combination . se)
|
||||
; transformer used internally
|
||||
; slot-names : the internal names of the slots
|
||||
; combination : procedure being specialized, followed by its arguments
|
||||
; se : slots-or-exprs, the qualifiers of the macro
|
||||
|
||||
(define-syntax srfi-26-internal-cut
|
||||
(syntax-rules (<> <...>)
|
||||
|
||||
;; construct fixed- or variable-arity procedure:
|
||||
;; (begin proc) throws an error if proc is not an <expression>
|
||||
((srfi-26-internal-cut (slot-name ...) (proc arg ...))
|
||||
(lambda (slot-name ...) ((begin proc) arg ...)))
|
||||
((srfi-26-internal-cut (slot-name ...) (proc arg ...) <...>)
|
||||
(lambda (slot-name ... . rest-slot) (apply proc arg ... rest-slot)))
|
||||
|
||||
;; process one slot-or-expr
|
||||
((srfi-26-internal-cut (slot-name ...) (position ...) <> . se)
|
||||
(srfi-26-internal-cut (slot-name ... x) (position ... x) . se))
|
||||
((srfi-26-internal-cut (slot-name ...) (position ...) nse . se)
|
||||
(srfi-26-internal-cut (slot-name ...) (position ... nse) . se))))
|
||||
|
||||
; (srfi-26-internal-cute slot-names nse-bindings combination . se)
|
||||
; transformer used internally
|
||||
; slot-names : the internal names of the slots
|
||||
; nse-bindings : let-style bindings for the non-slot expressions.
|
||||
; combination : procedure being specialized, followed by its arguments
|
||||
; se : slots-or-exprs, the qualifiers of the macro
|
||||
|
||||
(define-syntax srfi-26-internal-cute
|
||||
(syntax-rules (<> <...>)
|
||||
|
||||
;; If there are no slot-or-exprs to process, then:
|
||||
;; construct a fixed-arity procedure,
|
||||
((srfi-26-internal-cute
|
||||
(slot-name ...) nse-bindings (proc arg ...))
|
||||
(let nse-bindings (lambda (slot-name ...) (proc arg ...))))
|
||||
;; or a variable-arity procedure
|
||||
((srfi-26-internal-cute
|
||||
(slot-name ...) nse-bindings (proc arg ...) <...>)
|
||||
(let nse-bindings (lambda (slot-name ... . x) (apply proc arg ... x))))
|
||||
|
||||
;; otherwise, process one slot:
|
||||
((srfi-26-internal-cute
|
||||
(slot-name ...) nse-bindings (position ...) <> . se)
|
||||
(srfi-26-internal-cute
|
||||
(slot-name ... x) nse-bindings (position ... x) . se))
|
||||
;; or one non-slot expression
|
||||
((srfi-26-internal-cute
|
||||
slot-names nse-bindings (position ...) nse . se)
|
||||
(srfi-26-internal-cute
|
||||
slot-names ((x nse) . nse-bindings) (position ... x) . se))))
|
||||
|
||||
; exported syntax
|
||||
|
||||
(define-syntax cut
|
||||
(syntax-rules ()
|
||||
((cut . slots-or-exprs)
|
||||
(srfi-26-internal-cut () () . slots-or-exprs))))
|
||||
|
||||
(define-syntax cute
|
||||
(syntax-rules ()
|
||||
((cute . slots-or-exprs)
|
||||
(srfi-26-internal-cute () () () . slots-or-exprs))))
|
||||
|
|
@ -0,0 +1,569 @@
|
|||
; MODULE DEFINITION FOR SRFI-27, C/SCHEME-IMPLEMENTATION
|
||||
; ======================================================
|
||||
;
|
||||
; Sebastian.Egner@philips.com, Mar-2002, in Scheme 48 0.57
|
||||
;
|
||||
; This file contains the top-level definition for the C-code
|
||||
; implementation of SRFI-27 for the Scheme 48 0.57 system.
|
||||
;
|
||||
; 1. The core generator is implemented in 'mrg32k3a-b.c'.
|
||||
; 2. The generic parts of the interface are in 'mrg32k3a.scm'.
|
||||
; 3. The non-generic parts (record type, time, error, C-bindings) are here.
|
||||
;
|
||||
; creating the module:
|
||||
; copy mrg32k3a-b.c into $SCHEME48/c/srfi-27/mrg32k3a-b.c
|
||||
; edit $SCHEME48/Makefile.in
|
||||
; add c/srfi-27/mrg32k3a-b.o to EXTERNAL_OBJECTS
|
||||
; add mrg32k3a_init to EXTERNAL_INITIALIZERS
|
||||
; add the make line c/srfi-27/mrg32k3a-b.o: c/scheme48.h
|
||||
; cd $SCHEME48
|
||||
; make clean
|
||||
; configure
|
||||
; make
|
||||
; cd $SRFI27
|
||||
; ,config ,load srfi-27-b.scm
|
||||
;
|
||||
; loading the module, once created:
|
||||
; ,open srfi-27
|
||||
;
|
||||
; history of this file:
|
||||
; SE, 22-Mar-2002: initial version
|
||||
; SE, 25-Mar-2002: initial version
|
||||
; MG, September 2002: merged in mrg32k2a.scm, move package definitons to
|
||||
; more-packages.scm
|
||||
|
||||
(define-record-type :random-source
|
||||
(:random-source-make
|
||||
state-ref
|
||||
state-set!
|
||||
randomize!
|
||||
pseudo-randomize!
|
||||
make-integers
|
||||
make-reals)
|
||||
:random-source?
|
||||
(state-ref :random-source-state-ref)
|
||||
(state-set! :random-source-state-set!)
|
||||
(randomize! :random-source-randomize!)
|
||||
(pseudo-randomize! :random-source-pseudo-randomize!)
|
||||
(make-integers :random-source-make-integers)
|
||||
(make-reals :random-source-make-reals))
|
||||
|
||||
; We have neither scsh nor posix...
|
||||
; (define (:random-source-current-time)
|
||||
; (time-seconds (current-time)))
|
||||
(import-lambda-definition :random-source-current-time () "current_time")
|
||||
|
||||
; interface to core generator
|
||||
|
||||
(import-lambda-definition mrg32k3a-pack-state1 (state))
|
||||
(import-lambda-definition mrg32k3a-unpack-state1 (state))
|
||||
(import-lambda-definition mrg32k3a-random-range ())
|
||||
(import-lambda-definition mrg32k3a-random-integer (state range))
|
||||
(import-lambda-definition mrg32k3a-random-real (state))
|
||||
|
||||
(define (mrg32k3a-pack-state state)
|
||||
(mrg32k3a-pack-state1
|
||||
(list->vector
|
||||
(apply append
|
||||
(map (lambda (x)
|
||||
(list (modulo x 65536)
|
||||
(quotient x 65536)))
|
||||
(vector->list state))))))
|
||||
|
||||
(define (mrg32k3a-unpack-state state)
|
||||
(let ((s (mrg32k3a-unpack-state1 state)) (w 65536))
|
||||
(vector
|
||||
(+ (vector-ref s 0) (* (vector-ref s 1) w))
|
||||
(+ (vector-ref s 2) (* (vector-ref s 3) w))
|
||||
(+ (vector-ref s 4) (* (vector-ref s 5) w))
|
||||
(+ (vector-ref s 6) (* (vector-ref s 7) w))
|
||||
(+ (vector-ref s 8) (* (vector-ref s 9) w))
|
||||
(+ (vector-ref s 10) (* (vector-ref s 11) w)))))
|
||||
; GENERIC PART OF MRG32k3a-GENERATOR FOR SRFI-27
|
||||
; ==============================================
|
||||
;
|
||||
; Sebastian.Egner@philips.com, 2002.
|
||||
;
|
||||
; This is the generic R5RS-part of the implementation of the MRG32k3a
|
||||
; generator to be used in SRFI-27. It is based on a separate implementation
|
||||
; of the core generator (presumably in native code) and on code to
|
||||
; provide essential functionality not available in R5RS (see below).
|
||||
;
|
||||
; compliance:
|
||||
; Scheme R5RS with integer covering at least {-2^53..2^53-1}.
|
||||
; In addition,
|
||||
; SRFI-23: error
|
||||
;
|
||||
; history of this file:
|
||||
; SE, 22-Mar-2002: refactored from earlier versions
|
||||
; SE, 25-Mar-2002: pack/unpack need not allocate
|
||||
; SE, 27-Mar-2002: changed interface to core generator
|
||||
; SE, 10-Apr-2002: updated spec of mrg32k3a-random-integer
|
||||
|
||||
; Generator
|
||||
; =========
|
||||
;
|
||||
; Pierre L'Ecuyer's MRG32k3a generator is a Combined Multiple Recursive
|
||||
; Generator. It produces the sequence {(x[1,n] - x[2,n]) mod m1 : n}
|
||||
; defined by the two recursive generators
|
||||
;
|
||||
; x[1,n] = ( a12 x[1,n-2] + a13 x[1,n-3]) mod m1,
|
||||
; x[2,n] = (a21 x[2,n-1] + a23 x[2,n-3]) mod m2,
|
||||
;
|
||||
; where the constants are
|
||||
; m1 = 4294967087 = 2^32 - 209 modulus of 1st component
|
||||
; m2 = 4294944443 = 2^32 - 22853 modulus of 2nd component
|
||||
; a12 = 1403580 recursion coefficients
|
||||
; a13 = -810728
|
||||
; a21 = 527612
|
||||
; a23 = -1370589
|
||||
;
|
||||
; The generator passes all tests of G. Marsaglia's Diehard testsuite.
|
||||
; Its period is (m1^3 - 1)(m2^3 - 1)/2 which is nearly 2^191.
|
||||
; L'Ecuyer reports: "This generator is well-behaved in all dimensions
|
||||
; up to at least 45: ..." [with respect to the spectral test, SE].
|
||||
;
|
||||
; The period is maximal for all values of the seed as long as the
|
||||
; state of both recursive generators is not entirely zero.
|
||||
;
|
||||
; As the successor state is a linear combination of previous
|
||||
; states, it is possible to advance the generator by more than one
|
||||
; iteration by applying a linear transformation. The following
|
||||
; publication provides detailed information on how to do that:
|
||||
;
|
||||
; [1] P. L'Ecuyer, R. Simard, E. J. Chen, W. D. Kelton:
|
||||
; An Object-Oriented Random-Number Package With Many Long
|
||||
; Streams and Substreams. 2001.
|
||||
; To appear in Operations Research.
|
||||
;
|
||||
; Arithmetics
|
||||
; ===========
|
||||
;
|
||||
; The MRG32k3a generator produces values in {0..2^32-209-1}. All
|
||||
; subexpressions of the actual generator fit into {-2^53..2^53-1}.
|
||||
; The code below assumes that Scheme's "integer" covers this range.
|
||||
; In addition, it is assumed that floating point literals can be
|
||||
; read and there is some arithmetics with inexact numbers.
|
||||
;
|
||||
; However, for advancing the state of the generator by more than
|
||||
; one step at a time, the full range {0..2^32-209-1} is needed.
|
||||
|
||||
|
||||
; Required: Backbone Generator
|
||||
; ============================
|
||||
;
|
||||
; At this point in the code, the following procedures are assumed
|
||||
; to be defined to execute the core generator:
|
||||
;
|
||||
; (mrg32k3a-pack-state unpacked-state) -> packed-state
|
||||
; (mrg32k3a-unpack-state packed-state) -> unpacked-state
|
||||
; pack/unpack a state of the generator. The core generator works
|
||||
; on packed states, passed as an explicit argument, only. This
|
||||
; allows native code implementations to store their state in a
|
||||
; suitable form. Unpacked states are #(x10 x11 x12 x20 x21 x22)
|
||||
; with integer x_ij. Pack/unpack need not allocate new objects
|
||||
; in case packed and unpacked states are identical.
|
||||
;
|
||||
; (mrg32k3a-random-range) -> m-max
|
||||
; (mrg32k3a-random-integer packed-state range) -> x in {0..range-1}
|
||||
; advance the state of the generator and return the next random
|
||||
; range-limited integer.
|
||||
; Note that the state is not necessarily advanced by just one
|
||||
; step because we use the rejection method to avoid any problems
|
||||
; with distribution anomalies.
|
||||
; The range argument must be an exact integer in {1..m-max}.
|
||||
; It can be assumed that range is a fixnum if the Scheme system
|
||||
; has such a number representation.
|
||||
;
|
||||
; (mrg32k3a-random-real packed-state) -> x in (0,1)
|
||||
; advance the state of the generator and return the next random
|
||||
; real number between zero and one (both excluded). The type of
|
||||
; the result should be a flonum if possible.
|
||||
|
||||
; Required: Record Data Type
|
||||
; ==========================
|
||||
;
|
||||
; At this point in the code, the following procedures are assumed
|
||||
; to be defined to create and access a new record data type:
|
||||
;
|
||||
; (:random-source-make a0 a1 a2 a3 a4 a5) -> s
|
||||
; constructs a new random source object s consisting of the
|
||||
; objects a0 .. a5 in this order.
|
||||
;
|
||||
; (:random-source? obj) -> bool
|
||||
; tests if a Scheme object is a :random-source.
|
||||
;
|
||||
; (:random-source-state-ref s) -> a0
|
||||
; (:random-source-state-set! s) -> a1
|
||||
; (:random-source-randomize! s) -> a2
|
||||
; (:random-source-pseudo-randomize! s) -> a3
|
||||
; (:random-source-make-integers s) -> a4
|
||||
; (:random-source-make-reals s) -> a5
|
||||
; retrieve the values in the fields of the object s.
|
||||
|
||||
; Required: Current Time as an Integer
|
||||
; ====================================
|
||||
;
|
||||
; At this point in the code, the following procedure is assumed
|
||||
; to be defined to obtain a value that is likely to be different
|
||||
; for each invokation of the Scheme system:
|
||||
;
|
||||
; (:random-source-current-time) -> x
|
||||
; an integer that depends on the system clock. It is desired
|
||||
; that the integer changes as fast as possible.
|
||||
|
||||
|
||||
; Accessing the State
|
||||
; ===================
|
||||
|
||||
(define (mrg32k3a-state-ref packed-state)
|
||||
(cons 'lecuyer-mrg32k3a
|
||||
(vector->list (mrg32k3a-unpack-state packed-state))))
|
||||
|
||||
(define (mrg32k3a-state-set external-state)
|
||||
|
||||
(define (check-value x m)
|
||||
(if (and (integer? x)
|
||||
(exact? x)
|
||||
(<= 0 x (- m 1)))
|
||||
#t
|
||||
(error "illegal value" x)))
|
||||
|
||||
(if (and (list? external-state)
|
||||
(= (length external-state) 7)
|
||||
(eq? (car external-state) 'lecuyer-mrg32k3a))
|
||||
(let ((s (cdr external-state)))
|
||||
(check-value (list-ref s 0) mrg32k3a-m1)
|
||||
(check-value (list-ref s 1) mrg32k3a-m1)
|
||||
(check-value (list-ref s 2) mrg32k3a-m1)
|
||||
(check-value (list-ref s 3) mrg32k3a-m2)
|
||||
(check-value (list-ref s 4) mrg32k3a-m2)
|
||||
(check-value (list-ref s 5) mrg32k3a-m2)
|
||||
(if (or (zero? (+ (list-ref s 0) (list-ref s 1) (list-ref s 2)))
|
||||
(zero? (+ (list-ref s 3) (list-ref s 4) (list-ref s 5))))
|
||||
(error "illegal degenerate state" external-state))
|
||||
(mrg32k3a-pack-state (list->vector s)))
|
||||
(error "malformed state" external-state)))
|
||||
|
||||
|
||||
; Pseudo-Randomization
|
||||
; ====================
|
||||
;
|
||||
; Reference [1] above shows how to obtain many long streams and
|
||||
; substream from the backbone generator.
|
||||
;
|
||||
; The idea is that the generator is a linear operation on the state.
|
||||
; Hence, we can express this operation as a 3x3-matrix acting on the
|
||||
; three most recent states. Raising the matrix to the k-th power, we
|
||||
; obtain the operation to advance the state by k steps at once. The
|
||||
; virtual streams and substreams are now simply parts of the entire
|
||||
; periodic sequence (which has period around 2^191).
|
||||
;
|
||||
; For the implementation it is necessary to compute with matrices in
|
||||
; the ring (Z/(m1*m1)*Z)^(3x3). By the Chinese-Remainder Theorem, this
|
||||
; is isomorphic to ((Z/m1*Z) x (Z/m2*Z))^(3x3). We represent such a pair
|
||||
; of matrices
|
||||
; [ [[x00 x01 x02],
|
||||
; [x10 x11 x12],
|
||||
; [x20 x21 x22]], mod m1
|
||||
; [[y00 y01 y02],
|
||||
; [y10 y11 y12],
|
||||
; [y20 y21 y22]] mod m2]
|
||||
; as a vector of length 18 of the integers as writen above:
|
||||
; #(x00 x01 x02 x10 x11 x12 x20 x21 x22
|
||||
; y00 y01 y02 y10 y11 y12 y20 y21 y22)
|
||||
;
|
||||
; As the implementation should only use the range {-2^53..2^53-1}, the
|
||||
; fundamental operation (x*y) mod m, where x, y, m are nearly 2^32,
|
||||
; is computed by breaking up x and y as x = x1*w + x0 and y = y1*w + y0
|
||||
; where w = 2^16. In this case, all operations fit the range because
|
||||
; w^2 mod m is a small number. If proper multiprecision integers are
|
||||
; available this is not necessary, but pseudo-randomize! is an expected
|
||||
; to be called only occasionally so we do not provide this implementation.
|
||||
|
||||
(define mrg32k3a-m1 4294967087) ; modulus of component 1
|
||||
(define mrg32k3a-m2 4294944443) ; modulus of component 2
|
||||
|
||||
(define mrg32k3a-initial-state ; 0 3 6 9 12 15 of A^16, see below
|
||||
'#( 1062452522
|
||||
2961816100
|
||||
342112271
|
||||
2854655037
|
||||
3321940838
|
||||
3542344109))
|
||||
|
||||
(define mrg32k3a-generators #f) ; computed when needed
|
||||
|
||||
(define (mrg32k3a-pseudo-randomize-state i j)
|
||||
|
||||
(define (product A B) ; A*B in ((Z/m1*Z) x (Z/m2*Z))^(3x3)
|
||||
|
||||
(define w 65536) ; wordsize to split {0..2^32-1}
|
||||
(define w-sqr1 209) ; w^2 mod m1
|
||||
(define w-sqr2 22853) ; w^2 mod m2
|
||||
|
||||
(define (lc i0 i1 i2 j0 j1 j2 m w-sqr) ; linear combination
|
||||
(let ((a0h (quotient (vector-ref A i0) w))
|
||||
(a0l (modulo (vector-ref A i0) w))
|
||||
(a1h (quotient (vector-ref A i1) w))
|
||||
(a1l (modulo (vector-ref A i1) w))
|
||||
(a2h (quotient (vector-ref A i2) w))
|
||||
(a2l (modulo (vector-ref A i2) w))
|
||||
(b0h (quotient (vector-ref B j0) w))
|
||||
(b0l (modulo (vector-ref B j0) w))
|
||||
(b1h (quotient (vector-ref B j1) w))
|
||||
(b1l (modulo (vector-ref B j1) w))
|
||||
(b2h (quotient (vector-ref B j2) w))
|
||||
(b2l (modulo (vector-ref B j2) w)))
|
||||
(modulo
|
||||
(+ (* (+ (* a0h b0h)
|
||||
(* a1h b1h)
|
||||
(* a2h b2h))
|
||||
w-sqr)
|
||||
(* (+ (* a0h b0l)
|
||||
(* a0l b0h)
|
||||
(* a1h b1l)
|
||||
(* a1l b1h)
|
||||
(* a2h b2l)
|
||||
(* a2l b2h))
|
||||
w)
|
||||
(* a0l b0l)
|
||||
(* a1l b1l)
|
||||
(* a2l b2l))
|
||||
m)))
|
||||
|
||||
(vector
|
||||
(lc 0 1 2 0 3 6 mrg32k3a-m1 w-sqr1) ; (A*B)_00 mod m1
|
||||
(lc 0 1 2 1 4 7 mrg32k3a-m1 w-sqr1) ; (A*B)_01
|
||||
(lc 0 1 2 2 5 8 mrg32k3a-m1 w-sqr1)
|
||||
(lc 3 4 5 0 3 6 mrg32k3a-m1 w-sqr1) ; (A*B)_10
|
||||
(lc 3 4 5 1 4 7 mrg32k3a-m1 w-sqr1)
|
||||
(lc 3 4 5 2 5 8 mrg32k3a-m1 w-sqr1)
|
||||
(lc 6 7 8 0 3 6 mrg32k3a-m1 w-sqr1)
|
||||
(lc 6 7 8 1 4 7 mrg32k3a-m1 w-sqr1)
|
||||
(lc 6 7 8 2 5 8 mrg32k3a-m1 w-sqr1)
|
||||
(lc 9 10 11 9 12 15 mrg32k3a-m2 w-sqr2) ; (A*B)_00 mod m2
|
||||
(lc 9 10 11 10 13 16 mrg32k3a-m2 w-sqr2)
|
||||
(lc 9 10 11 11 14 17 mrg32k3a-m2 w-sqr2)
|
||||
(lc 12 13 14 9 12 15 mrg32k3a-m2 w-sqr2)
|
||||
(lc 12 13 14 10 13 16 mrg32k3a-m2 w-sqr2)
|
||||
(lc 12 13 14 11 14 17 mrg32k3a-m2 w-sqr2)
|
||||
(lc 15 16 17 9 12 15 mrg32k3a-m2 w-sqr2)
|
||||
(lc 15 16 17 10 13 16 mrg32k3a-m2 w-sqr2)
|
||||
(lc 15 16 17 11 14 17 mrg32k3a-m2 w-sqr2)))
|
||||
|
||||
(define (power A e) ; A^e
|
||||
(cond
|
||||
((zero? e)
|
||||
'#(1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1))
|
||||
((= e 1)
|
||||
A)
|
||||
((even? e)
|
||||
(power (product A A) (quotient e 2)))
|
||||
(else
|
||||
(product (power A (- e 1)) A))))
|
||||
|
||||
(define (power-power A b) ; A^(2^b)
|
||||
(if (zero? b)
|
||||
A
|
||||
(power-power (product A A) (- b 1))))
|
||||
|
||||
(define A ; the MRG32k3a recursion
|
||||
'#( 0 1403580 4294156359
|
||||
1 0 0
|
||||
0 1 0
|
||||
527612 0 4293573854
|
||||
1 0 0
|
||||
0 1 0))
|
||||
|
||||
; check arguments
|
||||
(if (not (and (integer? i)
|
||||
(exact? i)
|
||||
(integer? j)
|
||||
(exact? j)))
|
||||
(error "i j must be exact integer" i j))
|
||||
|
||||
; precompute A^(2^127) and A^(2^76) only once
|
||||
|
||||
(if (not mrg32k3a-generators)
|
||||
(set! mrg32k3a-generators
|
||||
(list (power-power A 127)
|
||||
(power-power A 76)
|
||||
(power A 16))))
|
||||
|
||||
; compute M = A^(16 + i*2^127 + j*2^76)
|
||||
(let ((M (product
|
||||
(list-ref mrg32k3a-generators 2)
|
||||
(product
|
||||
(power (list-ref mrg32k3a-generators 0)
|
||||
(modulo i (expt 2 28)))
|
||||
(power (list-ref mrg32k3a-generators 1)
|
||||
(modulo j (expt 2 28)))))))
|
||||
(mrg32k3a-pack-state
|
||||
(vector
|
||||
(vector-ref M 0)
|
||||
(vector-ref M 3)
|
||||
(vector-ref M 6)
|
||||
(vector-ref M 9)
|
||||
(vector-ref M 12)
|
||||
(vector-ref M 15)))))
|
||||
|
||||
; True Randomization
|
||||
; ==================
|
||||
;
|
||||
; The value obtained from the system time is feed into a very
|
||||
; simple pseudo random number generator. This in turn is used
|
||||
; to obtain numbers to randomize the state of the MRG32k3a
|
||||
; generator, avoiding period degeneration.
|
||||
|
||||
(define (mrg32k3a-randomize-state state)
|
||||
|
||||
; G. Marsaglia's simple 16-bit generator with carry
|
||||
(define m 65536)
|
||||
(define x (modulo (:random-source-current-time) m))
|
||||
(define (random-m)
|
||||
(let ((y (modulo x m)))
|
||||
(set! x (+ (* 30903 y) (quotient x m)))
|
||||
y))
|
||||
(define (random n) ; m < n < m^2
|
||||
(modulo (+ (* (random-m) m) (random-m)) n))
|
||||
|
||||
; modify the state
|
||||
(let ((m1 mrg32k3a-m1)
|
||||
(m2 mrg32k3a-m2)
|
||||
(s (mrg32k3a-unpack-state state)))
|
||||
(mrg32k3a-pack-state
|
||||
(vector
|
||||
(+ 1 (modulo (+ (vector-ref s 0) (random (- m1 1))) (- m1 1)))
|
||||
(modulo (+ (vector-ref s 1) (random m1)) m1)
|
||||
(modulo (+ (vector-ref s 2) (random m1)) m1)
|
||||
(+ 1 (modulo (+ (vector-ref s 3) (random (- m2 1))) (- m2 1)))
|
||||
(modulo (+ (vector-ref s 4) (random m2)) m2)
|
||||
(modulo (+ (vector-ref s 5) (random m2)) m2)))))
|
||||
|
||||
|
||||
; Large Integers
|
||||
; ==============
|
||||
;
|
||||
; To produce large integer random deviates, for n > m-max, we first
|
||||
; construct large random numbers in the range {0..m-max^k-1} for some
|
||||
; k such that m-max^k >= n and then use the rejection method to choose
|
||||
; uniformly from the range {0..n-1}.
|
||||
|
||||
(define mrg32k3a-m-max
|
||||
(mrg32k3a-random-range))
|
||||
|
||||
(define (mrg32k3a-random-power state k) ; n = m-max^k, k >= 1
|
||||
(if (= k 1)
|
||||
(mrg32k3a-random-integer state mrg32k3a-m-max)
|
||||
(+ (* (mrg32k3a-random-power state (- k 1)) mrg32k3a-m-max)
|
||||
(mrg32k3a-random-integer state mrg32k3a-m-max))))
|
||||
|
||||
(define (mrg32k3a-random-large state n) ; n > m-max
|
||||
(do ((k 2 (+ k 1))
|
||||
(mk (* mrg32k3a-m-max mrg32k3a-m-max) (* mk mrg32k3a-m-max)))
|
||||
((>= mk n)
|
||||
(let* ((mk-by-n (quotient mk n))
|
||||
(a (* mk-by-n n)))
|
||||
(do ((x (mrg32k3a-random-power state k)
|
||||
(mrg32k3a-random-power state k)))
|
||||
((< x a) (quotient x mk-by-n)))))))
|
||||
|
||||
|
||||
; Multiple Precision Reals
|
||||
; ========================
|
||||
;
|
||||
; To produce multiple precision reals we produce a large integer value
|
||||
; and convert it into a real value. This value is then normalized.
|
||||
; The precision goal is unit <= 1/(m^k + 1), or 1/unit - 1 <= m^k.
|
||||
; If you know more about the floating point number types of the
|
||||
; Scheme system, this can be improved.
|
||||
|
||||
(define (mrg32k3a-random-real-mp state unit)
|
||||
(do ((k 1 (+ k 1))
|
||||
(u (- (/ 1 unit) 1) (/ u mrg32k3a-m1)))
|
||||
((<= u 1)
|
||||
(/ (exact->inexact (+ (mrg32k3a-random-power state k) 1))
|
||||
(exact->inexact (+ (expt mrg32k3a-m-max k) 1))))))
|
||||
|
||||
|
||||
; Provide the Interface as Specified in the SRFI
|
||||
; ==============================================
|
||||
;
|
||||
; An object of type random-source is a record containing the procedures
|
||||
; as components. The actual state of the generator is stored in the
|
||||
; binding-time environment of make-random-source.
|
||||
|
||||
(define (make-random-source)
|
||||
(let ((state (mrg32k3a-pack-state ; make a new copy
|
||||
(list->vector (vector->list mrg32k3a-initial-state)))))
|
||||
(:random-source-make
|
||||
(lambda ()
|
||||
(mrg32k3a-state-ref state))
|
||||
(lambda (new-state)
|
||||
(set! state (mrg32k3a-state-set new-state)))
|
||||
(lambda ()
|
||||
(set! state (mrg32k3a-randomize-state state)))
|
||||
(lambda (i j)
|
||||
(set! state (mrg32k3a-pseudo-randomize-state i j)))
|
||||
(lambda ()
|
||||
(lambda (n)
|
||||
(cond
|
||||
((not (and (integer? n) (exact? n) (positive? n)))
|
||||
(error "range must be exact positive integer" n))
|
||||
((<= n mrg32k3a-m-max)
|
||||
(mrg32k3a-random-integer state n))
|
||||
(else
|
||||
(mrg32k3a-random-large state n)))))
|
||||
(lambda args
|
||||
(cond
|
||||
((null? args)
|
||||
(lambda ()
|
||||
(mrg32k3a-random-real state)))
|
||||
((null? (cdr args))
|
||||
(let ((unit (car args)))
|
||||
(cond
|
||||
((not (and (real? unit) (< 0 unit 1)))
|
||||
(error "unit must be real in (0,1)" unit))
|
||||
((<= (- (/ 1 unit) 1) mrg32k3a-m1)
|
||||
(lambda ()
|
||||
(mrg32k3a-random-real state)))
|
||||
(else
|
||||
(lambda ()
|
||||
(mrg32k3a-random-real-mp state unit))))))
|
||||
(else
|
||||
(error "illegal arguments" args)))))))
|
||||
|
||||
(define random-source?
|
||||
:random-source?)
|
||||
|
||||
(define (random-source-state-ref s)
|
||||
((:random-source-state-ref s)))
|
||||
|
||||
(define (random-source-state-set! s state)
|
||||
((:random-source-state-set! s) state))
|
||||
|
||||
(define (random-source-randomize! s)
|
||||
((:random-source-randomize! s)))
|
||||
|
||||
(define (random-source-pseudo-randomize! s i j)
|
||||
((:random-source-pseudo-randomize! s) i j))
|
||||
|
||||
; ---
|
||||
|
||||
(define (random-source-make-integers s)
|
||||
((:random-source-make-integers s)))
|
||||
|
||||
(define (random-source-make-reals s . unit)
|
||||
(apply (:random-source-make-reals s) unit))
|
||||
|
||||
; ---
|
||||
|
||||
(define default-random-source
|
||||
(make-random-source))
|
||||
|
||||
(define random-integer
|
||||
(random-source-make-integers default-random-source))
|
||||
|
||||
(define random-real
|
||||
(random-source-make-reals default-random-source))
|
||||
|
|
@ -0,0 +1,33 @@
|
|||
(define format
|
||||
(lambda (format-string . objects)
|
||||
(let ((buffer (open-output-string)))
|
||||
(let loop ((format-list (string->list format-string))
|
||||
(objects objects))
|
||||
(cond ((null? format-list) (get-output-string buffer))
|
||||
((char=? (car format-list) #\~)
|
||||
(if (null? (cdr format-list))
|
||||
(error 'format "Incomplete escape sequence")
|
||||
(case (cadr format-list)
|
||||
((#\a)
|
||||
(if (null? objects)
|
||||
(error 'format "No value for escape sequence")
|
||||
(begin
|
||||
(display (car objects) buffer)
|
||||
(loop (cddr format-list) (cdr objects)))))
|
||||
((#\s)
|
||||
(if (null? objects)
|
||||
(error 'format "No value for escape sequence")
|
||||
(begin
|
||||
(write (car objects) buffer)
|
||||
(loop (cddr format-list) (cdr objects)))))
|
||||
((#\%)
|
||||
(newline buffer)
|
||||
(loop (cddr format-list) objects))
|
||||
((#\~)
|
||||
(write-char #\~ buffer)
|
||||
(loop (cddr format-list) objects))
|
||||
(else
|
||||
(error 'format "Unrecognized escape sequence")))))
|
||||
(else (write-char (car format-list) buffer)
|
||||
(loop (cdr format-list) objects)))))))
|
||||
|
||||
Loading…
Reference in New Issue