(define (ww m v) (write (list m v))(newline) v)
(define ex write) ; temp hack
(define (Do n p) (if (> n 0) (let ((u (- n 1))) (p u) (Do u p))))
; 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 (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))
          (if (= p 1) u (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)))))
          (fexpt (lambda (a p)(mexpt a p f))))
        (list f* f/ fexpt (- (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 (ffp p q) (let* ((T (gpa p))(f ((T 'gfip) q))(h ((T 'fops) f))
     (g (((T 'gap) q)))(p0 (g)))
  (list p0 (g) #(0 1) (T 'p+)(T 'p-)(car h)(cadr h))))

(define (tiv p q)(let* ((T (gpa p))
       (g (((T 'gap) q)))(f ((T 'gfip) q))(GF ((T 'fops) f))(f* (car GF))(f/ (cadr GF)))
       (g) ; dont invert 0!
    (let w () (let ((e (g))) (or (not e) (and (veq #(1) (f* e (f/ e))) (w)))))))

(define T (gpa 5))
(define pqr (T 'pqr))
(define p* (T 'p*))
(define f ((T 'gfip) 4))
; f ; => #5(2 0 0 0 1)
(define FS ((T 'fops) f))
(define f* (car FS))
(define f/ (cadr FS))
(f* #(2 0 4 3) (f/  #(2 0 4 3))) ; => #(1 0 0 0)
(f/ (f/ #(1 2 4 3))) ; => #4(1 2 4 3)

(((gpa 3) 'gfip) 4) ; => #5(2 1 0 0 1)
(((gpa 241) 'gfip) 6) ; => #7(7 0 0 0 0 0 1)

(define T (gpa 5))
(define pegcd (T 'pegcd))
(define f #(2 0 0 0 1))
(define z #(0 2))
(pegcd f z) ; => (#(3) . #(0 0 0 1))
(define p* (T 'p*))
(define p+ (T 'p+))
(p+ (p* #(3) f) (p* #(0 0 0 1) z)) ; => #(1 0 0 0 0) (good!!)
(define fopsp ((T 'fops) f))
(define f* (car fopsp))
(define f/ (cadr fopsp))

(define T (gpa 5))
(define p* (T 'p*))
(define d #(1 3 4))
(define pqr (T 'pqr))
(define pegcd (T 'pegcd))
(pegcd (p* d #(2 1 3)) d) ; => (#0() . #1(4))

(let ((T (gpa 5))) ((cadr ((T 'fops) ((T 'gfip) 4))) #(2 4 3 1 3))) ;=> #4(4 0 0 1)

(define T (gpa 5))
(define FS ((T 'fops) #5(2 0 0 0 1)))
(define f/ (cadr FS))
(f/ #(1 2 4 3)) ; => #2(2 1)

(define T (gpa 5))
(let* ((f*/ ((T 'fops) ((T 'gfip) 4)))(f* (car f*/))(f/ (cadr f*/)))
   (f/ (f/ #(0 2 0 0))))

(map ((gpa 3) 'tip) (list
  #(1 0 1) #(1 0 2) #(1 1 1) #(1 1 2) #(1 2 1) #(1 2 2))) ;=> (#t #f #f #t #f #t)

(((gpa 3) 'pgcmd) #(1 1 1 1) #(1 1)) ; => #(1 1)

(((gpa 3) 'pgcmd) #(1 1 1 1) #(1 1)) ; => #(1 1)

(((gpa 5) 'pegcd) #(1) #(2)) ;=> (#0() . #1(3))
(((gpa 5) 'pgcmd) #(1) #(2)) ; => #(1)
(((gpa 5) 'pegcd) #(2 0 0 0 1) #(2)) ;=> (#() . #(3))
(((gpa 5) 'pgcmd) #(2 0 0 0 1) #(2)) ; => #(1)

(((gpa 3) 'gip) 4) ; => (#5(2 1 2 2 1) #5(1 2 1 2 1) #5(2 1 1 2 1) #5(1 0 1 2 1) #5(1 1 0 2 1) #5(2 0 0 2 1) #5(2 2 2 1) #5(2 2 1) #5(1) #5(1 0 1) #5(1 2 0 1) #5(2 0 0 1) #5(2 0 2 0 1) #5(1 2 1 0 1) #5(1 1 1 0 1) #5(2 0 1 0 1) #5(2 2 0 0 1) #5(2 1 0 0 1))

(map ((gpa 2) 'tip) (list
   #(1 0 1) ; <- reducible
   #(1 1 1) #(1 1 0 1) #(1 1 0 0 1) #(1 0 1 0 0 1) #(1 1 0 0 0 0 1) #(1 1 0 0 0 1 1 0 1)
   #(1 0 0 0 1 0 0 0 0 1) #(1 0 0 1 0 0 0 0 0 0 1) #(1 0 1 0 0 0 0 0 0 0 0 1) ; <- 11
   #( 1 0 0 1 1 0 0 1 0 0 0 0 1)
   )) ; all but the first should be irreducible. (from HB App. Crypt.)

(let ((T (gpa 13))) ((T 'tip) ((T 'p*) #(2 5 4) #(8 3 7))))

(map ((gpa 13) 'trim) (list #(2 3 4) #( 3 2 0) #(5 6 0 0) #(0 0) #(0) #()))
; => (#3(2 3 4) #2(3 2) #2(5 6) #0() #0() #0())

(let ((tpqr (lambda (n d) (let* ((T (gpa 13))(p+ (T 'p+))(p- (T 'p-))(p* (T 'p*))
     (ln (vector-length n))(pn (make-vector ln)))
     (Do ln (lambda (j) (vector-set! pn j (vector-ref n j))))
     (let* ((P ((T 'pqr) n d))(q (car P))(r (cdr P)))
   (list (p+ (p* q d) r) pn q r))))))
   (list (tpqr #(9 9 9 9 9) #(1 1)) (tpqr #(3 3 3) #(1 1)) (tpqr #(1 1) #(3 3 3))))

(((gpa 11) 'p*) #(1 1 1) #(1 1 1)) ; => #5(1 2 3 2 1)

(map ((gpa 11) 'pqr)
   (list #(3 4 7 1 2 2) #(3 4 7 1 2 2) #(3 4 7 1 2 2) #(3 4 7 1 2 2) #(3 4 7 1 2 2))
   (list  #(1 0 1) #(2 0 2) #(1 0 2) #(2 0 1) #(2 0 1 5)))
 ; => ((#4(5 10 2) . #2(9 5)) (#4(8 5 1) . #2(9 5)) (#4(3 0 1) . #2(0 4)) (#4(3 8 2) . #2(8 10)) (#3(7 10 7) . #3(0 6 8)))
(((gpa 17) 'pqr) #(3 4 7 1 2) #(2 0 8 1 5)) ; => (#1(14) . #4(9 4 14 4))

(define (tntip p q r)
(let* ((T (gpa p)) (p+ (T 'p+))(p* (T 'p*))(tip (T 'tip))
       (i->p (T 'i->p))(gap (T 'gap))(f (Cart (gap q) (gap r)))(F (f))
       (monq (i->p (expt p q)))(monr (i->p (expt p r))))
      (let r () (let ((g (F))) (if g (let (
          (x (p+ monq (car g)))
          (y (p+ monr (cdr g))))
        (if (not (or (zero? (vector-ref x 0)) (zero? (vector-ref y 0))))
            (if (tip (p* x y)) (list g tip p* f) (r)) (r))))))))