Gauche でグレブナ基底計算

D加群セミナーの関係で、
グレブナ基底と代数多様体入門〈上〉イデアル・多様体・アルゴリズム
を読んだ副産物です。
ブッフベルガーのアルゴリズムで一番簡単なやつを実装しただけなので、おもちゃの域をでるものではありません。90式はブリキ缶だぜ!

(use srfi-1)
(use srfi-43)
(use util.combinations)

(define (poly . l)
  (apply hash-table (cons 'equal? (remove (lambda (x) (zero? (cdr x))) l))))

(define (order-lex a b)
  (let loop ((a a)
	     (b b))
    (cond 
     ((null? a) #f)
     ((null? b) #t)
     ((> (car a) (car b)) #t)
     ((< (car a) (car b)) #f)
     (else (loop (cdr a) (cdr b))))))

(define (order-grlex a b)
  (let ((sum-a (fold + 0 a))
	(sum-b (fold + 0 b)))
    (cond 
     ((> sum-a sum-b) #t)
     ((< sum-a sum-b) #f)
     ((= sum-a sum-b) (order-lex a b)))))

(define (order-grevlex a b)
  (let ((sum-a (fold + 0 a))
	(sum-b (fold + 0 b)))
    (cond 
     ((> sum-a sum-b) #t)
     ((< sum-a sum-b) #f)
     ((= sum-a sum-b) (order-lex (reverse b) (reverse a))))))
    
(define (multideg order f)
  (fold (lambda (a b) (if (order a b) a b)) '() (hash-table-keys f)))

(define (LC order f)
  (hash-table-get f (multideg order f)))

(define (LM order f)
  (cons (multideg order f) 1))

(define (LT order f)
  (cons (multideg order f) (LC order f)))

(define (divide-term a b)
  (cons (map - (car a) (car b)) (/ (cdr a) (cdr b))))  

(define (zero-poly? f)
  (every zero? (hash-table-values f)))

(define (+poly a b)
  (define ht (poly))
  (hash-table-for-each a (lambda (akey avalue)
			   (hash-table-put! ht akey avalue)))
  (hash-table-for-each b (lambda (bkey bvalue)
			   (hash-table-update! ht bkey (cut + <> bvalue) 0)))
  (hash-table-for-each ht (lambda (key value)
			    (when (zero? value) (hash-table-delete! ht key))))
  ht)

(define (-poly a b)
  (define ht (poly))
  (hash-table-for-each a (lambda (akey avalue)
			   (hash-table-put! ht akey avalue)))
  (hash-table-for-each b (lambda (bkey bvalue)
			   (hash-table-update! ht bkey (cut - <> bvalue) 0)))
  (hash-table-for-each ht (lambda (key value)
			    (when (zero? value) (hash-table-delete! ht key))))
  ht)

(define (*poly a b)
  (define ht (poly))
  (hash-table-for-each a (lambda (akey avalue)
			   (hash-table-for-each  b (lambda (bkey bvalue)
						     (hash-table-update! ht (map + akey bkey) (cut + <> (* avalue bvalue)) 0)))))
  (hash-table-for-each ht (lambda (key value)
			    (when (zero? value) (hash-table-delete! ht key))))
  ht)

(define (quotient&remainder order f g-list)
  (define quotients (make-vector (length g-list) (poly)))
  (let loop ((f f)
	     (intermediate-divident (poly)))
    (if (zero-poly? f)
	(values (vector->list quotients) intermediate-divident)
	(let ((ind (list-index (lambda (x)
				 (every >= (multideg order f) (multideg order x)))
			       g-list)))
	  (if ind
	      (let* ((g (list-ref g-list ind))
		     (p (poly (divide-term (LT order f) (LT order g))))) 
		(vector-set! quotients ind (+poly p (vector-ref quotients ind)))
		(loop (-poly f (*poly p g)) intermediate-divident))
	      (loop (-poly f (poly (LT order f))) (+poly (poly (LT order f)) intermediate-divident)))))))

(define (quotient-poly order f g-list)
  (receive (q r) (quotient&remainder order f g-list)
    q))

(define (remainder-poly order f g-list)
  (receive (q r) (quotient&remainder order f g-list)
    r))

(define (S-poly order f g)
  (define term-r (cons (map max (multideg order f) (multideg order g)) 1))
  (-poly (*poly (poly (divide-term term-r (LT order f))) f)
 	 (*poly (poly (divide-term term-r (LT order g))) g)))

(define (groebner order f-list)
  (let ((l (remove zero-poly? (map (lambda (x) (remainder-poly order (S-poly order (car x) (cadr x)) f-list))
				   (combinations f-list 2)))))
    (if (null? l)
	f-list
	(groebner order (append f-list l)))))

データ構造は、単項式 $xy^{2}z$ を'(1 2 1)、項 $4xy^{2}z$ を'((1 2 1) . 4)、多項式をハッシュテーブルで(キーに単項式、値にその係数)としました。
多項式 $f = 4xy^{2}z + 4z^{2} - 5x^{3} + 7x^{2}z^{2}$ について、辞書式順序の場合、以下が成り立つ。
$multideg(f) = (3,0,0)$
$LC(f) = -5$
$LM(f) = x^{3}$
$LT(f) = -5x^{3}$

(define f (poly '((1 2 1) . 4) '((0 0 2) . 4) '((3 0 0) . -5) '((2 0 2) . 7)))
(multideg order-lex f) == (3 0 0)
(LC order-lex f) == -5
(LM order-lex f) == ((3 0 0) . 1)
(LT order-lex f) == ((3 0 0) . -5)

S多項式はこうです。
$f = x^{3}y^{2} - x^{2}y^{3} + x, g = 3x^{4}y + y^{2}$ のときgrlex順序に関するS多項式は
$S(f,g) = -x^{3}y^{3} + x^{2} - (\frac{1}{3})y^{3}$

(S-poly order-grlex (poly '((3 2) . 1)
                          '((2 3) . -1)
                          '((1 0) . 1))
                    (poly '((4 1) . 3)
                          '((0 2) . 1)))
                          
== '(((3 3) . -1) ((2 0) . 1) ((0 3) . -1/3))

グレブナ基底はこうです。
$I = \left\langle f_{1},f_{2} \right\rangle = \left\langle x^{3}-2xy, x^{2}y - 2y^{2} + x \right\rangle$ のgrlex順序に関するグレブナ基底は、
$F=\left\{ f_{1},f_{2},f_{3},f_{4},f_{5} \right\} = \left\{ x^{3}-2xy, x^{2}y - 2y^{2} + x, -x^{2}, -2xy, -2y^{2} + x \right\}$

(groebner order-grlex
          (list (poly '((3 0) . 1)
                      '((1 1) . -2))
                (poly '((2 1) . 1)
                      '((0 2) . -2)
                      '((1 0) . 1))))
    
== '((((3 0) . 1) ((1 1) . -2))
     (((2 1) . 1) ((0 2) . -2) ((1 0) . 1))
     (((2 0) . -1))
     (((1 1) . -2))
     (((0 2) . -2) ((1 0) . 1)))