; http://cap-lore.com/MathPhys/Field/finite
(let ((Do ((fileVal "Do") 'Do))(egcd (fileVal "egcd")))
(lambda (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) (let ((y (rc x))) (vector-set! a x y) (vector-set! a y 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 (> (vector-ref tp 0) 0)(tip tp)) tp (r)))))))))
     (gsip (lambda (m) (let* ((pq (expt p m))(m20 (* m 20))
        (tp (lambda (k) (let ((x (i->p (+ pq k))))(and (tip x) x))))) (or
            (let seek ((n 2)) (and (< n m20) (or (tp n)) (seek (+ n 1))))
            (let seek ((n 1)) (or (tp (+ p n)) (seek (+ n 1))))))))
     (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))))
         (let ((pl (list (cons 'f* f*) (cons 'f/ f/) (cons 'fexpt fexpt))))
         (lambda (sy) (cdr (assq sy pl))))))))
(let ((veq (lambda (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))))))
(let ((pl (map cons 
     '(pqr p->i i->p gap gip gfip gsip veq m+ m- m* m/ p+ p- p* mexpt trim pgcmd tip pegcd fops p)
     (list pqr p->i i->p gap gip gfip gsip veq m+ m- m* m/ p+ p- p* mexpt trim pgcmd tip pegcd fops p))))
  (lambda (sy) (cdr (assq sy pl))))))))))
