SICP 読み (87) 2.5.3 例: 記号代数

あと 4 頁がなかなか進まん。来週以降は対応可能な時間も減る。

問題 2.91

問題 2.88 の解をコピッて作成する方向で。まず試験か。div-poly は商と剰余なリストを返す、という事に。

  ("2-91"
   (setup (lambda ()
	    (install-scheme-number-package)
	    (install-polynomial-package)))

   ("example"
    (assert-equal '(polynomial x (((3 1) (1 1)) ((1 1) (0 -1))))
		  (div (make-polynomial 'x '((5 1) (0 -1)))
		       (make-polynomial 'x '((2 1) (0 -1)))))
    )

   ("x^2 - 1/ x^2 - 1"
    (assert-equal '(polynomial x (((0 1)) ()))
		  (div (make-polynomial 'x '((2 1) (0 -1)))
		       (make-polynomial 'x '((2 1) (0 -1)))))
    )

   ("x^2 + 2x + 1 / x + 1"
    (assert-equal '(polynomial x (((1 1) (0 1)) ()))
		  (div (make-polynomial 'x '((2 1) (1 2) (0 1)))
		       (make-polynomial 'x '((1 1) (0 1)))))
    )
   )

なんか微妙。結果としては

'((polynomial x ((1 1) (0 1))) (polynomial x ()))

の方が多項式としては自然に見える。この方向で実装を検討してみる。

検討

処理としては

  1. 被除数の最高次の項を除数の最高次の項で割る。この結果が商の第一項。(make-term new-o new-c) すれば良い。
  2. この結果に除数を掛ける。
  3. 上の結果を被除数から引く。これと除数を div-term に再帰的に渡す。
  4. 除数の次数が被除数の次数を越えたら停止し、その被除数を剰余とする。

最初とケツのソレは提示されてますな。2 番目と 3 番目が_結果の残りを再帰的に (以下略)_な部分ですか。
で、rest-of-result には商と剰余のリストが入ってるハズなんで、完全な結果は

(list (cons (make-term new-o newc) (car rest-of-result)) (cadr rest-of-result))

で形成できるんかなぁ。

実装

上記を踏まえ、以下の実装を。require してる 2.5.1 は略。

(require "2.5.1")

(define (install-polynomial-package)
  (define (adjoin-term term term-list)
    (if (=zero? (coeff term))
	term-list
	(cons term term-list)))
  (define (the-empty-termlist) '())
  (define (first-term term-list) (car term-list))
  (define (rest-terms term-list) (cdr term-list))
  (define (empty-termlist? term-list) (null? term-list))
  (define (make-term order coeff) (list order coeff))
  (define (order term) (car term))
  (define (coeff term) (cadr term))

  (define (make-poly variable term-list)
    (cons variable term-list))
  (define (variable p) (car p))
  (define (term-list p) (cdr p))

  (define (variable? x) (symbol? x))
  (define (same-variable? v1 v2)
    (and (variable? v1) (variable? v2) (eq? v1 v2)))

  (define (add-poly p1 p2)
    (if (same-variable? (variable p1) (variable p2))
	(make-poly (variable p1)
		   (add-terms (term-list p1)
			      (term-list p2)))
	(error "Polys not in same var --- ADD-POLY"
	       (list p1 p2))))

  (define (add-terms L1 L2)
    (cond ((empty-termlist? L1) L2)
	  ((empty-termlist? L2) L1)
	  (else
	   (let ((t1 (first-term L1)) (t2 (first-term L2)))
	     (cond ((> (order t1) (order t2))
		    (adjoin-term
		     t1 (add-terms (rest-terms L1) L2)))
		   ((< (order t1) (order t2))
		    (adjoin-term
		     t2 (add-terms L1 (rest-terms L2))))
		   (else
		    (adjoin-term
		     (make-term (order t1)
				(add (coeff t1) (coeff t2)))
		     (add-terms (rest-terms L1)
				(rest-terms L2)))))))))

  (define (mul-poly p1 p2)
    (if (same-variable? (variable p1) (variable p2))
	(make-poly (variable p1)
		   (mul-terms (term-list p1)
			      (term-list p2)))
	(error "Polys not in same var --- MUL-POLY"
	       (list p1 p2))))

  (define (mul-terms L1 L2)
    (if (empty-termlist? L1)
	(the-empty-termlist)
	(add-terms (mul-term-by-all-terms (first-term L1) L2)
		   (mul-terms (rest-terms L1) L2))))

  (define (mul-term-by-all-terms t1 L)
    (if (empty-termlist? L)
	(the-empty-termlist)
	(let ((t2 (first-term L)))
	  (adjoin-term
	   (make-term (+ (order t1) (order t2))
		      (mul (coeff t1) (coeff t2)))
	   (mul-term-by-all-terms t1 (rest-terms L))))))

  (define (neg-poly p)
    (make-poly (variable p)
	       (neg-terms (term-list p))))

  (define (neg-terms L)
    (if (empty-termlist? L)
	(the-empty-termlist)
	(let ((t (first-term L)))
	  (adjoin-term
	   (make-term (order t) (- (coeff t)))
	   (neg-terms (rest-terms L)))))
     )

  (define (div-poly p1 p2)
    (define (quotient solution)
      (car solution)
      )
    (define (reminder solution)
      (cadr solution)
      )
    (if (same-variable? (variable p1) (variable p2))
	(let ((solution (div-terms (term-list p1) (term-list p2))))
	  (list (make-poly (variable p1)
			   (quotient solution))
		(make-poly (variable p1)
			   (reminder solution))))
	(error "Polys not in same var --- DIV-POLY"
	       (list p1 p2)))
    )

  (define (div-terms L1 L2)
    (if (empty-termlist? L1)
	(list (the-empty-termlist) (the-empty-termlist))
	(let ((t1 (first-term L1))
	      (t2 (first-term L2)))
	  (if (> (order t2) (order t1))
	      (list (the-empty-termlist) L1)
	      (let ((new-c (div (coeff t1) (coeff t2)))
		    (new-o (- (order t1) (order t2))))
		(let ((rest-of-result (div-terms 
				       (add-terms 
					(neg-terms 
					 (mul-term-by-all-terms 
					  (make-term new-o new-c)
					  L2))
					L1)
				       L2)))
		  (list (cons (make-term new-o new-c) (car rest-of-result))
			(cadr rest-of-result))
		  ))))))

  (define (tag p) (attach-tag 'polynomial p))
  (put 'div '(polynomial polynomial)
       (lambda (p1 p2)
	 (let ((result (div-poly p1 p2)))
	   (list (tag (car result)) (tag (cadr result))))))
  (put '=zero? '(polynomial)
       (lambda (p) (empty-termlist? (term-list p))))
  (put 'neg '(polynomial)
       (lambda (p) (tag (neg-poly p))))
  (put 'add '(polynomial polynomial)
       (lambda (p1 p2) (tag (add-poly p1 p2))))
  (put 'mul '(polynomial polynomial)
       (lambda (p1 p2) (tag (mul-poly p1 p2))))
  (put 'make 'polynomial
       (lambda (var terms) (tag (make-poly var terms))))
  'done)

(define (make-polynomial var terms)
  ((get 'make 'polynomial) var terms))
(define (neg x) (apply-generic 'neg x))

で試験してみたんですが、多項式のソレをまた間違えてた。正しくは

'((polynomial x ((1 1) (0 1))) (polynomial x ()))

ぢゃなくて

'((polynomial x (1 1) (0 1)) (polynomial x ()))

だな。

あ、違う。こうだ。

'((polynomial x (1 1) (0 1)) (polynomial x))

一応試験は通ったんですが、put してなくて、な不具合あり。試験は以下。

  ("2-91"
   (setup (lambda ()
	    (install-scheme-number-package)
	    (install-polynomial-package)))

   ("example"
    (assert-equal '((polynomial x (3 1) (1 1))
		    (polynomial x (1 1) (0 -1)))
		  (div (make-polynomial 'x '((5 1) (0 -1)))
		       (make-polynomial 'x '((2 1) (0 -1)))))
    )

   ("x^2 - 1/ x^2 - 1"
    (assert-equal '((polynomial x (0 1)) 
		    (polynomial x))
		  (div (make-polynomial 'x '((2 1) (0 -1)))
		       (make-polynomial 'x '((2 1) (0 -1)))))
    )

   ("x^2 + 2x + 1 / x + 1"
    (assert-equal '((polynomial x (1 1) (0 1))
		    (polynomial x))
		  (div (make-polynomial 'x '((2 1) (1 2) (0 1)))
		       (make-polynomial 'x '((1 1) (0 1)))))
    )
   )