Compute the Groebner Basis in Scheme.

I attended class of D-module and read "Ideals, Varieties, and Algorithms".
I wrote a Scheme program to compute the Groebner basis.
Used Gauche, a scheme implementation.

(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)))))

And now, let's test them.

Test.

(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)

Test.

(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))   

Test.

(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)))