;;; MATRIX -- Obtained from Andrew Wright. ; 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))) (zebra (first-row-perm 'child) (make-row->func 1 -1) (make-row->func -1 1) 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) (modulo 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 (quotient 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? (modulo 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 (quotient 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? (remainder blen 10)) ; (_-*- (quotient 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) (run-benchmark "matrix" matrix-iters (lambda (result) (equal? result '(((1 1 1 1 1) (1 1 1 1 -1) (1 1 1 -1 1) (1 1 -1 -1 -1) (1 -1 1 -1 -1) (1 -1 -1 1 1)) ((1 1 1 1 1) (1 1 1 1 -1) (1 1 1 -1 1) (1 1 -1 1 -1) (1 -1 1 -1 -1) (1 -1 -1 1 1)) ((1 1 1 1 1) (1 1 1 1 -1) (1 1 1 -1 1) (1 1 -1 1 -1) (1 -1 1 -1 1) (1 -1 -1 1 1)) ((1 1 1 1 1) (1 1 1 1 -1) (1 1 1 -1 1) (1 1 -1 1 1) (1 -1 1 1 -1) (1 -1 -1 -1 1)) ((1 1 1 1 1) (1 1 1 1 -1) (1 1 1 -1 1) (1 1 -1 1 1) (1 -1 1 1 1) (1 -1 -1 -1 -1))))) (lambda (number-of-cols inv-size) (lambda () (really-go number-of-cols inv-size))) 5 5))