版權申明：本文為博主窗戶(Colin Cai)原創，歡迎轉帖。如要轉貼，必須注明原文網址 http://www.cnblogs.com/Colin-Cai/p/9652019.html 作者：窗戶 QQ/微信：6679072 E-mail：6679072@qq.com
前一章，我們知道了使用素域的多項式環的商環構造任意的有限域的方法。這一章里，我們就用程序實現任意有限域里的運算。
我在這里還是同第一章一樣，選擇用Scheme來描述。
數據表示
首先，我們需要的一個參數就是域的特征，記為 p
根據上章分析，我們還需要一個不可分多項式，稱為生成多項式，記為 poly
上一章還有一個重要結論是，如果poly的次數是n, 那么所有的次數小於n的多項式共計pn個，分別屬於不同的商集，除此之外並無其他商集，那么也就對應着最后的域的元。
所以，我們這里就以這pn個小於n次的多項式來代表該域的所有元，這種表示是合理的。
現在，我們來考慮不可分多項式以及域里每個元的計法:
我們用一個長度為n+1的list來記錄生成多項式poly就可以了(n次多項式一共n+1個系數)，這個list的第0項、第1項、...第n-1項、第n項分別對應着poly的n次系數、n-1次系數、...1次系數、0次系數(常數項)。比如，特征2素域下的3次不可分多項式 x3+x+1，記作(1 0 1 1)。
而對於每個域里的元，我們也用上述的寫法來記，只是，我們把這個多項式補齊成n-1次，沒有的項系數都為0，比如在2特征、生成多項式x3+x+1的域下，元[x+1]記為(0 1 1)，也就是一個域下所有的元，記錄成list都是生成多項式的次數這么長。
函數定義
我們再來定義函數接口，我們就這樣去表示有限域里的加減乘除：
(add poly p a b) 代表p特征、生成多項式為poly時，a+b
(sub poly p a b) 代表p特征、生成多項式為poly時，a-b
(mul poly p a b) 代表p特征、生成多項式為poly時，a*b
(div poly p a b) 代表p特征、生成多項式為poly時，a/b
當然，除法涉及到求逆，
(inv poly p a) 代表p特征、生成多項式為poly時，a的逆
程序設計
有了上述准備，再來看第一章的素域里的加減乘除：
;加法 (define (addp p a b) (if (>= (+ a b) p) (- (+ a b) p) (+ a b) ) ) ;減法 (define (subp p a b) (addp p a (if (zero? b) 0 (- p b))) ) ;乘法 (define (mulp p a b) (remainder (* a b) p) ) ;求逆 (define (invp p a) (define (pow n a b) (cond ((zero? n) a) ((zero? (remainder n 2)) (pow (quotient n 2) a (mulp p b b))) (else (pow (quotient n 2) (mulp p a b) (mulp p b b))) ) ) (pow (- p 2) 1 a) ) ;除法 (define (divp p a b) (mulp p a (invp p b)) )
加法和減法很容易做，一個map，對應的都是相同次數的項的系數，然后做合並同類項即可。
(define (addp p a b) (if (>= (+ a b) p) (- (+ a b) p) (+ a b) ) ) (define (subp p a b) (addp p a (if (zero? b) 0 (- p b))) )
而對於inv求逆來說，因為pn階域里，除0之外的所有的元乘法構成一個群，所以任何一個元的pn-1 次冪就是自己的逆元。
(define (inv poly p a) (define (pow poly p a n) (if (zero? n) (append (make-list (- (length a) 1) 0) '(1)) ;指數為0的時候就是1元 (mul poly p a (pow poly p a (- n 1))) ;否則就是a和a的n-1次冪的乘積 ) ) (define (pow-n p n) (if (zero? n) 1 ;指數為0的時候就是1 (* p (pow-n p (- n 1))) ;否則就是p和p的n-1次冪的乘積 ) ) (pow poly p a (- (pow-n p (length a)) 2 )) )
於是，除法也就順理成章：
(define (div poly p a b)
(mul poly p a (inv poly p b))
)
當然，上面的求冪算法是個效率很差的算法。我們可以改個對數級的算法，此處不多解釋。理論依據我之前的文章《RSA簡介(二)——模冪算法》
於是我們效率高的求逆如下：
(define (inv poly p a) (define (pow poly p a n) (define (pow2 n r x) (cond ((zero? n) r) ((zero? (remainder n 2)) (pow2 (quotient n 2) r (mul poly p x x))) (else (pow2 (quotient n 2) (mul poly p r x) (mul poly p x x))) ) ) (pow2 n (append (make-list (- (length a) 1) 0) '(1)) a) ) (define (pow-n p n) (define (pow2 n r x) (cond ((zero? n) r) ((zero? (remainder n 2)) (pow2 (quotient n 2) r (* x x))) (else (pow2 (quotient n 2) (* r x) (* x x))) ) ) (pow2 n 1 p) ) (pow poly p a (- (pow-n p (length a)) 2 )) )
關鍵是乘法，做起來有些難，對lisp沒有太多基礎且不習慣函數式編程的人可能很難寫出來。
兩個元的乘法結果應該是其多項式乘法然后對生成多項式做帶余除法得到的余數。
從而，可以先構造兩個函數，一個函數是多項式乘法：
(define (poly-mul p a b) (map (lambda (x) (remainder x p)) (cdr (foldr (lambda (x y) (cons (+ 1 (car y)) (map + (append (make-list (- (length b) (car y) 1) 0) (map (lambda (c) (* c x)) a) (make-list (car y) 0)) (cdr y) ) ) ) (make-list (+ (length a) (length b)) 0) b ) ) ) )
一個用迭代實現的帶余除法，如下
(define (division a b)
(define (division-iter a b b_len times)
(cond
((zero? times) a)
(else
(division-iter
(append
(map
(lambda (x y) (modulo (- x (* y (divp p (car a) (car b)))) p))
(cdr (take a b_len))
(cdr b)
)
(drop a b_len)
)
b
b_len
(- times 1)
)
)
)
)
(division-iter a b (length b) (- (length a) (length b) -1))
)
有了上述基礎，就可以構造乘法，乘法就是
(define (mul poly p a b) (division (poly-mul p a b) poly ) )
多項式乘法poly-mul的實現中用到了foldr算子，這個算子並不是定義在Scheme任何標准中的，獨屬於racket。不過這個算子實際上就是著名的SICP里的accumulate算子，我也一直很奇怪如此重要的算子為什么不在標准中體現。以下是這個算子的定義。
(define (foldr op init lst) (if (null? lst) init (op (car lst) (foldr op init (cdr lst))) ) )
測試
以上，域的加減乘除都有了。
我們可以測試一下，階最小的不是素域的有限域是4階，特征為2。
特征2素域下的不可分2次多項式只有一個，x2+x+1
(define poly '(1 1 1)) ;生成多項式 (define p 2) ;特征 (define (list-all) (define (sub1 x) (foldr (lambda (m n) (if (zero? (apply * n)) (cons m n) (cons (- 1 m) n))) '() x) ) (define (func2 x y) (let* ((a (car x)) (b (cdr x)) (asum (apply + a)) (bsum (apply + b)) ) (cond ((zero? (+ asum bsum)) (cons x y)) ((zero? bsum) (func2 (cons (sub1 a) (sub1 b)) (cons x y))) (else (func2 (cons a (sub1 b)) (cons x y))) ) ) ) (let ((m (make-list (- (length poly) 1) (- p 1)))) (func2 (cons m m) '()) ) ) (for-each (lambda (x) (display (string-append (format "~a+~a=~a ~a-~a=~a ~a*~a=~a" (car x) (cdr x) (add p (car x) (cdr x)) (car x) (cdr x) (sub p (car x) (cdr x)) (car x) (cdr x) (mul poly p (car x) (cdr x)) ) (if (zero? ( foldr + 0 (cdr x))) "\n" (format " ~a/~a=~a\n" (car x) (cdr x) (div poly p (car x) (cdr x)) ) ) ) ) ) (list-all) )
上面測試代碼中poly和p可以任意改，這里只以4階域為測試對象。
以上運行結果如下：
(0 0)+(0 0)=(0 0) (0 0)-(0 0)=(0 0) (0 0)*(0 0)=(0 0)
(0 0)+(0 1)=(0 1) (0 0)-(0 1)=(0 1) (0 0)*(0 1)=(0 0) (0 0)/(0 1)=(0 0)
(0 0)+(1 0)=(1 0) (0 0)-(1 0)=(1 0) (0 0)*(1 0)=(0 0) (0 0)/(1 0)=(0 0)
(0 0)+(1 1)=(1 1) (0 0)-(1 1)=(1 1) (0 0)*(1 1)=(0 0) (0 0)/(1 1)=(0 0)
(0 1)+(0 0)=(0 1) (0 1)-(0 0)=(0 1) (0 1)*(0 0)=(0 0)
(0 1)+(0 1)=(0 0) (0 1)-(0 1)=(0 0) (0 1)*(0 1)=(0 1) (0 1)/(0 1)=(0 1)
(0 1)+(1 0)=(1 1) (0 1)-(1 0)=(1 1) (0 1)*(1 0)=(1 0) (0 1)/(1 0)=(1 1)
(0 1)+(1 1)=(1 0) (0 1)-(1 1)=(1 0) (0 1)*(1 1)=(1 1) (0 1)/(1 1)=(1 0)
(1 0)+(0 0)=(1 0) (1 0)-(0 0)=(1 0) (1 0)*(0 0)=(0 0)
(1 0)+(0 1)=(1 1) (1 0)-(0 1)=(1 1) (1 0)*(0 1)=(1 0) (1 0)/(0 1)=(1 0)
(1 0)+(1 0)=(0 0) (1 0)-(1 0)=(0 0) (1 0)*(1 0)=(1 1) (1 0)/(1 0)=(0 1)
(1 0)+(1 1)=(0 1) (1 0)-(1 1)=(0 1) (1 0)*(1 1)=(0 1) (1 0)/(1 1)=(1 1)
(1 1)+(0 0)=(1 1) (1 1)-(0 0)=(1 1) (1 1)*(0 0)=(0 0)
(1 1)+(0 1)=(1 0) (1 1)-(0 1)=(1 0) (1 1)*(0 1)=(1 1) (1 1)/(0 1)=(1 1)
(1 1)+(1 0)=(0 1) (1 1)-(1 0)=(0 1) (1 1)*(1 0)=(0 1) (1 1)/(1 0)=(1 0)
(1 1)+(1 1)=(0 0) (1 1)-(1 1)=(0 0) (1 1)*(1 1)=(1 0) (1 1)/(1 1)=(0 1)