cyclone/test-matrix.scm
Justin Ethier 0113982e87 WIP
2018-10-29 13:11:32 -04:00

881 lines
31 KiB
Scheme

;;; MATRIX -- Obtained from Andrew Wright.
(import (scheme base) (scheme read) (scheme write) (scheme time))
;;;; We need R6RS div and mod for this benchmark.
;
;(define (div x y)
; (cond ((and (exact-integer? x)
; (exact-integer? y)
; (>= x 0))
; (quotient x y))
; ((< y 0)
; ;; x < 0, y < 0
; (let* ((q (quotient x y))
; (r (- x (* q y))))
; (if (= r 0)
; q
; (+ q 1))))
; (else
; ;; x < 0, y > 0
; (let* ((q (quotient x y))
; (r (- x (* q y))))
; (if (= r 0)
; q
; (- q 1))))))
;
;(define (mod x y)
; (cond ((and (exact-integer? x)
; (exact-integer? y)
; (>= x 0))
; (remainder x y))
; ((< y 0)
; ;; x < 0, y < 0
; (let* ((q (quotient x y))
; (r (- x (* q y))))
; (if (= r 0)
; 0
; (- r y))))
; (else
; ;; x < 0, y > 0
; (let* ((q (quotient x y))
; (r (- x (* q y))))
; (if (= r 0)
; 0
; (+ r y))))))
;
;;; Chez-Scheme compatibility stuff:
;
;(define (chez-box x) (cons x '()))
;(define (chez-unbox x) (car x))
;(define (chez-set-box! x y) (set-car! x y))
;
;;; Test that a matrix with entries in {+1, -1} is maximal among the matricies
;;; obtainable by
;;; re-ordering the rows
;;; re-ordering the columns
;;; negating any subset of the columns
;;; negating any subset of the rows
;;; Where we compare two matricies by lexicographically comparing the first row,
;;; then the next to last, etc., and we compare a row by lexicographically
;;; comparing the first entry, the second entry, etc., and we compare two
;;; entries by +1 > -1.
;;; Note, this scheme obeys the useful fact that if (append mat1 mat2) is
;;; maximal, then so is mat1. Thus, we can build up maximal matricies
;;; row by row.
;;;
;;; Once you have chosen the row re-ordering so that you know which row goes
;;; last, the set of columns to negate is fixed (since the last row must be
;;; all +1's).
;;;
;;; Note, the column ordering is really totally determined as follows:
;;; all columns for which the second row is +1 must come before all
;;; columns for which the second row is -1.
;;; among columns for which the second row is +1, all columns for which
;;; the third row is +1 come before those for which the third is
;;; -1, and similarly for columns in which the second row is -1.
;;; etc
;;; Thus, each succeeding row sorts columns withing refinings equivalence
;;; classes.
;;;
;;; Maximal? assumes that mat has atleast one row, and that the first row
;;; is all +1's.
;(define maximal?
; (lambda (mat)
; (let pick-first-row
; ((first-row-perm
; (gen-perms mat)))
; (if first-row-perm
; (and (zunda first-row-perm mat)
; (pick-first-row (first-row-perm 'brother)))
; #t))))
(define zunda
(lambda (first-row-perm mat)
(let* ((first-row
(first-row-perm 'now))
(number-of-cols
(length first-row))
(make-row->func
(lambda (if-equal if-different)
(lambda (row)
(let ((vec
(make-vector number-of-cols)))
(do ((i 0 (+ i 1))
(first first-row
(cdr first))
(row row
(cdr row)))
((= i number-of-cols))
(vector-set! vec
i
(if (= (car first) (car row))
if-equal
if-different)))
(lambda (i)
(vector-ref vec i))))))
(mat
(cdr mat)))
;(make-row->func 1 -1)
(zebra (first-row-perm 'child)
(make-row->func 1 -1)
(make-row->func -1 1)
mat
number-of-cols))))
;; TODO: with this test code, why is the fast-eq inlined????
(write (zunda 1 -1))
(define zebra
(lambda (row-perm row->func+ row->func- mat number-of-cols)
(write (list row-perm row->func+ row->func- mat number-of-cols))))
;(define zebra
; (lambda (row-perm row->func+ row->func- mat number-of-cols)
; (let _-*-
; ((row-perm
; row-perm)
; (mat
; mat)
; (partitions
; (list (miota number-of-cols))))
; (or (not row-perm)
; (and
; (zulu (car mat)
; (row->func+ (row-perm 'now))
; partitions
; (lambda (new-partitions)
; (_-*- (row-perm 'child)
; (cdr mat)
; new-partitions)))
; (zulu (car mat)
; (row->func- (row-perm 'now))
; partitions
; (lambda (new-partitions)
; (_-*- (row-perm 'child)
; (cdr mat)
; new-partitions)))
; (let ((new-row-perm
; (row-perm 'brother)))
; (or (not new-row-perm)
; (_-*- new-row-perm
; mat
; partitions))))))))
;
;
;(define zulu
; (let ((cons-if-not-null
; (lambda (lhs rhs)
; (if (null? lhs)
; rhs
; (cons lhs rhs)))))
; (lambda (old-row new-row-func partitions equal-cont)
; (let _-*-
; ((p-in
; partitions)
; (old-row
; old-row)
; (rev-p-out
; '()))
; (let _-split-
; ((partition
; (car p-in))
; (old-row
; old-row)
; (plus
; '())
; (minus
; '()))
; (if (null? partition)
; (let _-minus-
; ((old-row
; old-row)
; (m
; minus))
; (if (null? m)
; (let ((rev-p-out
; (cons-if-not-null
; minus
; (cons-if-not-null
; plus
; rev-p-out)))
; (p-in
; (cdr p-in)))
; (if (null? p-in)
; (equal-cont (reverse rev-p-out))
; (_-*- p-in old-row rev-p-out)))
; (or (= 1 (car old-row))
; (_-minus- (cdr old-row)
; (cdr m)))))
; (let ((next
; (car partition)))
; (case (new-row-func next)
; ((1)
; (and (= 1 (car old-row))
; (_-split- (cdr partition)
; (cdr old-row)
; (cons next plus)
; minus)))
; ((-1)
; (_-split- (cdr partition)
; old-row
; plus
; (cons next minus)))))))))))
;
;(define all?
; (lambda (ok? lst)
; (let _-*-
; ((lst
; lst))
; (or (null? lst)
; (and (ok? (car lst))
; (_-*- (cdr lst)))))))
;
;(define gen-perms
; (lambda (objects)
; (let _-*-
; ((zulu-future
; objects)
; (past
; '()))
; (if (null? zulu-future)
; #f
; (lambda (msg)
; (case msg
; ((now)
; (car zulu-future))
; ((brother)
; (_-*- (cdr zulu-future)
; (cons (car zulu-future)
; past)))
; ((child)
; (gen-perms
; (fold past cons (cdr zulu-future))))
; ((puke)
; (cons (car zulu-future)
; (fold past cons (cdr zulu-future))))
; (else
; (error 'gen-perms "Bad msg: ~a" msg))))))))
;
;(define fold
; (lambda (lst folder state)
; (let _-*-
; ((lst
; lst)
; (state
; state))
; (if (null? lst)
; state
; (_-*- (cdr lst)
; (folder (car lst)
; state))))))
;
;(define miota
; (lambda (len)
; (let _-*-
; ((i 0))
; (if (= i len)
; '()
; (cons i
; (_-*- (+ i 1)))))))
;
;(define proc->vector
; (lambda (size proc)
; (let ((res
; (make-vector size)))
; (do ((i 0
; (+ i 1)))
; ((= i size))
; (vector-set! res
; i
; (proc i)))
; res)))
;
;;; Given a prime number P, return a procedure which, given a `maker' procedure,
;;; calls it on the operations for the field Z/PZ.
;(define make-modular
; (lambda (modulus)
; (let* ((reduce
; (lambda (x)
; (mod x modulus)))
; (coef-zero?
; (lambda (x)
; (zero? (reduce x))))
; (coef-+
; (lambda (x y)
; (reduce (+ x y))))
; (coef-negate
; (lambda (x)
; (reduce (- x))))
; (coef-*
; (lambda (x y)
; (reduce (* x y))))
; (coef-recip
; (let ((inverses
; (proc->vector (- modulus 1)
; (lambda (i)
; (extended-gcd (+ i 1)
; modulus
; (lambda (gcd inverse ignore)
; inverse))))))
; ;; Coef-recip.
; (lambda (x)
; (let ((x
; (reduce x)))
; (vector-ref inverses (- x 1)))))))
; (lambda (maker)
; (maker 0;; coef-zero
; 1;; coef-one
; coef-zero?
; coef-+
; coef-negate
; coef-*
; coef-recip)))))
;
;;; Extended Euclidean algorithm.
;;; (extended-gcd a b cont) computes the gcd of a and b, and expresses it
;;; as a linear combination of a and b. It returns calling cont via
;;; (cont gcd a-coef b-coef)
;;; where gcd is the GCD and is equal to a-coef * a + b-coef * b.
;(define extended-gcd
; (let ((n->sgn/abs
; (lambda (x cont)
; (if (>= x 0)
; (cont 1 x)
; (cons -1 (- x))))))
; (lambda (a b cont)
; (n->sgn/abs a
; (lambda (p-a p)
; (n->sgn/abs b
; (lambda (q-b q)
; (let _-*-
; ((p
; p)
; (p-a
; p-a)
; (p-b
; 0)
; (q
; q)
; (q-a
; 0)
; (q-b
; q-b))
; (if (zero? q)
; (cont p p-a p-b)
; (let ((mult
; (div p q)))
; (_-*- q
; q-a
; q-b
; (- p (* mult q))
; (- p-a (* mult q-a))
; (- p-b (* mult q-b)))))))))))))
;
;;; Given elements and operations on the base field, return a procedure which
;;; computes the row-reduced version of a matrix over that field. The result
;;; is a list of rows where the first non-zero entry in each row is a 1 (in
;;; the coefficient field) and occurs to the right of all the leading non-zero
;;; entries of previous rows. In particular, the number of rows is the rank
;;; of the original matrix, and they have the same row-space.
;;; The items related to the base field which are needed are:
;;; coef-zero additive identity
;;; coef-one multiplicative identity
;;; coef-zero? test for additive identity
;;; coef-+ addition (two args)
;;; coef-negate additive inverse
;;; coef-* multiplication (two args)
;;; coef-recip multiplicative inverse
;;; Note, matricies are stored as lists of rows (i.e., lists of lists).
;(define make-row-reduce
; (lambda (coef-zero coef-one coef-zero? coef-+ coef-negate coef-* coef-recip)
; (lambda (mat)
; (let _-*-
; ((mat
; mat))
; (if (or (null? mat)
; (null? (car mat)))
; '()
; (let _-**-
; ((in
; mat)
; (out
; '()))
; (if (null? in)
; (map
; (lambda (x)
; (cons coef-zero x))
; (_-*- out))
; (let* ((prow
; (car in))
; (pivot
; (car prow))
; (prest
; (cdr prow))
; (in
; (cdr in)))
; (if (coef-zero? pivot)
; (_-**- in
; (cons prest out))
; (let ((zap-row
; (map
; (let ((mult
; (coef-recip pivot)))
; (lambda (x)
; (coef-* mult x)))
; prest)))
; (cons (cons coef-one zap-row)
; (map
; (lambda (x)
; (cons coef-zero x))
; (_-*-
; (fold in
; (lambda (row mat)
; (cons
; (let ((first-col
; (car row))
; (rest-row
; (cdr row)))
; (if (coef-zero? first-col)
; rest-row
; (map
; (let ((mult
; (coef-negate first-col)))
; (lambda (f z)
; (coef-+ f
; (coef-* mult z))))
; rest-row
; zap-row)))
; mat))
; out))))))))))))))
;
;
;;; Given elements and operations on the base field, return a procedure which
;;; when given a matrix and a vector tests to see if the vector is in the
;;; row-space of the matrix. This returned function is curried.
;;; The items related to the base field which are needed are:
;;; coef-zero additive identity
;;; coef-one multiplicative identity
;;; coef-zero? test for additive identity
;;; coef-+ addition (two args)
;;; coef-negate additive inverse
;;; coef-* multiplication (two args)
;;; coef-recip multiplicative inverse
;;; Note, matricies are stored as lists of rows (i.e., lists of lists).
;(define make-in-row-space?
; (lambda (coef-zero coef-one coef-zero? coef-+ coef-negate coef-* coef-recip)
; (let ((row-reduce
; (make-row-reduce coef-zero
; coef-one
; coef-zero?
; coef-+
; coef-negate
; coef-*
; coef-recip)))
; (lambda (mat)
; (let ((mat
; (row-reduce mat)))
; (lambda (row)
; (let _-*-
; ((row
; row)
; (mat
; mat))
; (if (null? row)
; #t
; (let ((r-first
; (car row))
; (r-rest
; (cdr row)))
; (cond ((coef-zero? r-first)
; (_-*- r-rest
; (map cdr
; (if (or (null? mat)
; (coef-zero? (caar mat)))
; mat
; (cdr mat)))))
; ((null? mat)
; #f)
; (else
; (let* ((zap-row
; (car mat))
; (z-first
; (car zap-row))
; (z-rest
; (cdr zap-row))
; (mat
; (cdr mat)))
; (if (coef-zero? z-first)
; #f
; (_-*-
; (map
; (let ((mult
; (coef-negate r-first)))
; (lambda (r z)
; (coef-+ r
; (coef-* mult z))))
; r-rest
; z-rest)
; (map cdr mat)))))))))))))))
;
;
;;; Given a prime number, return a procedure which takes integer matricies
;;; and returns their row-reduced form, modulo the prime.
;(define make-modular-row-reduce
; (lambda (modulus)
; ((make-modular modulus)
; make-row-reduce)))
;
;
;(define make-modular-in-row-space?
; (lambda (modulus)
; ((make-modular modulus)
; make-in-row-space?)))
;
;
;
;;; Usual utilities.
;
;
;
;;; Given a bound, find a prime greater than the bound.
;(define find-prime
; (lambda (bound)
; (let* ((primes
; (list 2))
; (last
; (chez-box primes))
; (is-next-prime?
; (lambda (trial)
; (let _-*-
; ((primes
; primes))
; (or (null? primes)
; (let ((p
; (car primes)))
; (or (< trial (* p p))
; (and (not (zero? (mod trial p)))
; (_-*- (cdr primes))))))))))
; (if (> 2 bound)
; 2
; (let _-*-
; ((trial
; 3))
; (if (is-next-prime? trial)
; (let ((entry
; (list trial)))
; (set-cdr! (chez-unbox last) entry)
; (chez-set-box! last entry)
; (if (> trial bound)
; trial
; (_-*- (+ trial 2))))
; (_-*- (+ trial 2))))))))
;
;;; Given the size of a square matrix consisting only of +1's and -1's,
;;; return an upper bound on the determinant.
;(define det-upper-bound
; (lambda (size)
; (let ((main-part
; (expt size
; (div size 2))))
; (if (even? size)
; main-part
; (* main-part
; (do ((i 0 (+ i 1)))
; ((>= (* i i) size)
; i)))))))
;
;;; Fold over all maximal matrices.
;(define go
; (lambda (number-of-cols inv-size folder state)
; (let* ((in-row-space?
; (make-modular-in-row-space?
; (find-prime
; (det-upper-bound inv-size))))
; (make-tester
; (lambda (mat)
; (let ((tests
; (let ((old-mat
; (cdr mat))
; (new-row
; (car mat)))
; (fold-over-subs-of-size old-mat
; (- inv-size 2)
; (lambda (sub tests)
; (cons
; (in-row-space?
; (cons new-row sub))
; tests))
; '()))))
; (lambda (row)
; (let _-*-
; ((tests
; tests))
; (and (not (null? tests))
; (or ((car tests) row)
; (_-*- (cdr tests)))))))))
; (all-rows;; all rows starting with +1 in decreasing order
; (fold
; (fold-over-rows (- number-of-cols 1)
; cons
; '())
; (lambda (row rows)
; (cons (cons 1 row)
; rows))
; '())))
; (let _-*-
; ((number-of-rows
; 1)
; (rev-mat
; (list
; (car all-rows)))
; (possible-future
; (cdr all-rows))
; (state
; state))
; (let ((zulu-future
; (remove-in-order
; (if (< number-of-rows inv-size)
; (in-row-space? rev-mat)
; (make-tester rev-mat))
; possible-future)))
; (if (null? zulu-future)
; (folder (reverse rev-mat)
; state)
; (let _-**-
; ((zulu-future
; zulu-future)
; (state
; state))
; (if (null? zulu-future)
; state
; (let ((rest-of-future
; (cdr zulu-future)))
; (_-**- rest-of-future
; (let* ((first
; (car zulu-future))
; (new-rev-mat
; (cons first rev-mat)))
; (if (maximal? (reverse new-rev-mat))
; (_-*- (+ number-of-rows 1)
; new-rev-mat
; rest-of-future
; state)
; state))))))))))))
;
;(define go-folder
; (lambda (mat bsize.blen.blist)
; (let ((bsize
; (car bsize.blen.blist))
; (size
; (length mat)))
; (if (< size bsize)
; bsize.blen.blist
; (let ((blen
; (cadr bsize.blen.blist))
; (blist
; (cddr bsize.blen.blist)))
; (if (= size bsize)
; (let ((blen
; (+ blen 1)))
; ;; (if
; ;; (let _-*-
; ;; ((blen
; ;; blen))
; ;; (or (< blen 10)
; ;; (and (zero? (mod blen 10))
; ;; (_-*- (div blen 10)))))
; ;;
; ;; (begin
; ;; (display blen)
; ;; (display " of size ")
; ;; (display bsize)
; ;; (newline)))
;
; (cons bsize
; (cons blen
; (cond ((< blen 3000)
; (cons mat blist))
; ((= blen 3000)
; (cons "..." blist))
; (else
; blist)))))
; ;; (begin
; ;; (newline)
; ;; (display "First of size ")
; ;; (display size)
; ;; (display ":")
; ;; (newline)
; ;; (for-each
; ;; (lambda (row)
; ;; (display " ")
; ;; (for-each
; ;; (lambda (e)
; ;; (case e
; ;; ((1)
; ;; (display " 1"))
; ;; ((-1)
; ;; (display " -1"))))
; ;; row)
; ;; (newline))
; ;; mat)
;
; (list size 1 mat)))))))
;
;(define really-go
; (lambda (number-of-cols inv-size)
; (cddr
; (go number-of-cols
; inv-size
; go-folder
; (list -1 -1)))))
;
;(define remove-in-order
; (lambda (remove? lst)
; (reverse
; (fold lst
; (lambda (e lst)
; (if (remove? e)
; lst
; (cons e lst)))
; '()))))
;
;;; The first fold-over-rows is slower than the second one, but folds
;;; over rows in lexical order (large to small).
;(define fold-over-rows
; (lambda (number-of-cols folder state)
; (if (zero? number-of-cols)
; (folder '()
; state)
; (fold-over-rows (- number-of-cols 1)
; (lambda (tail state)
; (folder (cons -1 tail)
; state))
; (fold-over-rows (- number-of-cols 1)
; (lambda (tail state)
; (folder (cons 1 tail)
; state))
; state)))))
;
;;; Fold over subsets of a given size.
;(define fold-over-subs-of-size
; (lambda (universe size folder state)
; (let ((usize
; (length universe)))
; (if (< usize size)
; state
; (let _-*-
; ((size
; size)
; (universe
; universe)
; (folder
; folder)
; (csize
; (- usize size))
; (state
; state))
; (cond ((zero? csize)
; (folder universe state))
; ((zero? size)
; (folder '() state))
; (else
; (let ((first-u
; (car universe))
; (rest-u
; (cdr universe)))
; (_-*- size
; rest-u
; folder
; (- csize 1)
; (_-*- (- size 1)
; rest-u
; (lambda (tail state)
; (folder (cons first-u tail)
; state))
; csize
; state))))))))))
;
;(define (main)
; (let* ((count (read))
; (input1 (read))
; (input2 (read))
; (output (read))
; (s3 (number->string count))
; (s2 (number->string input2))
; (s1 (number->string input1))
; (name "matrix"))
; (run-r7rs-benchmark
; (string-append name ":" s1 ":" s2 ":" s3)
; count
; (lambda () (really-go (hide count input1) (hide count input2)))
; (lambda (result) (equal? result output)))))
;
;;;; The following code is appended to all benchmarks.
;
;;;; Given an integer and an object, returns the object
;;;; without making it too easy for compilers to tell
;;;; the object will be returned.
;
;(define (hide r x)
; (call-with-values
; (lambda ()
; (values (vector values (lambda (x) x))
; (if (< r 100) 0 1)))
; (lambda (v i)
; ((vector-ref v i) x))))
;
;;;; Given the name of a benchmark,
;;;; the number of times it should be executed,
;;;; a thunk that runs the benchmark once,
;;;; and a unary predicate that is true of the
;;;; correct results the thunk may return,
;;;; runs the benchmark for the number of specified iterations.
;
;(define (run-r7rs-benchmark name count thunk ok?)
;
; ;; Rounds to thousandths.
; (define (rounded x)
; (/ (round (* 1000 x)) 1000))
;
; (display "Running ")
; (display name)
; (newline)
; (flush-output-port (current-output-port))
; (let* ((j/s (jiffies-per-second))
; (t0 (current-second))
; (j0 (current-jiffy)))
; (let loop ((i 0)
; (result #f))
; (cond ((< i count)
; (loop (+ i 1) (thunk)))
; ((ok? result)
; (let* ((j1 (current-jiffy))
; (t1 (current-second))
; (jifs (- j1 j0))
; (secs (inexact (/ jifs j/s)))
; (secs2 (rounded (- t1 t0))))
; (display "Elapsed time: ")
; (write secs)
; (display " seconds (")
; (write secs2)
; (display ") for ")
; (display name)
; (newline)
; (display "+!CSVLINE!+")
; (display (this-scheme-implementation-name))
; (display ",")
; (display name)
; (display ",")
; (display secs)
; (newline)
; (flush-output-port (current-output-port)))
; result)
; (else
; (display "ERROR: returned incorrect result: ")
; (write result)
; (newline)
; (flush-output-port (current-output-port))
; result)))))
;(define (this-scheme-implementation-name)
; (string-append "cyclone-" (Cyc-version)))
;(main)