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 <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);
+}
+
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>" (symbol->string thing)))
+    ((char? thing) "#<char>")
+    ((string? thing) "#<string>")
+    ((list? thing) (string-append "#" (number->string (length thing))
+                                  "<list>"))
+                                  
+    ((pair? thing) "#<pair>")
+    ((array? thing) "#<array>")
+    ((vector? thing) (string-append "#" (number->string
+                                         (vector-length thing))
+                                    "<vector>"))
+    ((procedure? thing) "#<procedure>")
+    (else
+     (case thing
+       ((()) "()")
+       ((#t) "#t")
+       ((#f) "#f")
+       (else
+        "#<whatsit>")))))
+
+;;; 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 <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))))
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)))))))
+