有限域(4)——程序实现有限域的运算

　　版权申明：本文为博主窗户(Colin Cai)原创，欢迎转帖。如要转贴，必须注明原文网址

　　http://www.cnblogs.com/Colin-Cai/p/9652019.html 

　　作者：窗户

　　QQ/微信：6679072

　　E-mail：6679072@qq.com

　　前一章，我们知道了使用素域的多项式环的商环构造任意的有限域的方法。这一章里，我们就用程序实现任意有限域里的运算。

　　我在这里还是同第一章一样，选择用Scheme来描述。

　　

　　数据表示

　　首先，我们需要的一个参数就是域的特征，记为 p

　　根据上章分析，我们还需要一个不可分多项式，称为生成多项式，记为 poly

　　上一章还有一个重要结论是，如果poly的次数是n, 那么所有的次数小于n的多项式共计pⁿ个，分别属于不同的商集，除此之外并无其他商集，那么也就对应着最后的域的元。

　　所以，我们这里就以这pⁿ个小于n次的多项式来代表该域的所有元，这种表示是合理的。

　　

　　现在，我们来考虑不可分多项式以及域里每个元的计法:

　　我们用一个长度为n+1的list来记录生成多项式poly就可以了(n次多项式一共n+1个系数)，这个list的第0项、第1项、...第n-1项、第n项分别对应着poly的n次系数、n-1次系数、...1次系数、0次系数(常数项)。比如，特征2素域下的3次不可分多项式 x³+x+1，记作(1  0 1 1)。

　　而对于每个域里的元，我们也用上述的写法来记，只是，我们把这个多项式补齐成n-1次，没有的项系数都为0，比如在2特征、生成多项式x³+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求逆来说，因为pⁿ阶域里，除0之外的所有的元乘法构成一个群，所以任何一个元的pⁿ-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次多项式只有一个，x²+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)))
                        "
"
                        (format " ~a/~a=~a
"
                         (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)