(define (ww m v) (write (list m v))(newline) v)
(define ex write) ; temp hack
; http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
; (egcd a b) => (x . y) such that ax+by = (gcd a b)
(define (egcd a b) (if (= b 0) '(1 . 0) (let* ((q (quotient a b)) (r (- a (* q b))))
   (if (zero? r) '(0 . 1) (let ((c (egcd b r))) (cons (cdr c) (- (car c) (* (cdr c) q))))))))
; (define (t a b) (let ((e (egcd a b))) (list e (gcd a b) (+ (* a (car e)) (* b (cdr e))))))
(define (veq a b)(let ((la (vector-length a))(lb (vector-length b))
    (rc (lambda (a la b lb) (and (let w ((c la))(let ((C (- c 1))) (or (= c 0)
          (and (= (vector-ref a C)(vector-ref b C)) (w C)))))
       (let w ((c (- lb la))(d la)) (or (= c 0) (and (zero? (vector-ref b d))(w (- c 1)(+ d 1)))))))))
    (if (< la lb) (rc a la b lb) (rc b lb a la))))
(define (tlu t sy) (let r ((v (car t))(n (cddr t))) (let s ((v v)(n n)) (if (null? n) #f
     (if (pair? (car n)) (or (r (car v)(cdar n)) (s (cdr v)(cdr n)))
     (if (eq? (car n) sy) (car v) (s (cdr v)(cdr n))))))))
(define (Do n p) (if (> n 0) (let ((u (- n 1))) (p u) (Do u p))))

; Cartesian product of two streams, streamed.
(define (Cart g h) (lambda () (let* ((G (g))(H (h))(gs (G)))
   (lambda () (let ((n (H)))
      (if (not n) (begin (set! H (h)) (set! gs (G)) (set! n (H))))
      (and gs (cons gs n)))))))
(define (SG->lst sg) (reverse (let ((g (sg))) (let w ((e (g))(lst '()))
          (if e (w (g) (cons e lst)) lst)))))
(define (lst->SG l) (lambda () (let ((c l)) (lambda () (and (pair? c)
     (let ((v (car c))) (begin (set! c (cdr c)) v)))))))

(define (gpa p) (let (
    (m+ (lambda (a b) (let ((x (+ a b))) (if (< x p) x (- x p)))))
    (m- (lambda (a) (if (zero? a) 0 (- p a))))
    (m* (lambda (a b) (modulo (* a b) p)))
    (m/ (let ((rc (lambda (j) (modulo (cdr (egcd p j)) p))))
        (if (< p 10000) (let ((a (make-vector p 0)))
          (lambda (x) (let ((z (vector-ref a x))) (if (zero? z)
            (begin (vector-set! a x (rc x)) (vector-ref a x)) (vector-ref a x))))) rc))))
 (let (
    (p+ (lambda (a b) (let* ((la (vector-length a))(lb (vector-length b))(as (< la lb))
            (s (make-vector (if as lb la))))
      (Do (if as la lb) (lambda (j) (vector-set! s j (m+ (vector-ref a j) (vector-ref b j)))))
      (if as (Do (- lb la) (lambda (j) (vector-set! s (+ la j) (vector-ref b (+ la j)))))
               (Do (- la lb) (lambda (j) (vector-set! s (+ lb j) (vector-ref a (+ lb j)))))) s)))
    (p- (lambda (a) (let* ((la (vector-length a))(s (make-vector la)))
         (Do la (lambda (j) (vector-set! s j (m- (vector-ref a j))))) s))))
 (letrec (
     (trim (lambda (a) (let m ((n (vector-length a))) (if (zero? n) #()
       (if (zero? (vector-ref a (- n 1))) (if (= 1 n) #() (m (- n 1)))
        (let ((r (make-vector n))) (Do n (lambda (j) (vector-set! r j (vector-ref a j)))) r))))))
     (pqr (lambda (N d)
     (let* ((n (trim N))(ln (vector-length n)) (ldm (- (vector-length d) 1)))
        (if (< ln ldm) (cons #() n) (let ((ht (vector-ref d ldm)))
           (if (zero? ht) (ex "divide check")) (let ((htr (m/ ht))
              (nd (make-vector ldm))(q (make-vector (- ln ldm))))
                 (Do ldm (lambda (i) (vector-set! nd i (m* htr (vector-ref d i)))))
                 (Do (- ln ldm) (lambda (j) (let ((t (vector-ref n (+ j ldm)))) (vector-set! q j (m* htr t))
                       (Do ldm (lambda (i) (vector-set! n (+ i j) (m+ (vector-ref n (+ i j))
                           (m- (m* t (vector-ref nd i))))))))))
        (let ((r (make-vector ldm)))
           (Do ldm (lambda (j) (vector-set! r j (vector-ref n j)))) (cons q r))))))))
     (p* (lambda (a b) (let* ((la (vector-length a))(lb (vector-length b)))
        (if (= (+ la lb) 0) #() (let ((p (make-vector (+ la lb -1) 0)))
           (Do la (lambda (i) (Do lb (lambda (j)
              (vector-set! p (+ i j) (m+ (vector-ref p (+ i j))
                 (m* (vector-ref a i) (vector-ref b j)))))))) p)))))
     (pegcd (lambda (a b) (let* ((b (trim b))(qr (pqr a b))(q (car qr))(r (trim (cdr qr))))
        (if (zero? (vector-length r)) (cons #() (vector (m/ (vector-ref b (- (vector-length b) 1)))))
           (let ((c (pegcd b r))) (cons (cdr c) (p+ (car c) (p- (p* (cdr c) q)))))))))
     (pgcmd (lambda (a b) (let* ((A (trim a))(B (trim b))
        (la (vector-length A))(lb (vector-length B)))
           (letrec ((d (lambda (l s) (if (zero? (vector-length s)) l (d s (trim (cdr (pqr l s))))))))
              (let* ((a (if (< la lb) (d B A) (d A B))) (l (vector-length a)))
                 (let ((f (vector-ref a (- l 1)))) (if (> f 1) (let ((r (m/ f)))
                    (Do l (lambda (j) (vector-set! a j (m* r (vector-ref a j))))))) a))))))
     (mexpt (lambda (u p f) ; compute u^p mod f
       (let ((l (vector-length f))) (let pl ((u u)(p p))
          (if (zero? p) (let ((a (make-vector (- l 1) 0))) (vector-set! a 0 1) a)
            (if (even? p) (pl (cdr (pqr (p* u u) f)) (/ p 2))
              (cdr (pqr (p* u (pl u (- p 1))) f))))))))
     (tip (lambda (f)  ; f is vector polynomial, list of coefficients, constant first.
        ;Testing a Polynomial for Irreducibility
        ; From Algorithm 4.69 of HB. of App. Cryptog.
        (let tr ((u #(0 1)) (m (vector-length f))) (or (< m 3)
           (let* ((up (mexpt u p f)) (d (pgcmd f (p+ up (vector 0 (- p 1))))))
             (and (= 1 (vector-length d)) (tr up (- m 2))))))))
     (p->i (lambda (P) (let ((w (vector-length P))) (let m ((n 0))
        (if (= w n) 0 (+ (vector-ref P n) (* p (m (+ n 1)))))))))
     (i->p (lambda (n) (let r ((k n)(s 0)) (if (zero? k) (make-vector s)
         (let ((x (r (quotient k p)(+ s 1)))) (vector-set! x s (remainder k p)) x)))))
     (gap (lambda (m) (let ((N (let P ((z m)) (if (zero? z) 1 (* p (P (- z 1)))))))
        (lambda () (let ((n 0)) (lambda () (and (< n N)
           (let ((x (i->p n))) (set! n (+ n 1)) x))))))))
     (gip (lambda (m) (let ((g ((gap m)))(m (i->p (expt p m))))
        (let r ((l '())) (let ((x (g))) (if x (let ((tp (p+ m x)))
          (if (and (positive? (vector-ref tp 0))(tip tp)) (r (cons tp l)) (r l))) l))))))
     (gfip (lambda (m) (let ((g ((gap m)))(m (i->p (expt p m))))
        (let r () (let ((x (g))) (if x (let ((tp (p+ m x)))
          (if (and (positive? (vector-ref tp 0))(tip tp)) tp (r)))))))))
     (fops (lambda (f) (let (
          (f* (lambda (a b) (cdr (pqr (p* a b) f))))
          (f/ (lambda (a) (trim (cdr (pegcd f a))))))
        (list f* f/ (- (vector-length f) 1))))))

(let ((pl
(cons     (list (list pqr p->i i->p gap gip gfip) (list m+ m- m* m/) (list p+ p- p*)
            (list mexpt trim pgcmd tip pegcd fops p))
  (quote (list (list pqr p->i i->p gap gip gfip) (list m+ m- m* m/) (list p+ p- p*)
            (list mexpt trim pgcmd tip pegcd fops ch))))))
  (lambda (sy) (let ((a (tlu pl sy))) (or a pl))))))))

(define matpak (lambda (zero zer? one fa fs fm fi) (letrec (
   (tf (lambda (x zp) (if (null? x) (zp (list one))
     (if (zer? (caar x)) (let ((z (tf (cdr x) zp)))
     (cons (car z) (cons (car x) (cdr z)))) x))))
   (tr (trnxx zero zer?))
   (sm (lambda (s a)(map (lambda (w) (fm s w)) a)))
   (fms (lambda (a s b) (if (null? a) (sm (fs zero s) b) (if (null? b) a
     (cons (fs (car a) (fm s (car b))) (fms (cdr a) s (cdr b)))))))
   (deter (lambda (a) (letrec ((p #f)
     (tf (lambda (x) (if (null? x) (list one)
       (if (zer? (caar x)) (let ((z (tf (cdr x))))
       (set! p (not p))
       (cons (car z) (cons (car x) (cdr z)))) x)))))
     (let inx ((d one)(a a))
          (if (null? a) (if p (fs zero d) d) (let* (
            (A (tf a))
            (i (fi (caar A)))
            (b (map (lambda (z) (fm z i)) (cdar A))))
         (inx (fm (caar A) d) (map (lambda (x w) (fms x w b))
          (map cdr (cdr A)) (map car (cdr A))))))))))
   (inv (lambda (a nxp) (let ol ((ut (let inx (
      (a (let pad ((x a)(e (list one)))
         (if (null? x) '() (cons (let ap ((z (car x))(ln a))
           (if (null? ln) e (if (null? z) (cons zero (ap z (cdr ln)))
             (cons (car z)(ap (cdr z)(cdr ln))))))
          (pad (cdr x) (cons zero e))))))
      (np nxp))
          (if (null? a) '() (let* (
            (A (tf a np))
            (i (fi (caar A)))
            (b (map (lambda (z) (fm z i)) (cdar A))))
         (cons b (inx (map (lambda (x w) (fms x w b))
          (map cdr (cdr A)) (map car (cdr A)))
          (lambda (w) (np (cons (fs zero (ip w b)) w))))))))))
      (if (null? ut) '() (cons
        (let eg ((top (car ut))(bod (cdr ut))) (if (null? bod) top
          (eg (fms (cdr top) (car top) (car bod))(cdr bod))))
      (ol (cdr ut)))))))
   (ip (lambda (x y)(if (or (null? x) (null? y)) zero (fa (fm (car x)(car y))
        (ip (cdr x)(cdr y))))))
   (mp (lambda (a b)(let ((b (tr b)))
       (map (lambda(x) (map (lambda (y) (ip x y)) b)) a)))))

(list mp inv ip tr deter))))

(define (trnxx zer zer?) (lambda (x) (if (null? x) '()
  (let ((z ((trnxx zer zer?) (cdr x)))) (let m ((u (car x))(v z))
    (if (null? u) (map (lambda (q) (if (null? q) '() (cons zer q))) v)
     (if (null? v) (let y ((q u)) (if (null? q) '() (let ((p (y (cdr q))))
         (if (and (null? p) (zer? (car q))) '()
           (cons (if (zer? (car q)) '() (list (car q))) p)))))
       (cons (cons (car u)(car v)) (m (cdr u)(cdr v))))))))))

(define (DoL n p) (let l ((v '())(n (- n 1))) (if (< n 0) v (l (cons (p n) v)(- n 1)))))

(define (Ffield T q) (let* ((f ((T 'gfip) q))(fops ((T 'fops) f))(i->p (T 'i->p))
               (f* (car fops))(f/ (cadr fops)) (zer (i->p 0)) (p+ (T 'p+)) (p- (T 'p-)))
        (matpak zer (lambda(x) (veq x zer)) (i->p 1) p+ (lambda (a b) (p+ a (p- b))) f* f/)))

(define (rm T q n) (let ((i->p (T 'i->p))(sz (expt (T 'ch) q)))
   (DoL n (lambda (d) (DoL n (lambda (d) (i->p (inexact->exact
          (floor (* (random) sz))))))))))

(define (tmz m) (let ((trz
    (lambda (l) (let rz ((w l)) (or (null? w)
       (and (veq #() (car w)) (rz (cdr w))))))))
  (or (null? m) (let ((c (map car m))(R (map cdr m)))
   (and (veq #(1) (car c))
        (trz (car R))
        (trz (cdr c))
        (tmz (cdr R)))))))

(define (mtt p q n) (let* ((T (gpa p))(m (rm T q n))(mt (Ffield T q))
    (minv (cadr mt))(mmul (car mt))(ic (mmul m (minv m cc))))
    (list (tmz ic) m ic)))

(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 cc (mer "Singular:"))