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

pseudoremainder-terms を盛り込んでみた。試験 NG。
こんな感じ。

$ test/run-test.scm -vv
- (test suite) 2.5.3
-- (test case) 2-95: .F
 expected:<(polynomial x (2 18954) (1 -37908) (0 18954))>
 but was:<(polynomial x (2 1458) (1 -2916) (0 1458))> in pseudoremainder-terms

-- (test case) 2-94: ....
-- (test case) include rational package: .....
-- (test case) 2-88: ....
-- (test case) 2-87: ...
-- (test case) make-polynomial test: ...
-- (test case) add test (1): ........
-- (test case) mul test: .........

38 tests, 80 assertions, 79 successes, 1 failures, 0 errors
Testing time: 0.033306
$

うーん。どーゆー事だ、と言いつつ確認。

gosh> (/ 18954 1458)
13
gosh> (/ 37908 2916)
13
gosh>

13 ってどっかで見た数字だよなぁ。何か勘違いしてるのかしらん。

問題 2.95 を再度

で、もう一度置き換えにトライしたりなんかしてたんですが、結果は微妙。実装をよくよく見てみると

                         (let ((rest-of-result (div-terms substract L2)))

ってなっている。ちなみに手続き名は div-terms-p で前のバージョンの div-terms も残っている状態。要は div-terms-p を div-terms のコピーで作ったのは良いが、再帰呼び出しの所 (上記) は修正を忘れていた、という訳で。

もしかして ... と思い div-terms-p を呼び出す形に修正したら試験パス。な、なんでよ、と言いつつテキストを見ると昨晩のエントリのようにいちいち整数化因子を掛けなくてもできるように読める。これは別途置き換えにトライする必要がありそうやっさ。

でも実装が微妙だなぁ。それは良いとして、とりあえず置き換えてみる

gcd-terms-p から

(gcd-terms-p '((4 11) (3 -22) (2 18) (1 -14) (0 7))
	     '((3 13) (2 -21) (1 3) (0 5)))

(gcd-terms-p '((3 13) (2 -21) (1 3) (0 5))
	     (reminder-terms-p '((4 11) (3 -22) (2 18) (1 -14) (0 7))
			       '((3 13) (2 -21) (1 3) (0 5))))

(gcd-terms-p '((3 13) (2 -21) (1 3) (0 5))
	     (cadr (div-terms-p '((4 11) (3 -22) (2 18) (1 -14) (0 7))
				'((3 13) (2 -21) (1 3) (0 5)))))

div-terms-p

で、div-terms-p を

(div-terms-p '((4 11) (3 -22) (2 18) (1 -14) (0 7))
	     '((3 13) (2 -21) (1 3) (0 5)))

;; int-factor の準備
(make-term 0 (expt 13 (+ 1 (- 4 3))))
;; -> '(0 169)

;; pL の準備
(mul-term-by-all-terms '(0 169) '((4 11) (3 -22) (2 18) (1 -14) (0 7)))
;; -> '((4 1859) (3 -3718) (2 3042) (1 -2366) (0 1183))

;; t1 は '(4 1859)
;; t2 は '(3 13)
;; (> (order t2) (order t1)) は #f

;; new-o は 1 new-c は 143
(mul-term-by-all-terms '(1 143) '((3 13) (2 -21) (1 3) (0 5)))
;; -> '((4 1859) (3 -3003) (2 429) (0 715))
;; neg すると '((4 -1859) (3 3003) (2 -429) (1 -715))
;; pL と add-terms
(add-terms '((4 -1859) (3 3003) (2 -429) (1 -715))
	   '((4 1859) (3 -3718) (2 3042) (1 -2366) (0 1183)))
;; -> '((3 -715) (2 2613) (1 -3081) (0 1183))

(div-terms '((3 -715) (2 2613) (1 -3081) (0 1183))
	   '((3 13) (2 -21) (1 3) (0 5)))

div-terms 呼び出し

次は div-terms になる

;; t1 は '(3 -715)
;; t2 は '(3 13)
;; (> (order t2) (order t1)) は #f

;; new-o は 0 new-c は -55
(mul-term-by-all-terms '(0 -55) '((3 13) (2 -21) (1 3) (0 5)))
;; -> '((3 -715) (2 1155) (1 -165) (0 -275))
;; neg すると '((3 715) (2 -1155) (1 165) (0 275))
;; L1 と add-terms
(add-terms '((3 715) (2 -1155) (1 165) (0 275))
	   '((3 -715) (2 2613) (1 -3081) (0 1183)))
;; -> '((2 1458) (1 - 2916) (0 1458))

(div-terms '((2 1458) (1 - 2916) (0 1458))
	   '((3 13) (2 -21) (1 3) (0 5)))

div-terms 呼び出し (2)

;; t1 は '(2 1458)
;; t2 は '(3 13)
;; (> (order t2) (order t1)) は真

;; (list '() '((2 1458) (1 -2916) (0 1458))) を返却

戻り

どこに戻るか、というと直前の dir-terms か

;; rest-of-result に '(() ((2 1458) (1 -2916) (0 1458))) が設定
;; new-o は 0 で new-c は -55
(list (cons '(0 -55) '()) '((2 1458) (1 -2916) (0 1458)))
;; -> '(((0 -55)) ((2 1458) (1 -2916) (0 1458))) が戻る

さらに戻り (div-terms-p に戻る)

;; rest-of-result に '(((0 -55)) ((2 1458) (1 -2916) (0 1458))) が設定
;; new-o は 1 で new-c は 143
(list (cons '(1 143) '((0 -55))) '((2 1458) (1 -2916) (0 1458)))
;; -> '(((1 143) (0 -55)) ((2 1458) (1 -2916) (0 1458)))

gcd-terms-p 方面に戻る

(gcd-terms-p '((3 13) (2 -21) (1 3) (0 5))
	     (cadr (div-terms-p '((4 11) (3 -22) (2 18) (1 -14) (0 7))
				'((3 13) (2 -21) (1 3) (0 5)))))

(gcd-terms-p '((3 13) (2 -21) (1 3) (0 5))
	     (cadr '(((1 143) (0 -55)) ((2 1458) (1 -2916) (0 1458)))))

(gcd-terms-p '((3 13) (2 -21) (1 3) (0 5))
	     '((2 1458) (1 -2916) (0 1458)))

(gcd-terms-p '((2 1458) (1 -2916) (0 1458))
	     (reminder-terms-p '((3 13) (2 -21) (1 3) (0 5))
			       '((2 1458) (1 -2916) (0 1458))))

(gcd-terms-p '((2 1458) (1 -2916) (0 1458))
	     (cadr (div-terms-p '((3 13) (2 -21) (1 3) (0 5))
				'((2 1458) (1 -2916) (0 1458)))))

もう一度、div-terms-p へ

答えは見えてるんだけどな。

(div-terms-p '((3 13) (2 -21) (1 3) (0 5))
	     '((2 1458) (1 -2916) (0 1458)))

;; int-factor の準備
(make-term 0 (expt 1458 (+ 1 (- 3 2))))
;; -> '(0 2125764)

;; pL の準備
(mul-term-by-all-terms '(0 2125764) '((3 13) (2 -21) (1 3) (0 5))
;; -> '((3 27634932) (2 -44641044) (1 6377292) (0 10628820))

;; t1 は '(3 27634932)
;; t2 は '(2 1458)
;; (> (order t2) (order t1)) は #f

;; new-o は 1 new-c は 18954
(mul-term-by-all-terms '(1 18954) '((2 1458) (1 -2916) (0 1458)))
;; -> '((3 27634932) (2 -55269864) (1 27634932))
;; neg すると '((3 -27634932) (2 55269864) (1 -27634932))
;; pL と add-terms
(add-terms '((3 -27634932) (2 55269864) (1 -27634932))
	   '((3 27634932) (2 -44641044) (1 6377292) (0 10628820)))
;; -> '((2 10628820) (1 -21257640) (0 10628820))

(div-terms '((2 10628820) (1 -21257640) (0 10628820))
	   '((2 1458) (1 -2916) (0 1458)))

div-terms 呼び出し

次は div-terms

;; t1 は '(2 10628820)
;; t2 は '(2 1458)
;; (> (order t2) (order t1)) は #f

;; new-o は 0 new-c は 7290
(mul-term-by-all-terms '(0 7290) '((2 1458) (1 -2916) (0 1458)))
;; -> '((2 10628820) (1 -21257640) (0 10628820))
;; neg すると '((2 -10628820) (1 21257640) (0 -10628820))
;; L1 と add-terms
(add-terms '((2 -10628820) (1 21257640) (0 -10628820))
	   '((2 10628820) (1 -21257640) (0 10628820))
;; -> '()

(div-terms '()
	   '((2 1458) (1 - 2916) (0 1458)))

div-terms 呼び出し (多分最後)

L1 が '() なので (list '() '()) を返却。

戻り

;; rest-of-result は '(()())
;; new-o が 0 で new-c は 7290
(list (cons '(0 7290) '()) '())
;; -> '(((0 7290)) ())

さらに戻り

;; rest-of-result は '(((0 7290)) ())
;; new-o が 1 で new-c が 18954
(list (cons '(1 18954) '((0 7290))) '())
;; -> '(((1 18954) (0 7290)) ())

gcd-terms-p 方面へ

(gcd-terms-p '((2 1458) (1 -2916) (0 1458))
	     (cadr (div-terms-p '((3 13) (2 -21) (1 3) (0 5))
				'((2 1458) (1 -2916) (0 1458)))))

(gcd-terms-p '((2 1458) (1 -2916) (0 1458))
	     (cadr '(((1 18954) (0 7290)) ())))

(gcd-terms-p '((2 1458) (1 -2916) (0 1458)) '())

で、gcd-terms-p は L2 が '() だったら L1 を戻すので解は

'(polynomial x (2 1458) (1 -2916) (0 1458))

ってなるのか。最初の一発だけで OK だったとは (を

問題 2.96 の実装

たまたま div-terms-p って別名な手続きを作成して、その中で div-term を呼んでいて気がついた、とは言え実装は纏めたい。どうしたものか。

で、微妙な頭をヒネってひり出したのが以下。リファクタリングの余地あり。まず実装。 (一部のみ)

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

  (define (gcd-terms L1 L2)
    (if (empty-termlist? L2)
	L1
	(gcd-terms L2 (pseudoremainder-terms L1 L2))))

  (define (pseudoremainder-terms L1 L2)
    (if (empty-termlist? L1)
	(the-tmpty-termlist)
	(let ((int-factor (make-term 0 
				     (expt (coeff (first-term L2))
					   (+ 1 (- (order (first-term L1))
						   (order (first-term L2))))))))
	  (let ((pL (mul-term-by-all-terms int-factor L1)))
	    (cadr (div-terms pL L2))))))

  (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 ((mul-result-by-divisor (mul-term-by-all-terms
					      (make-term new-o new-c)
					      L2)))
		  (let ((negate (neg-terms mul-result-by-divisor)))
		    (let ((substract (add-terms negate L1)))
		      (let ((rest-of-result (div-terms substract L2)))
			(list (cons (make-term new-o new-c)
				    (car rest-of-result))
			      (cadr rest-of-result)))))))))))

試験が以下。

  ("2-95"
   (setup (lambda ()
	    (install-rational-package)
	    (install-scheme-number-package)
	    (install-polynomial-package)))
   ("not equall"
    (let ((p1 (make-polynomial 'x '((2 1) (1 -2) (0 1))))
	  (p2 (make-polynomial 'x '((2 11) (0 7))))
	  (p3 (make-polynomial 'x '((1 13) (0 5)))))
      (let ((q1 (mul p1 p2))
	    (q2 (mul p1 p3)))
	(assert-equal q1
		      '(polynomial x (4 11) (3 -22) (2 18) (1 -14) (0 7)))
	(assert-equal q2
		      '(polynomial x (3 13) (2 -21) (1 3) (0 5)))
	(assert-not-equal p1 (greatest-common-divisor q1 q2))))
    )

   ("pseudoremainder-terms"
    (let ((p1 (make-polynomial 'x '((2 1) (1 -2) (0 1))))
	  (p2 (make-polynomial 'x '((2 11) (0 7))))
	  (p3 (make-polynomial 'x '((1 13) (0 5)))))
      (let ((q1 (mul p1 p2))
	    (q2 (mul p1 p3)))
	(assert-equal '(polynomial x (2 1458) (1 -2916) (0 1458))
		      (greatest-common-divisor q1 q2))))
    )
   )

もう少しで '(polynomial x (2 1) (1 2) (0 1)) になるんだけどな。問題 2.96 の b については scheme-number な gcd でイケるはず。

全係数を最大公約数で割る

gcd-poly ん中の gcd-terms 呼び出しを_最大公約数で割る手続き_で wrap すれば良い? がしかし、テキストには gcd-terms を modify せぇって書いてあるな。どうやりゃ良いのか。