;;; Primitives for Siskind's EM benchmarks,
;;; as revised by Will Clinger and Lars Hansen.
;;; Functional code is not allowed to use the primitives whose
;;; names contain an exclamation mark.

;;; Portability.
;;;
;;; A run-benchmark procedure must be supplied separately.
;;; (run-benchmark <string> <thunk> n) reports the name of the
;;; benchmark as <string> and also the time it takes to execute
;;; n iterations of the <thunk>.  If n is omitted, it should
;;; default to 1.  A sample definition of run-benchmark, for
;;; Chez Scheme, is provided in "run-benchmark.chez".
;;;
;;; The rest of this code should run in any R4RS-compatible system
;;; that uses IEEE floating point.  The most common portability
;;; problem is that some systems report an error for (log 0.0).
;;; Siskind's code therefore uses ieee-log instead, which can be
;;; defined in terms of log as here.  We use (exp 1000.0) to
;;; create a floating point infinity; this seems to work in most
;;; systems.

(define ieee-infinity (exp 1000.0))

(define (ieee-log z)
  (if (= z 0.0)
      (- ieee-infinity) ; computed only once during benchmark
      (log z)))

;;; End of code that depends upon IEEE floating point.

;;; The constants are hardwired to be inexact for efficiency.

(begin
 (define make-model vector)
 (define (model-pi model) (vector-ref model 0))
 (define (set-model-pi! model pi) (vector-set! model 0 pi))
 (define (model-mu model) (vector-ref model 1))
 (define (model-sigma model) (vector-ref model 2))
 (define (model-log-pi model) (vector-ref model 3))
 (define (set-model-log-pi! model log-pi) (vector-set! model 3 log-pi))
 (define (model-sigma-inverse model) (vector-ref model 4))
 (define (model-log-determinant-sigma model) (vector-ref model 5))
 (define (set-model-log-determinant-sigma! model log-determinant-sigma)
   (vector-set! model 5 log-determinant-sigma))
 )

(define (void) #f)

(define panic
  (lambda (s)
    "The panic procedure is assigned dynamically."))

(define (hex-string->number s)
 (let loop ((s (string->list s)) (c 0))
  (if (null? s)
      c
      (loop (cdr s) (+ (* 16 c)
                       (if (char-numeric? (car s))
                           (- (char->integer (car s)) (char->integer #\0))
                           (+ (- (char->integer (car s)) (char->integer #\a))
                              10)))))))

;;; The following code is a modified version of code taken from SLIB.
;;; Copyright (C) 1991, 1993 Aubrey Jaffer.
;
;Permission to copy this software, to redistribute it, and to use it
;for any purpose is granted, subject to the following restrictions and
;understandings.
;
;1.  Any copy made of this software must include this copyright notice
;in full.
;
;2.  I have made no warrantee or representation that the operation of
;this software will be error-free, and I am under no obligation to
;provide any services, by way of maintenance, update, or otherwise.
;
;3.  In conjunction with products arising from the use of this
;material, there shall be no use of my name in any advertising,
;promotional, or sales literature without prior written consent in
;each case.

(define *most-positive-fixnum* 65535)

(define (logical:logxor n1 n2)
 (cond ((= n1 n2) 0)
       ((zero? n1) n2)
       ((zero? n2) n1)
       (else (+ (* (logical:logxor (logical:ash-4 n1) (logical:ash-4 n2)) 16)
                (vector-ref (vector-ref logical:boole-xor (modulo n1 16))
                            (modulo n2 16))))))

(define (logical:ash-4 x)
 (if (negative? x) (+ -1 (quotient (+ 1 x) 16)) (quotient x 16)))

(define logical:boole-xor
 '#(#(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15)
    #(1 0 3 2 5 4 7 6 9 8 11 10 13 12 15 14)
    #(2 3 0 1 6 7 4 5 10 11 8 9 14 15 12 13)
    #(3 2 1 0 7 6 5 4 11 10 9 8 15 14 13 12)
    #(4 5 6 7 0 1 2 3 12 13 14 15 8 9 10 11)
    #(5 4 7 6 1 0 3 2 13 12 15 14 9 8 11 10)
    #(6 7 4 5 2 3 0 1 14 15 12 13 10 11 8 9)
    #(7 6 5 4 3 2 1 0 15 14 13 12 11 10 9 8)
    #(8 9 10 11 12 13 14 15 0 1 2 3 4 5 6 7)
    #(9 8 11 10 13 12 15 14 1 0 3 2 5 4 7 6)
    #(10 11 8 9 14 15 12 13 2 3 0 1 6 7 4 5)
    #(11 10 9 8 15 14 13 12 3 2 1 0 7 6 5 4)
    #(12 13 14 15 8 9 10 11 4 5 6 7 0 1 2 3)
    #(13 12 15 14 9 8 11 10 5 4 7 6 1 0 3 2)
    #(14 15 12 13 10 11 8 9 6 7 4 5 2 3 0 1)
    #(15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0)))

(define random:tap 24)

(define random:size 55)

(define (random:size-int l)
 (let ((trial (hex-string->number (make-string l #\f))))
  (if (and (exact? trial) (positive? trial) (>= *most-positive-fixnum* trial))
      l
      (random:size-int (- l 1)))))

(define random:chunk-size (* 4 (random:size-int 8)))

(define random:MASK
 (hex-string->number (make-string (quotient random:chunk-size 4) #\f)))

(define *random-state* '#())

(define (random:initialize)
  (let ((random-strings
         '#("d909ef3e" "fd330ab3" "e33f7843" "76783fbd" "f3675fb3"
            "b54ef879" "0be45590" "a6794679" "0bcd56d3" "fabcdef8"
            "9cbd3efd" "3fd3efcd" "e064ef27" "dddecc08" "34444292"
            "85444454" "4c519210" "c0366273" "54734567" "70abcddc"
            "1bbdac53" "616c5a86" "a982efa9" "105996a0" "5f0cccba"
            "1ea055e1" "fe2acd8d" "1891c1d4" "e6690270" "6912bccc"
            "2678e141" "61222224" "907abcbb" "4ad6829b" "9cdd1404"
            "57798841" "5b892496" "871c9cd1" "d1e67bda" "8b0a3233"
            "578ef23f" "28274ef6" "823ef5ef" "845678c5" "e67890a5"
            "5890abcb" "851fa9ab" "13efa13a" "b12278d6" "daf805ab"
            "a0befc36" "0068a7b5" "e024fd90" "a7b690e2" "27f3571a"
            )))
    (set! *random-state* (make-vector (+ random:size 1) 0))
    (let ((nibbles (quotient random:chunk-size 4)))
      (do ((i 0 (+ i 1)))
          ((= i random:size))
          (vector-set!
           *random-state*
           i
           (hex-string->number
            (substring (vector-ref random-strings i) 0 nibbles)))))))

(random:initialize)

;;; random:chunk returns an integer in the range of
;;; 0 to (- (expt 2 random:chunk-size) 1)
(define (random:chunk v)
 (let* ((p (vector-ref v random:size))
        (ans (logical:logxor
              (vector-ref v (modulo (- p random:tap) random:size))
              (vector-ref v p))))
  (vector-set! v p ans)
  (vector-set! v random:size (modulo (- p 1) random:size))
  ans))

(define (rand)
 (do ((ilen 0 (+ 1 ilen))
      (s random:MASK (+ random:MASK (* (+ 1 random:MASK) s))))
   ((>= s (- *most-positive-fixnum* 1))
    (let ((slop (modulo (+ s (- 1 *most-positive-fixnum*))
                        *most-positive-fixnum*)))
     (let loop ((n ilen) (r (random:chunk *random-state*)))
      (cond ((not (zero? n))
             (loop (+ -1 n)
                   (+ (* r (+ 1 random:MASK)) (random:chunk *random-state*))))
            ((>= r slop) (modulo r *most-positive-fixnum*))
            (else (loop ilen (random:chunk *random-state*)))))))))

;;; End of code taken from SLIB

(define log-math-precision 35.0)

(define minus-infinity (ieee-log 0.0))

(define first car)

(define second cadr)

(define rest cdr)

'(begin
  (define-syntax map-vector
    (syntax-rules
     ()
     ((map-vector ff vv)
      (let* ((f ff)
             (v vv)
             (n (vector-length v))
             (u (make-vector n)))
        (do ((i 0 (+ i 1)))
            ((= i n)
             u)
            (vector-set! u i (f (vector-ref v i))))))))
  
  (define-syntax map-vector-two
    (syntax-rules
     ()
     ((map-vector-two ff vv1 vv2)
      (let* ((f ff)
             (v1 vv1)
             (v2 vv2)
             (n (vector-length v1))
             (u (make-vector n)))
        (do ((i 0 (+ i 1)))
            ((= i n)
             u)
            (vector-set! u i (f (vector-ref v1 i) (vector-ref v2 i))))))))
  
  (define (reduce f l i)
    (cond ((null? l) i)
          ((null? (rest l)) (first l))
          (else (let loop ((l (rest l)) (c (first l)))
                  (if (null? l) c (loop (rest l) (f c (first l))))))))
  
  (define-syntax reduce-vector
    (syntax-rules
     ()
     ((reduce-vector ff vv iinit)
      (let* ((f ff)
             (v vv)
             (init iinit)
             (n (vector-length v)))
        (cond ((zero? n) init)
              ((= n 1) (vector-ref v 0))
              (else (let loop ((i 1) (c (vector-ref v 0)))
                      (if (= i n) c (loop (+ i 1) (f c (vector-ref v i)))))))))))
  )

(begin
 (define (map-vector f v)
   (let* ((n (vector-length v))
          (u (make-vector n)))
     (do ((i 0 (+ i 1)))
         ((= i n)
          u)
         (vector-set! u i (f (vector-ref v i))))))
 
 (define (map-vector-two f v1 v2)
   (let* ((n (vector-length v1))
          (u (make-vector n)))
     (do ((i 0 (+ i 1)))
         ((= i n)
          u)
         (vector-set! u i (f (vector-ref v1 i) (vector-ref v2 i))))))
 
 (define (reduce f l i)
   (cond ((null? l) i)
         ((null? (rest l)) (first l))
         (else (let loop ((l (rest l)) (c (first l)))
                 (if (null? l) c (loop (rest l) (f c (first l))))))))
 
 (define (reduce-vector f v init)
   (let* ((n (vector-length v)))
     (cond ((zero? n) init)
           ((= n 1) (vector-ref v 0))
           (else (let loop ((i 1) (c (vector-ref v 0)))
                   (if (= i n) c (loop (+ i 1) (f c (vector-ref v i)))))))))
 )

(define (every-n p n)
 (let loop ((i 0)) (or (>= i n) (and (p i) (loop (+ i 1))))))

(define (sum f n)
 (let loop ((n (- n 1)) (c 0.0))
  (if (negative? n) c (loop (- n 1) (+ c (f n))))))

'(begin
  (define-syntax dot
    (syntax-rules
     ()
     ((dot uu vv)
      (let* ((u uu)
             (v vv)
             (n (vector-length u)))
        (do ((i 0 (+ i 1))
             (sum 0.0 (+ sum
                         (* (vector-ref u i)
                            (vector-ref v i)))))
            ((= i n) sum))))))
  
  (define-syntax v+
    (syntax-rules
     ()
     ((v+ u v)
      (map-vector-two + u v))))
  
  (define-syntax v-
    (syntax-rules
     ()
     ((v- u v)
      (map-vector-two - u v))))
  )

(begin
 (define (dot u v)
   (let* ((n (vector-length u)))
     (do ((i 0 (+ i 1))
          (sum 0.0 (+ sum
                      (* (vector-ref u i)
                         (vector-ref v i)))))
         ((= i n) sum))))
 
 (define (v+ u v)
   (map-vector-two + u v))
 
 (define (v- u v)
   (map-vector-two - u v))
 )

(define (k*v k u) (map-vector (lambda (x) (* k x)) u))

(define (add-exp e1 e2)
 (let* ((e-max (max e1 e2))
        (e-min (min e1 e2))
        (factor (floor e-min)))
  (if (= e-max minus-infinity)
      minus-infinity
      (if (> (- e-max factor) log-math-precision)
          e-max
          (+ (ieee-log (+ (exp (- e-max factor)) (exp (- e-min factor))))
             factor)))))

(define (map-n f n)
 ;; needs work: To eliminate REVERSE.
 (let loop ((i 0) (c '()))
  (if (< i n) (loop (+ i 1) (cons (f i) c)) (reverse c))))

(define (map-n-vector f n)
 (let ((v (make-vector n)))
  (let loop ((i 0))
   (if (< i n)
       (begin (vector-set! v i (f i))
              (loop (+ i 1)))))
  v))

(define (remove-if-not p l)
 ;; needs work: To eliminate REVERSE.
 (let loop ((l l) (c '()))
  (cond ((null? l) (reverse c))
        ((p (first l)) (loop (rest l) (cons (first l) c)))
        (else (loop (rest l) c)))))

(define (positionv x l)
 (let loop ((l l) (i 0))
  (cond ((null? l) #f)
        ((eqv? x (first l)) i)
        (else (loop (rest l) (+ i 1))))))

(define (make-matrix m n) (map-n-vector (lambda (i) (make-vector n)) m))

(define (make-matrix-initial m n initial)
 (map-n-vector (lambda (i) (make-vector n initial)) m))

(define (matrix-copy m0)
  ; older code
  '(let* ((m (matrix-rows m0))
          (n (matrix-columns m0))
          (m1 (make-matrix m n)))
         (do ((i 0 (+ i 1)))
             ((= i m) m1)
             (do ((j 0 (+ j 1)))
                 ((= j n))
                 (matrix-set! m1 i j (matrix-ref m0 i j)))))
  (map-n-vector
   (lambda (i)
     (let ((mi (vector-ref m0 i)))
       (map-n-vector
        (lambda (j)
          (vector-ref mi j))
        (matrix-columns m0))))
   (matrix-rows m0)))

(define (matrix-rows a) (vector-length a))

(define (matrix-columns a) (vector-length (vector-ref a 0)))

(define (matrix-ref a i j) (vector-ref (vector-ref a i) j))

(define (matrix-set! a i j x) (vector-set! (vector-ref a i) j x))

(define (matrix-row-ref a i) (vector-ref a i))

(define (matrix-column-ref a j) (map-vector (lambda (v) (vector-ref v j)) a))

(define (matrix-row-set! a i v) (vector-set! a i v))

(define (m+ a b) (map-vector-two (lambda (u v) (v+ u v))
                                 a b))

(define (m- a b) (map-vector-two (lambda (u v) (v- u v))
                                 a b))

(define (m*v a v) (map-vector (lambda (u) (dot u v)) a))

(define (transpose a)
 (map-n-vector (lambda (j) (matrix-column-ref a j)) (matrix-columns a)))

(define (outer-product f u v)
 (map-vector (lambda (ui) (map-vector (lambda (vj) (f ui vj)) v)) u))

(define (self-outer-product f v) (outer-product f v v))

(define (m* a b) (outer-product (lambda (u v)
                                  (dot u v))
                                a
                                (transpose b)))

(define (k*m k m)
 (map-vector (lambda (row) (map-vector (lambda (e) (* k e)) row)) m))

(define (determinant a)
 (if (not (= (matrix-rows a) (matrix-columns a)))
     (panic "Can only find determinant of a square matrix"))
 (call-with-current-continuation
  (lambda (return)
   (let* ((n (matrix-rows a))
          (b (make-matrix n n))
          (d 1.0))
    (do ((i 0 (+ i 1))) ((= i n))
     (do ((j 0 (+ j 1))) ((= j n)) (matrix-set! b i j (matrix-ref a i j))))
    (do ((i 0 (+ i 1))) ((= i n))
     ;; partial pivoting reduces rounding errors
     (let ((greatest (abs (matrix-ref b i i)))
           (index i))
      (do ((j (+ i 1) (+ j 1))) ((= j n))
       (let ((x (abs (matrix-ref b j i))))
        (if (> x greatest) (begin (set! index j) (set! greatest x)))))
      (if (= greatest 0.0) (return 0.0))
      (if (not (= index i))
          (let ((v (matrix-row-ref b i)))
           (matrix-row-set! b i (matrix-row-ref b index))
           (matrix-row-set! b index v)
           (set! d (- d))))
      (let ((c (matrix-ref b i i)))
       (set! d (* d c))
       (do ((j i (+ j 1))) ((= j n))
        (matrix-set! b i j (/ (matrix-ref b i j) c)))
       (do ((j (+ i 1) (+ j 1))) ((= j n))
        (let ((e (matrix-ref b j i)))
         (do ((k (+ i 1) (+ k 1))) ((= k n))
          (matrix-set!
           b j k (- (matrix-ref b j k) (* e (matrix-ref b i k))))))))))
    d))))

(define (invert-matrix! a b)
  (if (not (= (matrix-rows a) (matrix-columns a)))
      (panic "Can only invert a square matrix"))
  (let* ((n (matrix-rows a))
         (c (make-matrix n n)))
    (do ((i 0 (+ i 1)))
        ((= i n))
        (do ((j 0 (+ j 1)))
            ((= j n))
            (matrix-set! b i j 0.0)
            (matrix-set! c i j (matrix-ref a i j))))
    (do ((i 0 (+ i 1)))
        ((= i n)) (matrix-set! b i i 1.0))
    (do ((i 0 (+ i 1)))
        ((= i n))
        (if (zero? (matrix-ref c i i))
            (call-with-current-continuation
             (lambda (return)
               (do ((j 0 (+ j 1)))
                   ((= j n))
                   (if (and (> j i) (not (zero? (matrix-ref c j i))))
                       (begin (let ((e (vector-ref c i)))
                                (vector-set! c i (vector-ref c j))
                                (vector-set! c j e))
                              (let ((e (vector-ref b i)))
                                (vector-set! b i (vector-ref b j))
                                (vector-set! b j e))
                              (return (void)))))
               (panic "Matrix is singular"))))
        (let ((d (/ (matrix-ref c i i))))
          (do ((j 0 (+ j 1)))
              ((= j n))
              (matrix-set! c i j (* d (matrix-ref c i j)))
              (matrix-set! b i j (* d (matrix-ref b i j))))
          (do ((k 0 (+ k 1)))
              ((= k n))
              (let ((d (- (matrix-ref c k i))))
                (if (not (= k i))
                    (do ((j 0 (+ j 1)))
                        ((= j n))
                        (matrix-set!
                         c k j (+ (matrix-ref c k j) (* d (matrix-ref c i j))))
                        (matrix-set!
                         b k j (+ (matrix-ref b k j) (* d (matrix-ref b i j))))))))))))

;;; Functional version.

(define (invert-matrix a)
  (let* ((n (matrix-rows a))
         (b (make-matrix-initial n n 0.0)))
    (invert-matrix! a b)
    b))

(define (jacobi! a)
  (if (not (and (= (matrix-rows a) (matrix-columns a))
                (every-n (lambda (i)
                           (every-n (lambda (j)
                                      (= (matrix-ref a i j) (matrix-ref a j i)))
                                    (matrix-rows a)))
                         (matrix-rows a))))
      (panic "Can only compute eigenvalues/eigenvectors of a symmetric matrix"))
  (let* ((n (matrix-rows a))
         (d (make-vector n))
         (v (make-matrix-initial n n 0.0))
         (b (make-vector n))
         (z (make-vector n 0.0)))
    (do ((ip 0 (+ ip 1)))
        ((= ip n))
        (matrix-set! v ip ip 1.0)
        (vector-set! b ip (matrix-ref a ip ip))
        (vector-set! d ip (matrix-ref a ip ip)))
    (let loop ((i 0))
      (if (> i 50) (panic "Too many iterations in JACOBI!"))
      (let ((sm (sum (lambda (ip)
                       (sum (lambda (ir)
                              (let ((iq (+ ip ir 1)))
                                (abs (matrix-ref a ip iq))))
                            (- n ip 1)))
                     (- n 1))))
        (if (not (zero? sm))
            (begin
             (let ((tresh (if (< i 3) (/ (* 0.2 sm) (* n n)) 0.0)))
               (do ((ip 0 (+ ip 1)))
                   ((= ip (- n 1)))
                   (do ((ir 0 (+ ir 1)))
                       ((= ir (- n ip 1)))
                       (let* ((iq (+ ip ir 1))
                              (g (* 100.0 (abs (matrix-ref a ip iq)))))
                         (cond
                          ((and (> i 3)
                                (= (+ (abs (vector-ref d ip)) g)
                                   (abs (vector-ref d ip)))
                                (= (+ (abs (vector-ref d iq)) g)
                                   (abs (vector-ref d iq))))
                           (matrix-set! a ip iq 0.0))
                          ((> (abs (matrix-ref a ip iq)) tresh)
                           (let* ((h (- (vector-ref d iq) (vector-ref d ip)))
                                  (t (if (= (+ (abs h) g) (abs h))
                                         (/ (matrix-ref a ip iq) h)
                                         (let ((theta (/ (* 0.5 h)
                                                         (matrix-ref a ip iq))))
                                           (if (negative? theta)
                                               (- (/ (+ (abs theta)
                                                        (sqrt (+ (* theta theta)
                                                                 1.0)))))
                                               (/ (+ (abs theta)
                                                     (sqrt (+ (* theta theta) 1.0))))))))
                                  (c (/ (sqrt (+ (* t t) 1.0))))
                                  (s (* t c))
                                  (tau (/ s (+ c 1.0)))
                                  (h (* t (matrix-ref a ip iq))))
                             (define (rotate a i j k l)
                               (let ((g (matrix-ref a i j))
                                     (h (matrix-ref a k l)))
                                 (matrix-set! a i j (- g (* s (+ h (* g tau)))))
                                 (matrix-set! a k l (+ h (* s (- g (* h tau)))))))
                             (vector-set! z ip (- (vector-ref z ip) h))
                             (vector-set! z iq (+ (vector-ref z iq) h))
                             (vector-set! d ip (- (vector-ref d ip) h))
                             (vector-set! d iq (+ (vector-ref d iq) h))
                             (matrix-set! a ip iq 0.0)
                             (do ((j 0 (+ j 1)))
                                 ((= j n))
                                 (cond ((< j ip) (rotate a j ip j iq))
                                       ((< ip j iq) (rotate a ip j j iq))
                                       ((< iq j) (rotate a ip j iq j)))
                                 (rotate v j ip j iq)))))))))
             (do ((ip 0 (+ ip 1)))
                 ((= ip n))
                 (vector-set! b ip (+ (vector-ref b ip) (vector-ref z ip)))
                 (vector-set! d ip (vector-ref b ip))
                 (vector-set! z ip 0.0))
             (loop (+ i 1))))))
    (do ((i 0 (+ i 1)))
        ((= i (- n 1)))
        (let ((k i)
              (p (vector-ref d i)))
          (do ((l 0 (+ l 1)))
              ((= l (- n i 1)))
              (let* ((j (+ i l 1)))
                (if (>= (vector-ref d j) p)
                    (begin (set! k j) (set! p (vector-ref d j))))))
          (if (not (= k i))
              (begin (vector-set! d k (vector-ref d i))
                     (vector-set! d i p)
                     (do ((j 0 (+ j 1)))
                         ((= j n))
                         (let ((p (matrix-ref v j i)))
                           (matrix-set! v j i (matrix-ref v j k))
                           (matrix-set! v j k p)))))))
    (list d v)))

;;; Functional version.

(define (jacobi a)
  (jacobi! (matrix-copy a)))

(define (vector->diagonal-matrix v)
  (let* ((n (vector-length v))
         (m (make-matrix-initial n n 0.0)))
    (do ((i 0 (+ i 1)))
        ((= i n)
         m)
        (matrix-set! m i i (vector-ref v i)))))

(define (clip-eigenvalues! a v)
  (let* ((j (jacobi! a))
         (l (first j))
         (e (second j)))
    (do ((k1 0 (+ k1 1)))
        ((= k1 (vector-length a)))
        (let ((a-k1 (vector-ref a k1))
              (e-k1 (vector-ref e k1)))
          (do ((k2 0 (+ k2 1)))
              ((= k2 (vector-length a-k1)))
              (let ((e-k2 (vector-ref e k2))
                    (s 0.0))
                (do ((k 0 (+ k 1)))
                    ((= k (vector-length a)))
                    (set! s (+ s (* (vector-ref e-k1 k)
                                    (max (vector-ref v k) (vector-ref l k))
                                    (vector-ref e-k2 k)))))
                (vector-set! a-k1 k2 s)))))))

;;; Functional version.

(define (clip-eigenvalues a v)
  (let* ((j (jacobi a))
         (l (first j))
         (e (second j)))
    (map-n-vector
     (lambda (k1)
       (let ((a-k1 (vector-ref a k1))
             (e-k1 (vector-ref e k1)))
         (map-n-vector
          (lambda (k2)
            (let ((e-k2 (vector-ref e k2))
                  (s 0.0))
              (do ((k 0 (+ k 1))
                   (s 0.0 (+ s (* (vector-ref e-k1 k)
                                  (max (vector-ref v k) (vector-ref l k))
                                  (vector-ref e-k2 k)))))
                  ((= k (vector-length a))
                   s))))
          (vector-length a-k1))))
     (vector-length a))))

(define (noise epsilon)
  (- (* 2.0 epsilon (/ (exact->inexact (rand)) *most-positive-fixnum*))
     epsilon))