private/galois.rkt
#lang racket/base

(require racket/vector)

(provide (all-defined-out))

;;
;; Addition (subtraction) of polynomials with coefficients in GF(2)
;; is simply XOR
;;
(define (gf+ a b)
  (bitwise-xor a b))

;;
;; Multiply two polynomials using the peasant's algoithm, reducing
;; modulo (x^8 + x^4 + x^3 + x^2 + 1)
;;
(define (gf* a b)
  (let ([m #x11D])
    (do ([x 0 (if (bitwise-bit-set? b 0)
                  (bitwise-xor x a)
                  x)]
         [a a (if (bitwise-bit-set? a 7)
                  (bitwise-xor (arithmetic-shift a 1) m)
                  (arithmetic-shift a 1))]
         [b b (arithmetic-shift b -1)])
      ((zero? b) (bitwise-and x #xFF)))))

;;
;; Compute the residue of a modulo m
;;
(define (gf-mod a m)
  (let ([a-len (integer-length a)]
        [m-len (integer-length m)])
    (do ([i (- a-len m-len) (sub1 i)]
         [a a (if (bitwise-bit-set? a (sub1 (+ m-len i)))
                  (bitwise-xor a (arithmetic-shift m i))
                  a)])
      ((< i 0) a))))

;;
;; Compute coefficients of polynomial g(z) = (z - α) ... (z - α^r)
;; which are polynomials in 2[X] modulo the polynomial above.
;;
(define (make-rs-generator r)
  (do ([i 0 (add1 i)]
       [a 1 (gf* a 2)]
       [gs (vector-append #(1) (make-vector (sub1 r) 0)) ; Start with g(x) = 1
          ;; Multiply g(x) by (x - 2^i) => x g(x) + 2^i g(x)
          (vector-map gf+
            (vector-append #(0) (vector-drop-right gs 1))
            (vector-map (lambda (g)
                          (gf* g a))
                        gs))])
    ((= i r) gs)))