(define exit (let* ((s (cons 0 0)) (z (call-with-current-continuation (lambda (e) (cons s e))))) (if (and (pair? z) (eq? (car z) s)) (cdr z) exit))) (define mer (lambda (x) (let ((r (call-with-current-continuation (lambda(d) (cons #f (lambda(y) (d (cons #t y)))))))) (if (car r) (begin (write (list "Exception:" x (cdr r))) (newline) (exit 0)) (cdr r))))) (define (merge-assim A B comp assim) (let mg ((A A)(B B)) (if (null? A) B (if (null? B) A (let ((t (comp (car A) (car B)))) (cond ((negative? t) (cons (car A) (mg (cdr A) B))) ((positive? t) (cons (car B) (mg A (cdr B)))) (else (cons (assim (car A)(car B)) (mg (cdr A) (cdr B)))))))))) (define (ms comp assim) (lambda (A) (let nm ((ll (map list A))) (if (pair? ll) (if (pair? (cdr ll)) (nm (let srt ((a ll)) (if (null? a) '() (if (null? (cdr a)) a (cons (merge-assim (car a) (cadr a) comp assim) (srt (cddr a))))))) (car ll)) '())))) (define (Do k proc)(let loop ((j 0)) (if (< j k) (begin (proc j) (loop (+ j 1)))))) (define (vg m n) (let* ( (L (* m n)) (s (make-vector L '())) (si (ms - (mer "Twins"))) (f 1/4) (mf (- f)) (term (list (cons L 0))) (ndx (lambda (i j)(+ (* n i) j))) (d (lambda n (cons (car n) (let sd ((n (si (cdr n)))) (if (null? n) term (cons (cons (car n) f) (sd (cdr n))))))))) (vector-set! s (ndx 0 0) (d 0 (ndx 0 1) (ndx 1 0))) (vector-set! s (ndx 0 (- n 1)) (d 0 (ndx 1 (- n 1)) (ndx 0 (- n 2)))) (vector-set! s (ndx (- m 1) 0) (d 0 (ndx (- m 2) 0) (ndx (- m 1) 1))) (vector-set! s (ndx (- m 1)(- n 1)) (d 0 (ndx (- m 2)(- n 1))(ndx (- m 1)(- n 2)))) (Do (- n 2) (lambda (jm) (let ((j (+ jm 1))) (vector-set! s (ndx 0 j) (d mf (ndx 1 j) (ndx 0 (+ j 1)) (ndx 0 (- j 1)))) (vector-set! s (ndx (- m 1) j) (d mf (ndx (- m 2) j) (ndx (- m 1) (+ j 1)) (ndx (- m 1) (- j 1))))))) (Do (- m 2)(lambda (jm) (let ((j (+ jm 1))) (vector-set! s (ndx j 0) (d f (ndx j 1) (ndx (+ j 1) 0) (ndx (- j 1) 0))) (vector-set! s (ndx j (- n 1)) (d f (ndx j (- n 2)) (ndx (+ j 1) (- n 1)) (ndx (- j 1) (- n 1))))))) (Do (- m 2) (lambda (im) (let ((i (+ im 1))) (Do (- n 2) (lambda (jm) (let ((j (+ jm 1))) (vector-set! s (ndx i j) (d 0 (ndx (- i 1) j) (ndx (+ i 1) j) (ndx i (- j 1)) (ndx i (+ j 1)))))))))) (Do L (lambda (x) (if (null? (vector-ref s x)) (write (list x))))) (cons (list m n) s) )) (define (sl a) ((ms (lambda (x y) (let ((s (- (car x) (car y)))) (if (zero? s) (- (cdr x) (cdr y)) s))) (mer "bump")) a)) (define (retire s n)(let* ((s (cdr s))(lc (vector-ref s n))(n2 (vector-length s))) (if (not (number? (car lc))) ((mer "Dead already") "Pooh")) (let qx ((Lc (cdr lc))) (let* ((N (car Lc))(id (car N))) (if (not (= id n2)) (begin (let* ( (nlc (vector-ref s id)) (f (cdr (assq n (cdr nlc)))) (g (cdr (assq id (cdr lc)))) (r (/ (- 1 (* f g))))) (vector-set! s id (cons (* r (+ (car nlc) (* f (car lc)))) (let xl ((lC (cdr lc))(nlC (cdr nlc))) (if (= (caar lC)(caar nlC)) (if (= (caar lC) n2) lC (cons (cons (caar lC) (* r (+ (* f (cdar lC)) (cdar nlC)))) (xl (cdr lC) (cdr nlC)))) (if (< (caar lC)(caar nlC)) (if (= (caar lC) id) (xl (cdr lC) nlC) (cons (cons (caar lC) (* r f (cdar lC))) (xl (cdr lC) nlC))) (if (= (caar nlC) n) (xl lC (cdr nlC)) (cons (cons (caar nlC) (* r (cdar nlC))) (xl lC (cdr nlC)))))))))) (qx (cdr Lc)))))) (vector-set! s n (cons "d" lc)))) (define (ln s er) (let ((c '())(n2 (vector-length (cdr s)))) (if (not (= n2 (apply * (car s)))) (er (list "Wrong Length" n2))) (Do n2 (lambda (j) (let ((x (vector-ref (cdr s) j))) (if (number? (car x)) (let zl ((el (cdr x))(p -1)) (let ((e (caar el))) (if (= e j) (er (cons "Reflexive!" j))) (if (not (= e n2)) (begin (if (>= p e) ((mer "sort") (cons e p))) (set! c (cons (cons j e) c)) (zl (cdr el) e))))) )))) (let ((d (sl c))(e (sl (map (lambda (x) (cons (cdr x)(car x))) c)))) (if (not (equal? d e)) (er (list "non-symmetric" d e)))))) (define (bs s n) (let ((lc (vector-ref (cdr s) n))(n2 (vector-length (cdr s)))) (if (string=? (car lc) "d") (begin (vector-set! (cdr s) n (cons "k" (let sum ((t (cddr lc))) (if (= (caar t) n2) (cadr lc) (+ (* (cdar t) (let* ((id (caar t))(olc (vector-ref (cdr s) id))) (if (string=? "k" (car olc)) (cdr olc) ((mer "Order") olc)))) (sum (cdr t))))))) lc) ((mer "order") lc)))) (define (vector-copy v) (let* ((l (vector-length v))(x (make-vector l))) (Do l (lambda (j) (vector-set! x j (vector-ref v j)))) x)) (define s (vg 10 10)) (define ge (mer "Problem!")) (define (solve s er) (let ((n (vector-length (cdr s))) (b (cons (car s)(vector-copy (cdr s))))) (ln s er) (Do n (lambda (j) (retire b j))) (Do n (lambda (j) (bs b (- (- n 1) j)))) b)) (define r (solve s ge)) ; and to test for symetry of solution; (define (vector-copy v) (let* ((l (vector-length v))(x (make-vector l))) (Do l (lambda (j) (vector-set! x j (vector-ref v j)))) x)) (define (state->float s)(cons (car s)(let* ((s (cdr s)) (nv (vector-copy s)) (l (vector-length s))) (Do l (lambda (j) (vector-set! nv j (let ((lc (vector-ref s j))) (cons (exact->inexact (car lc)) (let cl ((x (cdr lc))) (if (= (caar x) l) x (cons (cons (caar x) (exact->inexact (cdar x))) (cl (cdr x)))))))))) nv))) (define (flipH a)(let* ((ds (car a))(m (car ds))(n (cadr ds)) (ndx (lambda (i j)(+ (* n i) j)))(b (make-vector (* m n)))) (Do m (lambda (j) (Do n (lambda (k) (vector-set! b (ndx j k) (vector-ref (cdr a) (ndx j (- (- n 1) k)))))))) (cons ds b))) (define (flipV a)(let* ((ds (car a))(m (car ds))(n (cadr ds)) (ndx (lambda (i j)(+ (* n i) j)))(b (make-vector (* m n)))) (Do m (lambda (j) (Do n (lambda (k) (vector-set! b (ndx j k) (vector-ref (cdr a) (ndx (- (- m 1) j) k))))))) (cons ds b))) (define (flipD a)(let* ((ds (car a))(m (car ds))(n (cadr ds))(b (make-vector (* m n)))) (Do m (lambda (j) (Do n (lambda (k) (vector-set! b (+ (* m k) j) (vector-ref (cdr a) (+ k (* n j)))))))) (cons (cons n (cons m (cddr ds))) b))) (define (cmpv a b)(let ((L (vector-length (cdr a)))(c 0)) (Do L (lambda (j) (set! c (max c (abs (- (cdr (vector-ref (cdr a) j)) (cdr (vector-ref (cdr b) j)))))))) c)) (define (negv x)(let* ((a (cdr x))(L (vector-length a))(b (make-vector L))) (Do L (lambda (j) (vector-set! b j (let ((x (vector-ref a j))) (cons (car x)(- (cdr x))))))) (cons (car x) b))) ; r still holds the exact rational solution. (cmpv r (flipH r)) (cmpv r (flipV r)) (cmpv r (negv (flipD r))) (cmpv r (solve (state->float s) ge)) ; < (expt 2 -50) (cmpv r (negv (flipD (solve (vg (cadar s) (caar s)) ge)))) (define s (vg 3 14)) (define r (solve s ge)) (cmpv r (flipH r)) (cmpv r (flipV r)) (cmpv r (solve (state->float s) ge)) ; < (expt 2 -50) (cmpv r (negv (flipD (solve (vg (cadar s) (caar s)) ge)))) (solve (vg 4 25) ge) ; 160 sec (solve (vg 25 4) ge) ; 6 sec (solve (state->float (vg 4 25)) ge) (solve (state->float (vg 25 4)) ge) (cmpv (negv (solve (state->float (vg 4 25)) ge)) (flipD (solve (state->float (vg 25 4)) ge)))