From 024d938e39251b7bb5c9259693454760d1c5ad9b Mon Sep 17 00:00:00 2001 From: mainzelm Date: Wed, 11 Sep 2002 14:34:58 +0000 Subject: [PATCH] Added SRFI 25, 26, 27, 28, 30. --- c/srfi/srfi-27.c | 227 ++++++ scheme/more-interfaces.scm | 12 + scheme/more-packages.scm | 46 +- scheme/srfi/srfi-25.scm | 1380 ++++++++++++++++++++++++++++++++++++ scheme/srfi/srfi-26.scm | 96 +++ scheme/srfi/srfi-27.scm | 569 +++++++++++++++ scheme/srfi/srfi-28.scm | 33 + 7 files changed, 2360 insertions(+), 3 deletions(-) create mode 100644 c/srfi/srfi-27.c create mode 100644 scheme/srfi/srfi-25.scm create mode 100644 scheme/srfi/srfi-26.scm create mode 100644 scheme/srfi/srfi-27.scm create mode 100644 scheme/srfi/srfi-28.scm diff --git a/c/srfi/srfi-27.c b/c/srfi/srfi-27.c new file mode 100644 index 0000000..f864672 --- /dev/null +++ b/c/srfi/srfi-27.c @@ -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 + +/* 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); +} + diff --git a/scheme/more-interfaces.scm b/scheme/more-interfaces.scm index 9c69984..b8301db 100644 --- a/scheme/more-interfaces.scm +++ b/scheme/more-interfaces.scm @@ -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)) \ No newline at end of file diff --git a/scheme/more-packages.scm b/scheme/more-packages.scm index e1ca525..bc64276 100644 --- a/scheme/more-packages.scm +++ b/scheme/more-packages.scm @@ -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))) diff --git a/scheme/srfi/srfi-25.scm b/scheme/srfi/srfi-25.scm new file mode 100644 index 0000000..3ece129 --- /dev/null +++ b/scheme/srfi/srfi-25.scm @@ -0,0 +1,1380 @@ +;;; array +;;; 1997 - 2001 Jussi Piitulainen + + +;; This file is the result of +;; cat array.scm as-srfi-9-record.scm ix-ctor.scm op-ctor.scm > srfi-25.scm + + +;;; --- Intro --- + +;;; This interface to arrays is based on Alan Bawden's array.scm of +;;; 1993 (earlier version in the Internet Repository and another +;;; version in SLIB). This is a complete rewrite, to be consistent +;;; with the rest of Scheme and to make arrays independent of lists. + +;;; Some modifications are due to discussion in srfi-25 mailing list. + +;;; (array? obj) +;;; (make-array shape [obj]) changed arguments +;;; (shape bound ...) new +;;; (array shape obj ...) new +;;; (array-rank array) changed name back +;;; (array-start array dimension) new +;;; (array-end array dimension) new +;;; (array-ref array k ...) +;;; (array-ref array index) new variant +;;; (array-set! array k ... obj) changed argument order +;;; (array-set! array index obj) new variant +;;; (share-array array shape proc) changed arguments + +;;; All other variables in this file have names in "array:". + +;;; Should there be a way to make arrays with initial values mapped +;;; from indices? Sure. The current "initial object" is lame. +;;; +;;; Removed (array-shape array) from here. There is a new version +;;; in arlib though. + +;;; --- Representation type dependencies --- + +;;; The mapping from array indices to the index to the underlying vector +;;; is whatever array:optimize returns. The file "opt" provides three +;;; representations: +;;; +;;; mbda) mapping is a procedure that allows an optional argument +;;; tter) mapping is two procedures that takes exactly the indices +;;; ctor) mapping is a vector of a constant term and coefficients +;;; +;;; Choose one in "opt" to make the optimizer. Then choose the matching +;;; implementation of array-ref and array-set!. +;;; +;;; These should be made macros to inline them. Or have a good compiler +;;; and plant the package as a module. + +;;; 1. Pick an optimizer. +;;; 2. Pick matching index representation. +;;; 3. Pick a record implementation; as-procedure is generic; syntax inlines. +;;; 3. This file is otherwise portable. + +;;; --- Portable R5RS (R4RS and multiple values) --- + +;;; (array? obj) +;;; returns #t if `obj' is an array and #t or #f otherwise. + +(define (array? obj) + (array:array? obj)) + +;;; (make-array shape) +;;; (make-array shape obj) +;;; makes array of `shape' with each cell containing `obj' initially. + +(define (make-array shape . rest) + (or (array:good-shape? shape) + (error "make-array: shape is not a shape")) + (apply array:make-array shape rest)) + +(define (array:make-array shape . rest) + (let ((size (array:size shape))) + (array:make + (if (pair? rest) + (apply (lambda (o) (make-vector size o)) rest) + (make-vector size)) + (if (= size 0) + (array:optimize-empty + (vector-ref (array:shape shape) 1)) + (array:optimize + (array:make-index shape) + (vector-ref (array:shape shape) 1))) + (array:shape->vector shape)))) + +;;; (shape bound ...) +;;; makes a shape. Bounds must be an even number of exact, pairwise +;;; non-decreasing integers. Note that any such array can be a shape. + +(define (shape . bounds) + (let ((v (list->vector bounds))) + (or (even? (vector-length v)) + (error (string-append "shape: uneven number of bounds: " + (array:list->string bounds)))) + (let ((shp (array:make + v + (if (pair? bounds) + (array:shape-index) + (array:empty-shape-index)) + (vector 0 (quotient (vector-length v) 2) + 0 2)))) + (or (array:good-shape? shp) + (error (string-append "shape: bounds are not pairwise " + "non-decreasing exact integers: " + (array:list->string bounds)))) + shp))) + +;;; (array shape obj ...) +;;; is analogous to `vector'. + +(define (array shape . elts) + (or (array:good-shape? shape) + (error (string-append "array: shape " (array:thing->string shape) + " is not a shape"))) + (let ((size (array:size shape))) + (let ((vector (list->vector elts))) + (or (= (vector-length vector) size) + (error (string-append "array: an array of shape " + (array:shape-vector->string + (array:vector shape)) + " has " + (number->string size) + " elements but got " + (number->string (vector-length vector)) + " values: " + (array:list->string elts)))) + (array:make + vector + (if (= size 0) + (array:optimize-empty + (vector-ref (array:shape shape) 1)) + (array:optimize + (array:make-index shape) + (vector-ref (array:shape shape) 1))) + (array:shape->vector shape))))) + +;;; (array-rank array) +;;; returns the number of dimensions of `array'. + +(define (array-rank array) + (quotient (vector-length (array:shape array)) 2)) + +;;; (array-start array k) +;;; returns the lower bound index of array along dimension k. This is +;;; the least valid index along that dimension if the dimension is not +;;; empty. + +(define (array-start array d) + (vector-ref (array:shape array) (+ d d))) + +;;; (array-end array k) +;;; returns the upper bound index of array along dimension k. This is +;;; not a valid index. If the dimension is empty, this is the same as +;;; the lower bound along it. + +(define (array-end array d) + (vector-ref (array:shape array) (+ d d 1))) + +;;; (share-array array shape proc) +;;; makes an array that shares elements of `array' at shape `shape'. +;;; The arguments to `proc' are indices of the result. The values of +;;; `proc' are indices of `array'. + +;;; Todo: in the error message, should recognise the mapping and show it. + +(define (share-array array subshape f) + (or (array:good-shape? subshape) + (error (string-append "share-array: shape " + (array:thing->string subshape) + " is not a shape"))) + (let ((subsize (array:size subshape))) + (or (array:good-share? subshape subsize f (array:shape array)) + (error (string-append "share-array: subshape " + (array:shape-vector->string + (array:vector subshape)) + " does not map into supershape " + (array:shape-vector->string + (array:shape array)) + " under mapping " + (array:map->string + f + (vector-ref (array:shape subshape) 1))))) + (let ((g (array:index array))) + (array:make + (array:vector array) + (if (= subsize 0) + (array:optimize-empty + (vector-ref (array:shape subshape) 1)) + (array:optimize + (lambda ks + (call-with-values + (lambda () (apply f ks)) + (lambda ks (array:vector-index g ks)))) + (vector-ref (array:shape subshape) 1))) + (array:shape->vector subshape))))) + +;;; --- Hrmph --- + +;;; (array:share/index! ...) +;;; reuses a user supplied index object when recognising the +;;; mapping. The mind balks at the very nasty side effect that +;;; exposes the implementation. So this is not in the spec. +;;; But letting index objects in at all creates a pressure +;;; to go the whole hog. Arf. + +;;; Use array:optimize-empty for an empty array to get a +;;; clearly invalid vector index. + +;;; Surely it's perverse to use an actor for index here? But +;;; the possibility is provided for completeness. + +(define (array:share/index! array subshape proc index) + (array:make + (array:vector array) + (if (= (array:size subshape) 0) + (array:optimize-empty + (quotient (vector-length (array:shape array)) 2)) + ((if (vector? index) + array:optimize/vector + array:optimize/actor) + (lambda (subindex) + (let ((superindex (proc subindex))) + (if (vector? superindex) + (array:index/vector + (quotient (vector-length (array:shape array)) 2) + (array:index array) + superindex) + (array:index/array + (quotient (vector-length (array:shape array)) 2) + (array:index array) + (array:vector superindex) + (array:index superindex))))) + index)) + (array:shape->vector subshape))) + +(define (array:optimize/vector f v) + (let ((r (vector-length v))) + (do ((k 0 (+ k 1))) + ((= k r)) + (vector-set! v k 0)) + (let ((n0 (f v)) + (cs (make-vector (+ r 1))) + (apply (array:applier-to-vector (+ r 1)))) + (vector-set! cs 0 n0) + (let wok ((k 0)) + (if (< k r) + (let ((k1 (+ k 1))) + (vector-set! v k 1) + (let ((nk (- (f v) n0))) + (vector-set! v k 0) + (vector-set! cs k1 nk) + (wok k1))))) + (apply (array:maker r) cs)))) + +(define (array:optimize/actor f a) + (let ((r (array-end a 0)) + (v (array:vector a)) + (i (array:index a))) + (do ((k 0 (+ k 1))) + ((= k r)) + (vector-set! v (array:actor-index i k) 0)) + (let ((n0 (f a)) + (cs (make-vector (+ r 1))) + (apply (array:applier-to-vector (+ r 1)))) + (vector-set! cs 0 n0) + (let wok ((k 0)) + (if (< k r) + (let ((k1 (+ k 1)) + (t (array:actor-index i k))) + (vector-set! v t 1) + (let ((nk (- (f a) n0))) + (vector-set! v t 0) + (vector-set! cs k1 nk) + (wok k1))))) + (apply (array:maker r) cs)))) + +;;; --- Internals --- + +(define (array:shape->vector shape) + (let ((idx (array:index shape)) + (shv (array:vector shape)) + (rnk (vector-ref (array:shape shape) 1))) + (let ((vec (make-vector (* rnk 2)))) + (do ((k 0 (+ k 1))) + ((= k rnk) + vec) + (vector-set! vec (+ k k) + (vector-ref shv (array:shape-vector-index idx k 0))) + (vector-set! vec (+ k k 1) + (vector-ref shv (array:shape-vector-index idx k 1))))))) + +;;; (array:size shape) +;;; returns the number of elements in arrays of shape `shape'. + +(define (array:size shape) + (let ((idx (array:index shape)) + (shv (array:vector shape)) + (rnk (vector-ref (array:shape shape) 1))) + (do ((k 0 (+ k 1)) + (s 1 (* s + (- (vector-ref shv (array:shape-vector-index idx k 1)) + (vector-ref shv (array:shape-vector-index idx k 0)))))) + ((= k rnk) s)))) + +;;; (array:make-index shape) +;;; returns an index function for arrays of shape `shape'. This is a +;;; runtime composition of several variable arity procedures, to be +;;; passed to array:optimize for recognition as an affine function of +;;; as many variables as there are dimensions in arrays of this shape. + +(define (array:make-index shape) + (let ((idx (array:index shape)) + (shv (array:vector shape)) + (rnk (vector-ref (array:shape shape) 1))) + (do ((f (lambda () 0) + (lambda (k . ks) + (+ (* s (- k (vector-ref + shv + (array:shape-vector-index idx (- j 1) 0)))) + (apply f ks)))) + (s 1 (* s (- (vector-ref + shv + (array:shape-vector-index idx (- j 1) 1)) + (vector-ref + shv + (array:shape-vector-index idx (- j 1) 0))))) + (j rnk (- j 1))) + ((= j 0) + f)))) + + +;;; --- Error checking --- + +;;; (array:good-shape? shape) +;;; returns true if `shape' is an array of the right shape and its +;;; elements are exact integers that pairwise bound intervals `[lo..hi)´. + +(define (array:good-shape? shape) + (and (array:array? shape) + (let ((u (array:shape shape)) + (v (array:vector shape)) + (x (array:index shape))) + (and (= (vector-length u) 4) + (= (vector-ref u 0) 0) + (= (vector-ref u 2) 0) + (= (vector-ref u 3) 2)) + (let ((p (vector-ref u 1))) + (do ((k 0 (+ k 1)) + (true #t (let ((lo (vector-ref + v + (array:shape-vector-index x k 0))) + (hi (vector-ref + v + (array:shape-vector-index x k 1)))) + (and true + (integer? lo) + (exact? lo) + (integer? hi) + (exact? hi) + (<= lo hi))))) + ((= k p) true)))))) + +;;; (array:good-share? subv subsize mapping superv) +;;; returns true if the extreme indices in the subshape vector map +;;; into the bounds in the supershape vector. + +;;; If some interval in `subv' is empty, then `subv' is empty and its +;;; image under `f' is empty and it is trivially alright. One must +;;; not call `f', though. + +(define (array:good-share? subshape subsize f super) + (or (zero? subsize) + (letrec + ((sub (array:vector subshape)) + (dex (array:index subshape)) + (ck (lambda (k ks) + (if (zero? k) + (call-with-values + (lambda () (apply f ks)) + (lambda qs (array:good-indices? qs super))) + (and (ck (- k 1) + (cons (vector-ref + sub + (array:shape-vector-index + dex + (- k 1) + 0)) + ks)) + (ck (- k 1) + (cons (- (vector-ref + sub + (array:shape-vector-index + dex + (- k 1) + 1)) + 1) + ks))))))) + (let ((rnk (vector-ref (array:shape subshape) 1))) + (or (array:unchecked-share-depth? rnk) + (ck rnk '())))))) + +;;; Check good-share on 10 dimensions at most. The trouble is, +;;; the cost of this check is exponential in the number of dimensions. + +(define (array:unchecked-share-depth? rank) + (if (> rank 10) + (begin + (display `(warning: unchecked depth in share: + ,rank subdimensions)) + (newline) + #t) + #f)) + +;;; (array:check-indices caller indices shape-vector) +;;; (array:check-indices.o caller indices shape-vector) +;;; (array:check-index-vector caller index-vector shape-vector) +;;; return if the index is in bounds, else signal error. +;;; +;;; Shape-vector is the internal representation, with +;;; b and e for dimension k at 2k and 2k + 1. + +(define (array:check-indices who ks shv) + (or (array:good-indices? ks shv) + (error (array:not-in who ks shv)))) + +(define (array:check-indices.o who ks shv) + (or (array:good-indices.o? ks shv) + (error (array:not-in who (reverse (cdr (reverse ks))) shv)))) + +(define (array:check-index-vector who ks shv) + (or (array:good-index-vector? ks shv) + (error (array:not-in who (vector->list ks) shv)))) + +(define (array:check-index-actor who ks shv) + (let ((shape (array:shape ks))) + (or (and (= (vector-length shape) 2) + (= (vector-ref shape 0) 0)) + (error "not an actor")) + (or (array:good-index-actor? + (vector-ref shape 1) + (array:vector ks) + (array:index ks) + shv) + (array:not-in who (do ((k (vector-ref shape 1) (- k 1)) + (m '() (cons (vector-ref + (array:vector ks) + (array:actor-index + (array:index ks) + (- k 1))) + m))) + ((= k 0) m)) + shv)))) + +(define (array:good-indices? ks shv) + (let ((d2 (vector-length shv))) + (do ((kp ks (if (pair? kp) + (cdr kp))) + (k 0 (+ k 2)) + (true #t (and true (pair? kp) + (array:good-index? (car kp) shv k)))) + ((= k d2) + (and true (null? kp)))))) + +(define (array:good-indices.o? ks.o shv) + (let ((d2 (vector-length shv))) + (do ((kp ks.o (if (pair? kp) + (cdr kp))) + (k 0 (+ k 2)) + (true #t (and true (pair? kp) + (array:good-index? (car kp) shv k)))) + ((= k d2) + (and true (pair? kp) (null? (cdr kp))))))) + +(define (array:good-index-vector? ks shv) + (let ((r2 (vector-length shv))) + (and (= (* 2 (vector-length ks)) r2) + (do ((j 0 (+ j 1)) + (k 0 (+ k 2)) + (true #t (and true + (array:good-index? (vector-ref ks j) shv k)))) + ((= k r2) true))))) + +(define (array:good-index-actor? r v i shv) + (and (= (* 2 r) (vector-length shv)) + (do ((j 0 (+ j 1)) + (k 0 (+ k 2)) + (true #t (and true + (array:good-index? (vector-ref + v + (array:actor-index i j)) + shv + k)))) + ((= j r) true)))) + +;;; (array:good-index? index shape-vector 2d) +;;; returns true if index is within bounds for dimension 2d/2. + +(define (array:good-index? w shv k) + (and (integer? w) + (exact? w) + (<= (vector-ref shv k) w) + (< w (vector-ref shv (+ k 1))))) + +(define (array:not-in who ks shv) + (let ((index (array:list->string ks)) + (bounds (array:shape-vector->string shv))) + (error (string-append who + ": index " index + " not in bounds " bounds)))) + +(define (array:list->string ks) + (do ((index "" (string-append index (array:thing->string (car ks)) " ")) + (ks ks (cdr ks))) + ((null? ks) index))) + +(define (array:shape-vector->string shv) + (do ((bounds "" (string-append bounds + "[" + (number->string (vector-ref shv t)) + ".." + (number->string (vector-ref shv (+ t 1))) + ")" + " ")) + (t 0 (+ t 2))) + ((= t (vector-length shv)) bounds))) + +(define (array:thing->string thing) + (cond + ((number? thing) (number->string thing)) + ((symbol? thing) (string-append "#" (symbol->string thing))) + ((char? thing) "#") + ((string? thing) "#") + ((list? thing) (string-append "#" (number->string (length thing)) + "")) + + ((pair? thing) "#") + ((array? thing) "#") + ((vector? thing) (string-append "#" (number->string + (vector-length thing)) + "")) + ((procedure? thing) "#") + (else + (case thing + ((()) "()") + ((#t) "#t") + ((#f) "#f") + (else + "#"))))) + +;;; And to grok an affine map, vector->vector type. Column k of arr +;;; will contain coefficients n0 ... nm of 1 k1 ... km for kth value. +;;; +;;; These are for the error message when share fails. + +(define (array:index-ref ind k) + (if (vector? ind) + (vector-ref ind k) + (vector-ref + (array:vector ind) + (array:actor-index (array:index ind) k)))) + +(define (array:index-set! ind k o) + (if (vector? ind) + (vector-set! ind k o) + (vector-set! + (array:vector ind) + (array:actor-index (array:index ind) k) + o))) + +(define (array:index-length ind) + (if (vector? ind) + (vector-length ind) + (vector-ref (array:shape ind) 1))) + +(define (array:map->string proc r) + (let* ((m (array:grok/arguments proc r)) + (s (vector-ref (array:shape m) 3))) + (do ((i "" (string-append i c "k" (number->string k))) + (c "" ", ") + (k 1 (+ k 1))) + ((< r k) + (do ((o "" (string-append o c (array:map-column->string m r k))) + (c "" ", ") + (k 0 (+ k 1))) + ((= k s) + (string-append i " => " o))))))) + +(define (array:map-column->string m r k) + (let ((v (array:vector m)) + (i (array:index m))) + (let ((n0 (vector-ref v (array:vector-index i (list 0 k))))) + (let wok ((j 1) + (e (if (= n0 0) "" (number->string n0)))) + (if (<= j r) + (let ((nj (vector-ref v (array:vector-index i (list j k))))) + (if (= nj 0) + (wok (+ j 1) e) + (let* ((nj (if (= nj 1) "" + (if (= nj -1) "-" + (string-append (number->string nj) + " ")))) + (njkj (string-append nj "k" (number->string j)))) + (if (string=? e "") + (wok (+ j 1) njkj) + (wok (+ j 1) (string-append e " + " njkj)))))) + (if (string=? e "") "0" e)))))) + +(define (array:grok/arguments proc r) + (array:grok/index! + (lambda (vec) + (call-with-values + (lambda () + (array:apply-to-vector r proc vec)) + vector)) + (make-vector r))) + +(define (array:grok/index! proc in) + (let ((m (array:index-length in))) + (do ((k 0 (+ k 1))) + ((= k m)) + (array:index-set! in k 0)) + (let* ((n0 (proc in)) + (n (array:index-length n0))) + (let ((arr (make-array (shape 0 (+ m 1) 0 n)))) ; (*) + (do ((k 0 (+ k 1))) + ((= k n)) + (array-set! arr 0 k (array:index-ref n0 k))) ; (**) + (do ((j 0 (+ j 1))) + ((= j m)) + (array:index-set! in j 1) + (let ((nj (proc in))) + (array:index-set! in j 0) + (do ((k 0 (+ k 1))) + ((= k n)) + (array-set! arr (+ j 1) k (- (array:index-ref nj k) ; (**) + (array:index-ref n0 k)))))) + arr)))) +;; (*) Should not use `make-array' and `shape' here +;; (**) Should not use `array-set!' here +;; Should use something internal to the library instead: either lower +;; level code (preferable but complex) or alternative names to these same. +;;; array as-srfi-9-record +;;; 2001 Jussi Piitulainen + +;;; Untested. + +(define-record-type + array:srfi-9-record-type-descriptor + (array:make vec ind shp) + array:array? + (vec array:vector) + (ind array:index) + (shp array:shape)) +(define (array-ref a . xs) + (or (array:array? a) + (error "not an array")) + (let ((shape (array:shape a))) + (if (null? xs) + (array:check-indices "array-ref" xs shape) + (let ((x (car xs))) + (if (vector? x) + (array:check-index-vector "array-ref" x shape) + (if (integer? x) + (array:check-indices "array-ref" xs shape) + (if (array:array? x) + (array:check-index-actor "array-ref" x shape) + (error "not an index object")))))) + (vector-ref + (array:vector a) + (if (null? xs) + (vector-ref (array:index a) 0) + (let ((x (car xs))) + (if (vector? x) + (array:index/vector + (quotient (vector-length shape) 2) + (array:index a) + x) + (if (integer? x) + (array:vector-index (array:index a) xs) + (if (array:array? x) + (array:index/array + (quotient (vector-length shape) 2) + (array:index a) + (array:vector x) + (array:index x)) + (error "array-ref: bad index object"))))))))) + +(define (array-set! a x . xs) + (or (array:array? a) + (error "array-set!: not an array")) + (let ((shape (array:shape a))) + (if (null? xs) + (array:check-indices "array-set!" '() shape) + (if (vector? x) + (array:check-index-vector "array-set!" x shape) + (if (integer? x) + (array:check-indices.o "array-set!" (cons x xs) shape) + (if (array:array? x) + (array:check-index-actor "array-set!" x shape) + (error "not an index object"))))) + (if (null? xs) + (vector-set! (array:vector a) (vector-ref (array:index a) 0) x) + (if (vector? x) + (vector-set! (array:vector a) + (array:index/vector + (quotient (vector-length shape) 2) + (array:index a) + x) + (car xs)) + (if (integer? x) + (let ((v (array:vector a)) + (i (array:index a)) + (r (quotient (vector-length shape) 2))) + (do ((sum (* (vector-ref i 0) x) + (+ sum (* (vector-ref i k) (car ks)))) + (ks xs (cdr ks)) + (k 1 (+ k 1))) + ((= k r) + (vector-set! v (+ sum (vector-ref i k)) (car ks))))) + (if (array:array? x) + (vector-set! (array:vector a) + (array:index/array + (quotient (vector-length shape) 2) + (array:index a) + (array:vector x) + (array:index x)) + (car xs)) + (error (string-append + "array-set!: bad index object: " + (array:thing->string x))))))))) +(begin + (define array:opt-args '(ctor (4))) + (define (array:optimize f r) + (case r + ((0) (let ((n0 (f))) (array:0 n0))) + ((1) (let ((n0 (f 0))) (array:1 n0 (- (f 1) n0)))) + ((2) + (let ((n0 (f 0 0))) + (array:2 n0 (- (f 1 0) n0) (- (f 0 1) n0)))) + ((3) + (let ((n0 (f 0 0 0))) + (array:3 + n0 + (- (f 1 0 0) n0) + (- (f 0 1 0) n0) + (- (f 0 0 1) n0)))) + (else + (let ((v + (do ((k 0 (+ k 1)) (v '() (cons 0 v))) + ((= k r) v)))) + (let ((n0 (apply f v))) + (apply + array:n + n0 + (array:coefficients f n0 v v))))))) + (define (array:optimize-empty r) + (let ((x (make-vector (+ r 1) 0))) + (vector-set! x r -1) + x)) + (define (array:coefficients f n0 vs vp) + (case vp + ((()) '()) + (else + (set-car! vp 1) + (let ((n (- (apply f vs) n0))) + (set-car! vp 0) + (cons n (array:coefficients f n0 vs (cdr vp))))))) + (define (array:vector-index x ks) + (do ((sum 0 (+ sum (* (vector-ref x k) (car ks)))) + (ks ks (cdr ks)) + (k 0 (+ k 1))) + ((null? ks) (+ sum (vector-ref x k))))) + (define (array:shape-index) '#(2 1 0)) + (define (array:empty-shape-index) '#(0 0 -1)) + (define (array:shape-vector-index x r k) + (+ + (* (vector-ref x 0) r) + (* (vector-ref x 1) k) + (vector-ref x 2))) + (define (array:actor-index x k) + (+ (* (vector-ref x 0) k) (vector-ref x 1))) + (define (array:0 n0) (vector n0)) + (define (array:1 n0 n1) (vector n1 n0)) + (define (array:2 n0 n1 n2) (vector n1 n2 n0)) + (define (array:3 n0 n1 n2 n3) (vector n1 n2 n3 n0)) + (define (array:n n0 n1 n2 n3 n4 . ns) + (apply vector n1 n2 n3 n4 (append ns (list n0)))) + (define (array:maker r) + (case r + ((0) array:0) + ((1) array:1) + ((2) array:2) + ((3) array:3) + (else array:n))) + (define array:indexer/vector + (let ((em + (vector + (lambda (x i) (+ (vector-ref x 0))) + (lambda (x i) + (+ + (* (vector-ref x 0) (vector-ref i 0)) + (vector-ref x 1))) + (lambda (x i) + (+ + (* (vector-ref x 0) (vector-ref i 0)) + (* (vector-ref x 1) (vector-ref i 1)) + (vector-ref x 2))) + (lambda (x i) + (+ + (* (vector-ref x 0) (vector-ref i 0)) + (* (vector-ref x 1) (vector-ref i 1)) + (* (vector-ref x 2) (vector-ref i 2)) + (vector-ref x 3))) + (lambda (x i) + (+ + (* (vector-ref x 0) (vector-ref i 0)) + (* (vector-ref x 1) (vector-ref i 1)) + (* (vector-ref x 2) (vector-ref i 2)) + (* (vector-ref x 3) (vector-ref i 3)) + (vector-ref x 4))) + (lambda (x i) + (+ + (* (vector-ref x 0) (vector-ref i 0)) + (* (vector-ref x 1) (vector-ref i 1)) + (* (vector-ref x 2) (vector-ref i 2)) + (* (vector-ref x 3) (vector-ref i 3)) + (* (vector-ref x 4) (vector-ref i 4)) + (vector-ref x 5))) + (lambda (x i) + (+ + (* (vector-ref x 0) (vector-ref i 0)) + (* (vector-ref x 1) (vector-ref i 1)) + (* (vector-ref x 2) (vector-ref i 2)) + (* (vector-ref x 3) (vector-ref i 3)) + (* (vector-ref x 4) (vector-ref i 4)) + (* (vector-ref x 5) (vector-ref i 5)) + (vector-ref x 6))) + (lambda (x i) + (+ + (* (vector-ref x 0) (vector-ref i 0)) + (* (vector-ref x 1) (vector-ref i 1)) + (* (vector-ref x 2) (vector-ref i 2)) + (* (vector-ref x 3) (vector-ref i 3)) + (* (vector-ref x 4) (vector-ref i 4)) + (* (vector-ref x 5) (vector-ref i 5)) + (* (vector-ref x 6) (vector-ref i 6)) + (vector-ref x 7))) + (lambda (x i) + (+ + (* (vector-ref x 0) (vector-ref i 0)) + (* (vector-ref x 1) (vector-ref i 1)) + (* (vector-ref x 2) (vector-ref i 2)) + (* (vector-ref x 3) (vector-ref i 3)) + (* (vector-ref x 4) (vector-ref i 4)) + (* (vector-ref x 5) (vector-ref i 5)) + (* (vector-ref x 6) (vector-ref i 6)) + (* (vector-ref x 7) (vector-ref i 7)) + (vector-ref x 8))) + (lambda (x i) + (+ + (* (vector-ref x 0) (vector-ref i 0)) + (* (vector-ref x 1) (vector-ref i 1)) + (* (vector-ref x 2) (vector-ref i 2)) + (* (vector-ref x 3) (vector-ref i 3)) + (* (vector-ref x 4) (vector-ref i 4)) + (* (vector-ref x 5) (vector-ref i 5)) + (* (vector-ref x 6) (vector-ref i 6)) + (* (vector-ref x 7) (vector-ref i 7)) + (* (vector-ref x 8) (vector-ref i 8)) + (vector-ref x 9))))) + (it + (lambda (w) + (lambda (x i) + (+ + (* (vector-ref x 0) (vector-ref i 0)) + (* (vector-ref x 1) (vector-ref i 1)) + (* (vector-ref x 2) (vector-ref i 2)) + (* (vector-ref x 3) (vector-ref i 3)) + (* (vector-ref x 4) (vector-ref i 4)) + (* (vector-ref x 5) (vector-ref i 5)) + (* (vector-ref x 6) (vector-ref i 6)) + (* (vector-ref x 7) (vector-ref i 7)) + (* (vector-ref x 8) (vector-ref i 8)) + (* (vector-ref x 9) (vector-ref i 9)) + (do ((xi + 0 + (+ + (* (vector-ref x u) (vector-ref i u)) + xi)) + (u (- w 1) (- u 1))) + ((< u 10) xi)) + (vector-ref x w)))))) + (lambda (r) (if (< r 10) (vector-ref em r) (it r))))) + (define array:indexer/array + (let ((em + (vector + (lambda (x v i) (+ (vector-ref x 0))) + (lambda (x v i) + (+ + (* + (vector-ref x 0) + (vector-ref v (array:actor-index i 0))) + (vector-ref x 1))) + (lambda (x v i) + (+ + (* + (vector-ref x 0) + (vector-ref v (array:actor-index i 0))) + (* + (vector-ref x 1) + (vector-ref v (array:actor-index i 1))) + (vector-ref x 2))) + (lambda (x v i) + (+ + (* + (vector-ref x 0) + (vector-ref v (array:actor-index i 0))) + (* + (vector-ref x 1) + (vector-ref v (array:actor-index i 1))) + (* + (vector-ref x 2) + (vector-ref v (array:actor-index i 2))) + (vector-ref x 3))) + (lambda (x v i) + (+ + (* + (vector-ref x 0) + (vector-ref v (array:actor-index i 0))) + (* + (vector-ref x 1) + (vector-ref v (array:actor-index i 1))) + (* + (vector-ref x 2) + (vector-ref v (array:actor-index i 2))) + (* + (vector-ref x 3) + (vector-ref v (array:actor-index i 3))) + (vector-ref x 4))) + (lambda (x v i) + (+ + (* + (vector-ref x 0) + (vector-ref v (array:actor-index i 0))) + (* + (vector-ref x 1) + (vector-ref v (array:actor-index i 1))) + (* + (vector-ref x 2) + (vector-ref v (array:actor-index i 2))) + (* + (vector-ref x 3) + (vector-ref v (array:actor-index i 3))) + (* + (vector-ref x 4) + (vector-ref v (array:actor-index i 4))) + (vector-ref x 5))) + (lambda (x v i) + (+ + (* + (vector-ref x 0) + (vector-ref v (array:actor-index i 0))) + (* + (vector-ref x 1) + (vector-ref v (array:actor-index i 1))) + (* + (vector-ref x 2) + (vector-ref v (array:actor-index i 2))) + (* + (vector-ref x 3) + (vector-ref v (array:actor-index i 3))) + (* + (vector-ref x 4) + (vector-ref v (array:actor-index i 4))) + (* + (vector-ref x 5) + (vector-ref v (array:actor-index i 5))) + (vector-ref x 6))) + (lambda (x v i) + (+ + (* + (vector-ref x 0) + (vector-ref v (array:actor-index i 0))) + (* + (vector-ref x 1) + (vector-ref v (array:actor-index i 1))) + (* + (vector-ref x 2) + (vector-ref v (array:actor-index i 2))) + (* + (vector-ref x 3) + (vector-ref v (array:actor-index i 3))) + (* + (vector-ref x 4) + (vector-ref v (array:actor-index i 4))) + (* + (vector-ref x 5) + (vector-ref v (array:actor-index i 5))) + (* + (vector-ref x 6) + (vector-ref v (array:actor-index i 6))) + (vector-ref x 7))) + (lambda (x v i) + (+ + (* + (vector-ref x 0) + (vector-ref v (array:actor-index i 0))) + (* + (vector-ref x 1) + (vector-ref v (array:actor-index i 1))) + (* + (vector-ref x 2) + (vector-ref v (array:actor-index i 2))) + (* + (vector-ref x 3) + (vector-ref v (array:actor-index i 3))) + (* + (vector-ref x 4) + (vector-ref v (array:actor-index i 4))) + (* + (vector-ref x 5) + (vector-ref v (array:actor-index i 5))) + (* + (vector-ref x 6) + (vector-ref v (array:actor-index i 6))) + (* + (vector-ref x 7) + (vector-ref v (array:actor-index i 7))) + (vector-ref x 8))) + (lambda (x v i) + (+ + (* + (vector-ref x 0) + (vector-ref v (array:actor-index i 0))) + (* + (vector-ref x 1) + (vector-ref v (array:actor-index i 1))) + (* + (vector-ref x 2) + (vector-ref v (array:actor-index i 2))) + (* + (vector-ref x 3) + (vector-ref v (array:actor-index i 3))) + (* + (vector-ref x 4) + (vector-ref v (array:actor-index i 4))) + (* + (vector-ref x 5) + (vector-ref v (array:actor-index i 5))) + (* + (vector-ref x 6) + (vector-ref v (array:actor-index i 6))) + (* + (vector-ref x 7) + (vector-ref v (array:actor-index i 7))) + (* + (vector-ref x 8) + (vector-ref v (array:actor-index i 8))) + (vector-ref x 9))))) + (it + (lambda (w) + (lambda (x v i) + (+ + (* + (vector-ref x 0) + (vector-ref v (array:actor-index i 0))) + (* + (vector-ref x 1) + (vector-ref v (array:actor-index i 1))) + (* + (vector-ref x 2) + (vector-ref v (array:actor-index i 2))) + (* + (vector-ref x 3) + (vector-ref v (array:actor-index i 3))) + (* + (vector-ref x 4) + (vector-ref v (array:actor-index i 4))) + (* + (vector-ref x 5) + (vector-ref v (array:actor-index i 5))) + (* + (vector-ref x 6) + (vector-ref v (array:actor-index i 6))) + (* + (vector-ref x 7) + (vector-ref v (array:actor-index i 7))) + (* + (vector-ref x 8) + (vector-ref v (array:actor-index i 8))) + (* + (vector-ref x 9) + (vector-ref v (array:actor-index i 9))) + (do ((xi + 0 + (+ + (* + (vector-ref x u) + (vector-ref + v + (array:actor-index i u))) + xi)) + (u (- w 1) (- u 1))) + ((< u 10) xi)) + (vector-ref x w)))))) + (lambda (r) (if (< r 10) (vector-ref em r) (it r))))) + (define array:applier-to-vector + (let ((em + (vector + (lambda (p v) (p)) + (lambda (p v) (p (vector-ref v 0))) + (lambda (p v) + (p (vector-ref v 0) (vector-ref v 1))) + (lambda (p v) + (p + (vector-ref v 0) + (vector-ref v 1) + (vector-ref v 2))) + (lambda (p v) + (p + (vector-ref v 0) + (vector-ref v 1) + (vector-ref v 2) + (vector-ref v 3))) + (lambda (p v) + (p + (vector-ref v 0) + (vector-ref v 1) + (vector-ref v 2) + (vector-ref v 3) + (vector-ref v 4))) + (lambda (p v) + (p + (vector-ref v 0) + (vector-ref v 1) + (vector-ref v 2) + (vector-ref v 3) + (vector-ref v 4) + (vector-ref v 5))) + (lambda (p v) + (p + (vector-ref v 0) + (vector-ref v 1) + (vector-ref v 2) + (vector-ref v 3) + (vector-ref v 4) + (vector-ref v 5) + (vector-ref v 6))) + (lambda (p v) + (p + (vector-ref v 0) + (vector-ref v 1) + (vector-ref v 2) + (vector-ref v 3) + (vector-ref v 4) + (vector-ref v 5) + (vector-ref v 6) + (vector-ref v 7))) + (lambda (p v) + (p + (vector-ref v 0) + (vector-ref v 1) + (vector-ref v 2) + (vector-ref v 3) + (vector-ref v 4) + (vector-ref v 5) + (vector-ref v 6) + (vector-ref v 7) + (vector-ref v 8))))) + (it + (lambda (r) + (lambda (p v) + (apply + p + (vector-ref v 0) + (vector-ref v 1) + (vector-ref v 2) + (vector-ref v 3) + (vector-ref v 4) + (vector-ref v 5) + (vector-ref v 6) + (vector-ref v 7) + (vector-ref v 8) + (vector-ref v 9) + (do ((k r (- k 1)) + (r + '() + (cons (vector-ref v (- k 1)) r))) + ((= k 10) r))))))) + (lambda (r) (if (< r 10) (vector-ref em r) (it r))))) + (define array:applier-to-actor + (let ((em + (vector + (lambda (p a) (p)) + (lambda (p a) (p (array-ref a 0))) + (lambda (p a) + (p (array-ref a 0) (array-ref a 1))) + (lambda (p a) + (p + (array-ref a 0) + (array-ref a 1) + (array-ref a 2))) + (lambda (p a) + (p + (array-ref a 0) + (array-ref a 1) + (array-ref a 2) + (array-ref a 3))) + (lambda (p a) + (p + (array-ref a 0) + (array-ref a 1) + (array-ref a 2) + (array-ref a 3) + (array-ref a 4))) + (lambda (p a) + (p + (array-ref a 0) + (array-ref a 1) + (array-ref a 2) + (array-ref a 3) + (array-ref a 4) + (array-ref a 5))) + (lambda (p a) + (p + (array-ref a 0) + (array-ref a 1) + (array-ref a 2) + (array-ref a 3) + (array-ref a 4) + (array-ref a 5) + (array-ref a 6))) + (lambda (p a) + (p + (array-ref a 0) + (array-ref a 1) + (array-ref a 2) + (array-ref a 3) + (array-ref a 4) + (array-ref a 5) + (array-ref a 6) + (array-ref a 7))) + (lambda (p a) + (p + (array-ref a 0) + (array-ref a 1) + (array-ref a 2) + (array-ref a 3) + (array-ref a 4) + (array-ref a 5) + (array-ref a 6) + (array-ref a 7) + (array-ref a 8))))) + (it + (lambda (r) + (lambda (p a) + (apply + a + (array-ref a 0) + (array-ref a 1) + (array-ref a 2) + (array-ref a 3) + (array-ref a 4) + (array-ref a 5) + (array-ref a 6) + (array-ref a 7) + (array-ref a 8) + (array-ref a 9) + (do ((k r (- k 1)) + (r '() (cons (array-ref a (- k 1)) r))) + ((= k 10) r))))))) + (lambda (r) + "These are high level, hiding implementation at call site." + (if (< r 10) (vector-ref em r) (it r))))) + (define array:applier-to-backing-vector + (let ((em + (vector + (lambda (p ai av) (p)) + (lambda (p ai av) + (p (vector-ref av (array:actor-index ai 0)))) + (lambda (p ai av) + (p + (vector-ref av (array:actor-index ai 0)) + (vector-ref av (array:actor-index ai 1)))) + (lambda (p ai av) + (p + (vector-ref av (array:actor-index ai 0)) + (vector-ref av (array:actor-index ai 1)) + (vector-ref av (array:actor-index ai 2)))) + (lambda (p ai av) + (p + (vector-ref av (array:actor-index ai 0)) + (vector-ref av (array:actor-index ai 1)) + (vector-ref av (array:actor-index ai 2)) + (vector-ref av (array:actor-index ai 3)))) + (lambda (p ai av) + (p + (vector-ref av (array:actor-index ai 0)) + (vector-ref av (array:actor-index ai 1)) + (vector-ref av (array:actor-index ai 2)) + (vector-ref av (array:actor-index ai 3)) + (vector-ref av (array:actor-index ai 4)))) + (lambda (p ai av) + (p + (vector-ref av (array:actor-index ai 0)) + (vector-ref av (array:actor-index ai 1)) + (vector-ref av (array:actor-index ai 2)) + (vector-ref av (array:actor-index ai 3)) + (vector-ref av (array:actor-index ai 4)) + (vector-ref av (array:actor-index ai 5)))) + (lambda (p ai av) + (p + (vector-ref av (array:actor-index ai 0)) + (vector-ref av (array:actor-index ai 1)) + (vector-ref av (array:actor-index ai 2)) + (vector-ref av (array:actor-index ai 3)) + (vector-ref av (array:actor-index ai 4)) + (vector-ref av (array:actor-index ai 5)) + (vector-ref av (array:actor-index ai 6)))) + (lambda (p ai av) + (p + (vector-ref av (array:actor-index ai 0)) + (vector-ref av (array:actor-index ai 1)) + (vector-ref av (array:actor-index ai 2)) + (vector-ref av (array:actor-index ai 3)) + (vector-ref av (array:actor-index ai 4)) + (vector-ref av (array:actor-index ai 5)) + (vector-ref av (array:actor-index ai 6)) + (vector-ref av (array:actor-index ai 7)))) + (lambda (p ai av) + (p + (vector-ref av (array:actor-index ai 0)) + (vector-ref av (array:actor-index ai 1)) + (vector-ref av (array:actor-index ai 2)) + (vector-ref av (array:actor-index ai 3)) + (vector-ref av (array:actor-index ai 4)) + (vector-ref av (array:actor-index ai 5)) + (vector-ref av (array:actor-index ai 6)) + (vector-ref av (array:actor-index ai 7)) + (vector-ref av (array:actor-index ai 8)))))) + (it + (lambda (r) + (lambda (p ai av) + (apply + p + (vector-ref av (array:actor-index ai 0)) + (vector-ref av (array:actor-index ai 1)) + (vector-ref av (array:actor-index ai 2)) + (vector-ref av (array:actor-index ai 3)) + (vector-ref av (array:actor-index ai 4)) + (vector-ref av (array:actor-index ai 5)) + (vector-ref av (array:actor-index ai 6)) + (vector-ref av (array:actor-index ai 7)) + (vector-ref av (array:actor-index ai 8)) + (vector-ref av (array:actor-index ai 9)) + (do ((k r (- k 1)) + (r + '() + (cons + (vector-ref + av + (array:actor-index ai (- k 1))) + r))) + ((= k 10) r))))))) + (lambda (r) + "These are low level, exposing implementation at call site." + (if (< r 10) (vector-ref em r) (it r))))) + (define (array:index/vector r x v) + ((array:indexer/vector r) x v)) + (define (array:index/array r x av ai) + ((array:indexer/array r) x av ai)) + (define (array:apply-to-vector r p v) + ((array:applier-to-vector r) p v)) + (define (array:apply-to-actor r p a) + ((array:applier-to-actor r) p a))) diff --git a/scheme/srfi/srfi-26.scm b/scheme/srfi/srfi-26.scm new file mode 100644 index 0000000..bca3e9d --- /dev/null +++ b/scheme/srfi/srfi-26.scm @@ -0,0 +1,96 @@ +; REFERENCE IMPLEMENTATION FOR SRFI-26 "CUT" +; ========================================== +; +; Sebastian.Egner@philips.com, 5-Jun-2002. +; adapted from the posting by Al Petrofsky +; 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 + ((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)))) diff --git a/scheme/srfi/srfi-27.scm b/scheme/srfi/srfi-27.scm new file mode 100644 index 0000000..93c792d --- /dev/null +++ b/scheme/srfi/srfi-27.scm @@ -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)) diff --git a/scheme/srfi/srfi-28.scm b/scheme/srfi/srfi-28.scm new file mode 100644 index 0000000..3a2d765 --- /dev/null +++ b/scheme/srfi/srfi-28.scm @@ -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))))))) +