SICP 2.5.3 Example: Symbolic Algebra - プログラミング再入門

第２章も漸く大詰め。

ノート

再び登場、代数式による例。Symbolic Algebraはシンボル代数？記号代数？

Arithmetic on polynomials

多項式の計算。
ここでは変数がひとつの一元多項式に限定する。多項式は項の和で構成されて、項は係数（定数）、変数の冪乗、あるいは係数と変数の冪乗の積からなる。係数もまた代数式だが、式全体の変数は含まないが、別の変数は含むかも知れない。
$5x^2+3x+7$ と $5y^2+3y+7$ はひとつの文脈に現れて、xとyが別の変数である事を示している場合でなければ数学的には変数の名前は任意なのでどちらも同じ式を表している。けど、ここでは変数名が違う場合は違う式として扱う。つまり二元の多項式は扱わないので、演算に使う多項式の変数は同じである事を前提とする。また多項式同士の演算は加算と乗算のみを扱う。

（驚いた事に）多項式もひとつの型として前節のシステムに追加する。定数同士の演算に利用する。
型polyは変数（名）と項のリストからなる。項のリストがどう表されるかは後から出て来る。
加算は同じ次数の各項を足す。総称関数としてのaddを使う所がミソ。ここが多項式でも定数でも良いと言う事か。混ざった時に演算出来るのか？との疑問が生じるがやはり「適切なcoercion」を実装する模様。
乗算は片方の式の各項にもう片方の全ての項を乗算して全部足し合わせる。ここでも総称関数のmulを使う。片方の式のひとつの項にもう片方の式の全ての項を掛けた結果の多項式を順に加算して行く。各項の乗算は次数は加算して係数同士はmulで乗算する。

Representing term lists

保留されていた項リストの実装。
係数が０の項が沢山ある様な式の事を考慮して（次数係数）と言うペアのリストにする。（実はこうしないと次数が高い項から始まるリストは作れない）

実際に動かしてみる。 $x^5+2x^4+3x^2-2x-5$ 、 $(x^2+2x+3)+(3x^2+4x+5)=4x^2+6x+8$ と $(x^2+2x+3)\times(3x^2+4x+5)=3x^4+10x^3+22x^2+22x+15$

> (make-polynomial 'x '((5 1) (4 2) (2 3) (1 -2) (0 -5)))
{polynomial x {5 1} {4 2} {2 3} {1 -2} {0 -5}}
> (add (make-polynomial 'x '((2 1) (1 2) (0 3))) (make-polynomial 'x '((2 3) (1 4) (0 5))))
{polynomial x {2 {integer . 4}} {1 {integer . 6}} {0 {integer . 8}}}
> (mul (make-polynomial 'x '((2 1) (1 2) (0 3))) (make-polynomial 'x '((2 3) (1 4) (0 5))))
{polynomial x {4 {integer . 3}} {3 {integer . 10}} {2 {integer . 22}} {1 {integer . 22}} {0 {integer . 15}}}
>

Exercise 2.87

=zero?の実装。多項式に係数が０の項は絶対に含まれない事を前提とするなら式に一つも項がなければ０と言う事になる。

(define (install-polynomial-package)
…
  (put '=zero? '(polynomial)
       (lambda (p) (empty-termlist? (term-list p))))
  'done)

動作確認

> (=zero? (make-polynomial 'x '()))
#t
> (add (make-polynomial 'x (list (list '1 (make-polynomial 'y '((1 1) (0 2)))) (list 0 (make-polynomial 'y '((1 2))))))
       (make-polynomial 'x (list (list '1 (make-polynomial 'y '((1 -1) (0 -2)))) (list 0 (make-polynomial 'y '((1 -2)))))))
{polynomial x}
>

Exercise 2.88

まず各型にnegateを実装。ここでもinteger、rationalでは実装せずrealで実装。

(define (negate x) (apply-generic 'negate x))

(define (install-real-package)
…
  (put 'negate '(real)
       (lambda (x) (tag (- x))))
  'done)

動作確認

> (negate (make-integer 1))
{integer . -1}
> (negate (make-real -1.2))
{rational 5404319552844595.0 . 4503599627370496.0}
> (negate (make-rational 2 -3))
{rational 2 . 3}
> (negate (make-rational -2 3))
{rational 2 . 3}
>

complexの所は実際にはrectangular、polarで実装するのだが、complexから総称関数のnegateを呼んでしまうとrectangular、polarでdropしようとしておかしな事になってしまうので、ここではdropしてしまうapply-generalを呼ばないようにapply-operationを直接呼ぶ。

(define (install-rectangular-package)
…
  (put 'negate '(rectangular)
       (lambda (x) (tag (make-from-real-imag (negate (real-part x)) (negate (imag-part x))))))
  'done)

(define (install-polar-package)
…
  (put 'negate '(polar)
       (lambda (x) (tag (make-from-mag-ang (magnitude x) (negate (angle x))))))
  'done)

(define (install-complex-package)
…
  (put 'negate '(complex)
       (lambda (z) (tag (apply-operation 'negate z))))
  'done)

動作確認

> (negate (make-complex-from-real-imag (make-integer 4) (make-integer 3)))
{complex rectangular {integer . -4} integer . -3}
> (negate (make-complex-from-mag-ang 5 (/ pi 6)))
{complex polar 5 rational -4716158501352293.0 . 9007199254740992.0}
>

問題の多項式の部分は、add-poly、add-termを参考に：

(define (install-polynomial-package)
…
  (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) (negate (coeff t)))
                       (neg-terms (rest-terms L))))))
…

動作確認。(y + 2)x-3をnegateすると：

> (negate (make-polynomial 'x (list (list '1 (make-polynomial 'y '((1 1) (0 2)))) '(0 -3))))
{polynomial x {1 {polynomial y {1 {integer . -1}} {0 {integer . -2}}}} {0 {integer . 3}}}
>

(-y-2)x+3。
で、本題の多項式の引き算。

(define (install-polynomial-package)
…
  (put 'sub '(polynomial polynomial)
       (lambda (p1 p2) (tag (add-poly p1 (neg-poly p2)))))
  'done)

動作確認

> (sub (make-polynomial 'x '((2 3) (1 -4) (0 5))) (make-polynomial 'x '((2 1) (1 -4) (0 3))))
{polynomial x {2 {integer . 2}} {0 {integer . 2}}}
>

(3x²-4x+5)-(x²-4x+3)=2x²+2。
取り敢えずは良いか。

Exercise 2.89

次元が密集した多項式として別の実装をすると言う事はひとつひとつの項に次数を持たない、つまりデータの量としては密集していれば効率良く圧縮される事を目指すものと思われる。また今までの流れからして極力api等は変えない方針。
とするとterm-listの実装は変えるけどtermの実装を変えなければ、term-listからtermを取り出すfirst-termとtermをterm-listに追加するadjoin-termだけ変更すれば残りはそのまま使えるのでは。
そこでterm-listの実装だが、要するに次数を省略してデータを圧縮すると言う事。各項は次数は持たず係数だけ。次数は項リストの後ろから数える（多項式（polynomial）と言うからには次数は０以上）。

(define (install-polynomial-package)
…
  (define (first-term term-list)
    (make-term (degree term-list) (car term-list)))
    (define (adjoin-term term naked-termlist)
      (if (=zero? (coeff term))
          naked-termlist
          (cons (coeff term) (fill-terms (- (order term) 1) naked-termlist))))
…
; この実装用に追加する関数
  ; for dense term list
  (define (degree term-list)
    (if (empty-termlist? term-list)
        -1
        (length (cdr term-list))))
  (define (fill-terms order term-list)
    (if (<= order (degree term-list))
        term-list
        (fill-terms order (cons 0 term-list))))

実行結果

> (negate (make-polynomial 'x '(3 -4 5)))
{polynomial x {integer . -3} {integer . 4} {integer . -5}}
> (sub (make-polynomial 'x '(3 -4 5)) (make-polynomial 'x '(1 -4 3)))
{polynomial x {integer . 2} {integer . 0} {integer . 2}}
>

タグ付きになっちゃうのであまり圧縮は出来ていないけど。

Exercise 2.90

二つのterm-listの実装をrectangulraとpolarの様に同居させる。考えた事は以下の様な感じ：

sparseとdenseのterm-listのパッケージを作る。
termそのものは前問同様共通なのでtermを取り出す操作と追加する操作はそれぞれのterm-listの実装によって異なる。
termを操作する関数の置き場所はトップレベル（グローバル）なのもどうかなと思うので、polynomialパッケージの中に置く。
必然的にsparseとdenseのterm-listパッケージもpolynomialパッケージの中に定義。
メソッドとしてfirst-term、rest-terms、adjoin-term、empty-termlist?を総称関数にする。
adjoin-termはtermとterm-listを引数に取るので、現状だとtermにもタグが必要になる。でもsparseとdenseで使い分ける必要は無いし、タグが付いている状況、付いていない状況の扱いが面倒になるだけなので、term-listだけを引数に取る関数に変更。返り値としてtermを引数に取る関数を返す事にする。
the-empty-termlistがsparseかdenseかどちらかのタグ付きterm-listを作る必要があるが、どちらを作るのか判断する為にterm-listを引数に取る総称関数に変更。ただしこの引数自体を使用する事はない。
term-listのパッケージ内ではタグは外されているためempty-termlist?等が使えない。「タグ無し＝裸（naked）」と言う意味でempty-naked-termlist?とthe-empty-naked-termlistを定義する。この二つはsparseでもdenseでも同じなので共通にする。

総称関数部分

(define (first-term tl) (apply-operation 'first-term tl))
(define (rest-terms tl) (apply-operation 'rest-terms tl))
(define (adjoin-term t tl) ((apply-operation 'adjoin-term tl) t))
(define (empty-termlist? tl) (apply-operation 'empty-termlist? tl))
(define (the-empty-termlist tl) (apply-operation 'the-empty-termlist tl))

polynomial-package全体は以下の通り：

(define (install-polynomial-package)
  (define (order t) (car t))
  (define (coeff t) (cadr t))
  (define (make-term o c) (list o c))
  
  (define (empty-naked-termlist? ntl) (null? ntl))
  (define (the-empty-naked-termlist) '())

  (define (install-sparse-term-list-package)
    (define (tag tl) (cons 'sparse tl))
    (define (first-term naked-termlist) (car naked-termlist))
    (define (adjoin-term term naked-termlist)
      (if (=zero? (coeff term))
          naked-termlist
          (cons term naked-termlist)))
    ;
    (put 'first-term '(sparse)
         (lambda (ntl) (first-term ntl)))
    (put 'rest-terms '(sparse)
         (lambda (ntl) (tag (cdr ntl))))
    (put 'adjoin-term '(sparse)
         (lambda (ntl) (lambda (term) (tag (adjoin-term term ntl)))))
    (put 'empty-termlist? '(sparse) empty-naked-termlist?)
    (put 'make 'sparse
         (lambda () (tag (the-empty-naked-termlist))))
    (put 'the-empty-termlist '(sparse)
         (lambda (ntl) (tag (the-empty-naked-termlist))))
    'done)

  (define (install-dense-term-list-package)
    (define (tag tl) (cons 'dense tl))
    (define (first-term naked-termlist)
      (make-term (degree naked-termlist) (car naked-termlist)))
    (define (adjoin-term term naked-termlist)
      (if (=zero? (coeff term))
          naked-termlist
          (cons (coeff term) (fill-terms (- (order term) 1) naked-termlist))))
    ; for dense term list
    (define (degree naked-termlist)
      (if (empty-naked-termlist? naked-termlist)
          -1
          (length (cdr naked-termlist))))
    (define (fill-terms order naked-termlist)
      (if (<= order (degree naked-termlist))
          naked-termlist
          (fill-terms order (cons 0 naked-termlist))))
    ;
    (put 'first-term '(dense)
         (lambda (ntl) (first-term ntl)))
    (put 'rest-terms '(dense)
         (lambda (ntl) (tag (cdr ntl))))
    (put 'empty-termlist? '(dense) empty-naked-termlist?)
    (put 'adjoin-term '(dense)
         (lambda (ntl) (lambda (term) (tag (adjoin-term term ntl)))))
    (put 'make 'dense
         (lambda () (tag (the-empty-naked-termlist))))
    (put 'the-empty-termlist '(dense)
         (lambda (ntl) (tag (the-empty-naked-termlist))))
    'done)
  ;
  (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 (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 (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-terms L1 L2)
    (if (empty-termlist? L1)
        (the-empty-termlist L1)
        (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 L)
        (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 L)
        (let ((t (first-term L)))
          (adjoin-term (make-term (order t) (negate (coeff t)))
                       (neg-terms (rest-terms L))))))

  (define (tag p) (attach-tag 'polynomial p))
  ; inteface
  (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))))
  (put '=zero? '(polynomial)
       (lambda (p) (empty-termlist? (term-list p))))
  (put 'negate '(polynomial)
       (lambda (p) (tag (neg-poly p))))
  (put 'sub '(polynomial polynomial)
       (lambda (p1 p2) (tag (add-poly p1 (neg-poly p2)))))
  (install-sparse-term-list-package)
  (install-dense-term-list-package)
  'done)

(define (make-sparse-term-list)
  ((get 'make 'sparse)))
(define (make-dense-term-list)
  ((get 'make 'dense)))
(define (make-polynomial var terms)
  ((get 'make 'polynomial) var terms))

rest-termsはsparseとdenseで中身が同じに見えるので共通化できそうだが、ここで使っているtagはsparseとdenseでは異なる関数なので引数で渡すなりしないと共通化は出来ない。dynamic scopeはこの状況で使えるのかなと思った。

ここで、単項の時にメソッドが見つからなければ見つかるまでraiseする様にしていたが、実装の途中で無限再帰に陥る事甚だしく、やってらんなくなったので、raiseは直接apply-operationを呼ばずにその手前でメソッドが存在するかをチェックする様にした。

(define (raise x)
  (let ((proc (get 'raise (type-tag x))))
    (if proc
        (apply-operation 'raise x)
        (error "No raise method for" (type-tag x)))))

簡単に動作確認。

> (adjoin-term '(2 3) (adjoin-term '(1 -4) (adjoin-term '(0 5) (make-sparse-term-list))))
{sparse {2 3} {1 -4} {0 5}}
> (adjoin-term '(2 3) (adjoin-term '(1 -4) (adjoin-term '(0 5) (make-dense-term-list))))
{dense 3 -4 5}
> (sub (make-polynomial 'x (adjoin-term '(2 3) (adjoin-term '(1 -4) (adjoin-term '(0 5) (make-sparse-term-list)))))
       (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 -4) (adjoin-term '(0 3) (make-sparse-term-list))))))
{polynomial x sparse {2 {integer . 2}} {0 {integer . 2}}}
> (sub (make-polynomial 'x (adjoin-term '(2 3) (adjoin-term '(1 -4) (adjoin-term '(0 5) (make-dense-term-list)))))
       (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 -4) (adjoin-term '(0 3) (make-dense-term-list))))))
{polynomial x dense {integer . 2} {integer . 0} {integer . 2}}
>

Exercise 2.91

多項式の割り算。教科書はちょっと見にくいが $\frac{x^5-1}{x-1}=x^3+x, remainder x-1$ と書いてある。多項式の割り算は筆算の割り算と殆ど同じ。割られる数の最高次の項がが消える様に商となる多項式の各項を決めて行く。
$\hspace{42}x^3+x\\x^2-1\overline{\)x^5\hspace{53}-1}\\ \hspace{42}\underline{x^5-x^3\hspace{46}}\\ \hspace{75} x^3\hspace{22}-1\\ \hspace{75}\underline{x^3-x\hspace{22}}\\ \hspace{108}x-1$

(define (install-polynomial-package)
…
    (define (div-poly p1 p2)
    (if (same-variable? (variable p1) (variable p2))
        (make-poly (variable p1)
                   (div-terms (term-list p1)
                              (term-list p2)))
        (error "Polys not in same var -- DIV-POLY"
               (list p1 p2))))
…
  (define (div-terms L1 L2)
    (if (empty-termlist? L1)
        (list (the-empty-termlist L1) (the-empty-termlist L1)) ; need to add argument to determine sparse or dense.
        (let ((t1 (first-term L1))
              (t2 (first-term L2)))
          (if (> (order t2) (order t1))
              (list (the-empty-termlist L1) L1) ; need to add argument to determine sparse or dense.
              (let ((quotient (make-term (- (order t1) (order t2))
                                         (div (coeff t1) (coeff t2)))))
                (let ((rest-of-result (div-terms
                                       (add-terms L1 (neg-terms (mul-terms (adjoin-term quotient (the-empty-termlist L1))
                                                                           L2)))
                                       L2)))
                  (list (adjoin-term quotient (car rest-of-result)) (cadr rest-of-result))
                  ))))))
…
  (put 'div '(polynomial polynomial)
       (lambda (p1 p2) (tag (div-poly p1 p2))))
…
  'done)

div-polyはadd-polyと同じ（と言う事は纏められる事になるけど）。Exercise 2.90の結果からthe-empty-termlistには引数を追加。穴埋め問題にはなっているがnew-c、new-oと言うよりはその項が欲しいのでmake-termで項を作ってしまい、実装を少し問題文からは変える。その項で演算（mul-terms）する為にはterm-listにする必要があるのでadjoin-termで単項のterm-listに変換している。

実行結果

> (div (make-polynomial 'x (adjoin-term '(5 1) (adjoin-term '(0 -1) (make-sparse-term-list))))
       (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(0 -1) (make-sparse-term-list)))))
{polynomial x {sparse {3 {integer . 1}} {1 {integer . 1}}} {sparse {1 {integer . 1}} {0 -1}}}
> (div (make-polynomial 'x (adjoin-term '(5 1) (adjoin-term '(0 -1) (make-dense-term-list))))
       (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(0 -1) (make-dense-term-list)))))
{polynomial x {dense {integer . 1} 0 {integer . 1} 0} {dense {integer . 1} {integer . -1}}}
>

Hierarchies of types in symbolic algebra

記号代数に置ける型の階層。
xの多項式の係数にはyの多項式を含む事が出来るし、yの多項式の係数にはxの多項式を含む事が出来る。本来どちらが優先と言う事はないが、取り合えず変数に優先度を付ける事で多項式間で優先順位を付ける事が出来る。yの多項式をxの多項式に変換する為には式を展開してxで纏め直す作業が必要。

Exercise 2.92

話の流れから例えばyの多項式を係数として含むxの多項式とxの多項式を係数として含むyの多項式の足し算と掛け算を実装するものとする。xの方がyよりも優先度が高いとすると、xの多項式とyの多項式を渡された時（従来はadd-poly、mul-polyでエラーにしていたケース）にyの多項式をraiseしてxの多項式に変換すると言う事。

その前に、係数が別の変数の多項式を含む場合、その足し算や掛け算では係数同士の演算が必要となり、普通の定数と多項式の演算が出来なければならない。とすると定数は定数項のみの多項式（単項式？）にraiseできる必要があるし、定数項のみとなった多項式は定数にdrop出来なければならない。

定数から多項式へraiseする時に変数をどう決めるかが問題となる。ここでは変数名を表すシンボルの代わりにnil（Raketではnull）を使う事にする。complexからpolynomialへraiseを実装。従来と同じ方法でcomplexパッケージとcoercionパッケージにそれぞれメソッドを追加。coericonパッケージからpolynomialパッケージの中にはアクセス出来ないので、ここでmake-termはグローバルに出す。

(define (install-complex-package)
…
  (put 'raise '(complex)
       (lambda (x) ((get 'raise 'complex) x)))
…
  'done)

…
(define (make-term o c) (list o c))
…

(define (install-coercion)
…
  (put 'raise 'complex
       (lambda (x) (make-polynomial null (adjoin-term (make-term 0 (make-complex-from-real-imag (my-real-part x) (my-imag-part x))) (make-sparse-term-list)))))
…

動作確認。

> (raise (make-complex-from-real-imag 2 0))
{polynomial () sparse {0 {complex rectangular 2 . 0}}}
>

定数項が整数なのにcomplexになってしまう。raiseする過程で一旦complexになるが、元々それよりも低い型の場合は多項式の定数項を一旦dropする必要がありそう。

(define (install-coercion)
…
  (put 'raise 'complex
       (lambda (x) (make-polynomial null (adjoin-term (make-term 0 (make-complex-from-real-imag (my-real-part x) (my-imag-part x))) (make-sparse-term-list)))))
…

動作確認

> (raise (make-complex-from-real-imag 2 0))
{polynomial () sparse {0 2}}
>

dropの為にprojectとequ?を実装する必要がある。まず、projectをpolynomialパッケージに実装する。内部関数として定数項を取り出すconstant-termを実装する。

(define (install-polynomial-package)
…
  (define (constant-term p)
    (define (constant-term-aug tl)
      (cond ((empty-termlist? tl) 0)
            ((= (order (first-term tl)) 0) (coeff (first-term tl)))
            (else (constant-term-aug (rest-terms tl)))))
    (constant-term-aug (term-list p)))
…
  (put 'project '(polynomial)
       (lambda (p) (constant-term p)))
…
  'done)

動作確認。

> (project (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 4) (make-sparse-term-list))))))
4
> (project (raise (make-complex-from-real-imag 1 2)))
{complex rectangular 1 . 2}
>

次にequ?の実装。

(define (install-polynomial-package)
…
  (define (equ-poly? p1 p2)
    (define (equ-terms? t1 t2)
      (cond ((and (empty-termlist? t1) (empty-termlist? t2)) #t)
            ((or (empty-termlist? t1) (empty-termlist? t2)) #f)
            ((and (equ? (order (first-term t1)) (order (first-term t2)))
                  (equ? (coeff (first-term t1)) (coeff (first-term t2)))) (equ-terms? (rest-terms t1) (rest-terms t2)))
            (else #f)))
    (if (same-variable? (variable p1) (variable p2))
        (equ-terms? (term-list p1) (term-list p2))
        #f))
…
  (put 'equ? '(polynomial polynomial)
       (lambda (p1 p2) (equ-poly? p1 p2)))
…
  'done)

まず多項式同士の比較。

> (equ? (make-polynomial 'x (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list)))) (make-polynomial 'x (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
#t
> (equ? (make-polynomial 'x (adjoin-term '(1 2) (adjoin-term (list 0 (make-polynomial 'y (adjoin-term '(1 2) (make-sparse-term-list)))) (make-sparse-term-list))))
        (make-polynomial 'x (adjoin-term '(1 2) (adjoin-term (list 0 (make-polynomial 'y (adjoin-term '(1 2) (make-sparse-term-list)))) (make-sparse-term-list)))))
#t
> (equ? (make-polynomial 'x (adjoin-term '(1 2) (make-sparse-term-list))) (make-polynomial 'y (make-sparse-term-list)))
#f
>

係数が多項式になっていても大丈夫。
次に多項式と定数を比較する事を考える。定数をraiseした時にはシンボルはnilにしていた。定数項しか持たない多項式と定数が同値の場合にequ?が#tを返す為に、same-variable?でシンボルがnilの場合は（定数項しかない事を前提として）同じシンボルとみなす事とする。

  (define (same-variable? v1 v2)
    (if (and (variable? v1) (variable? v2))
        (eq? v1 v2)
        #t))

これで以下の比較が出来る様になる。

> (equ? (make-polynomial 'x (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list)))) (make-polynomial '() (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
#t
> (equ? (make-polynomial 'x (adjoin-term '(1 2) (adjoin-term '(0 4) (make-sparse-term-list)))) (make-polynomial '() (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
#f
>

ここまではraiseされたオブジェクトを想定して比較のテストをしていたが、実際にcomplex以下からraiseして比較する為にはcoerceのテーブルにpolynomialを追加してequ?等の演算にpolynomialとそれ以外の型が渡された時に自動的にcoerceする様にしなければならない。変換に必要なメソッドは全て揃っているので、テーブルにpolynomialを追加するだけ。

(define (compare-types a b)
…
  (comp-aux a b '(polynomial complex real rational integer)))

動作確認。

> (equ? (make-polynomial 'x (adjoin-term '(0 2) (make-sparse-term-list))) 2)
#t
> (equ? (make-polynomial 'x (adjoin-term '(0 3) (make-sparse-term-list))) 2)
#f
>

dropを実行してみる。

> (drop (make-polynomial 'x (adjoin-term '(0 2) (make-sparse-term-list))))
2
> (drop (make-polynomial 'x (adjoin-term (list 0 (make-complex-from-real-imag 2 0)) (make-sparse-term-list))))
2
>

この時点で定数と多項式の演算が出来る。

> (add 2 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (make-sparse-term-list)))))
{polynomial () sparse {2 1} {1 2} {0 2}}
> (add 2 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list))))))
{polynomial () sparse {2 1} {1 2} {0 {integer . 5}}}
> (mul 2 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list))))))
{polynomial () sparse {2 {integer . 2}} {1 {integer . 4}} {0 {integer . 6}}}
>

実はここまでが下準備。これからが本題。

変数が異なる二つの多項式を扱うとして、変数の優先度までテーブルでは管理出来ない。つまりapply-genericによって変数を変換する事は出来ない。なのでpolynomialパッケージの中の演算で変数の変換を行う。

まずは変数の優先度の判定。ここでは取り敢えず多項式の変数のシンボルのアルファベットの順番に優先度が高いものとする。つまりx、y、zではxが優先。シンボルを文字（列）に変換して比較する。空の変数の場合は優先度は最低とする。

(define (compare-symbols a b)
  (cond ((and (null? a) (null? b)) '=)
        ((null? a) '>)
        ((null? b) '<)
        (else
         (let ((as (symbol->string a))
               (bs (symbol->string b)))
           (cond ((string=? as bs) '=)
                 ((string<? as bs) '<)
                 ((string>? as bs) '>))))))

動作確認

> (compare-symbols 'x 'x)
=
> (compare-symbols 'x 'y)
<
> (compare-symbols 'y 'x)
>
> (compare-symbols 'x '())
<
> (compare-symbols '() 'x)
>
> (compare-symbols '() '())
=
>

add-poly、mul-polyのelse部分に実装するのではなく、add-poly、mul-polyに渡す前に変数を変換してしまう事にする。opとしてadd-polyかmul-polyを渡し、変数を変換してから演算するメソッドとする。

(define (coerce-poly-and-apply op p1 p2)
  (let ((symbols-are (compare-symbols (variable p1) (variable p2))))
    (cond ((eq? symbols-are '=) (op p1 p2))
          ((eq? symbols-are '<) (op p1 (coerce-poly (variable p1) p2)))
          ((eq? symbols-are '>) (op (coerce-poly (variable p2) p1) p2)))))

ちょっとスタブを作る。coerce-polyは何も変更せずに受け取ったものをそのまま返す。

(define (coerce-poly v p)
  (display "coercing to ")(display v)(display " from ")(display p)(newline)
  p)
(define (variable p)
  (car p))
(define (display-polies p1 p2)
  (display p1)(newline)
  (display p2)(newline))

動作確認

> (coerce-poly-and-apply display-polies
                         (contents (make-polynomial 'x (adjoin-term '(0 3) (make-sparse-term-list))))
                         (contents (make-polynomial 'x (adjoin-term '(0 3) (make-sparse-term-list)))))
(x sparse (0 3))
(x sparse (0 3))
> (coerce-poly-and-apply display-polies
                         (contents (make-polynomial 'x (adjoin-term '(0 3) (make-sparse-term-list))))
                         (contents (make-polynomial 'y (adjoin-term '(0 3) (make-sparse-term-list)))))
coercing to x from (y sparse (0 3))
(x sparse (0 3))
(y sparse (0 3))
> (coerce-poly-and-apply display-polies
                         (contents (make-polynomial 'y (adjoin-term '(0 3) (make-sparse-term-list))))
                         (contents (make-polynomial 'x (adjoin-term '(0 3) (make-sparse-term-list)))))
coercing to x from (y sparse (0 3))
(y sparse (0 3))
(x sparse (0 3))
> (coerce-poly-and-apply display-polies
                         (contents (make-polynomial 'y (adjoin-term '(0 3) (make-sparse-term-list))))
                         (contents (make-polynomial '() (adjoin-term '(0 3) (make-sparse-term-list)))))
coercing to y from (() sparse (0 3))
(y sparse (0 3))
(() sparse (0 3))
> (coerce-poly-and-apply display-polies
                         (contents (make-polynomial '() (adjoin-term '(0 3) (make-sparse-term-list))))
                         (contents (make-polynomial 'x (adjoin-term '(0 3) (make-sparse-term-list)))))
coercing to x from (() sparse (0 3))
(() sparse (0 3))
(x sparse (0 3))
>

一応期待通りに動いている。
compare-symbols、coerce-poly、coerce-poly-and-applyをpolynomialパッケージに入れて、addとmulを以下の様に変更する。

  (put 'add '(polynomial polynomial) 
       (lambda (p1 p2) (tag (coerce-poly-and-apply add-poly p1 p2))))
  (put 'mul '(polynomial polynomial) 
       (lambda (p1 p2) (tag (coerce-poly-and-apply mul-poly p1 p2))))

動作確認。

> (add (make-polynomial 'x (adjoin-term '(0 3) (make-sparse-term-list)))
       (make-polynomial 'y (adjoin-term '(0 3) (make-sparse-term-list))))
coercing to x from (y sparse (0 3))
. . Polys not in same var -- ADD-POLY {{x sparse {0 3}} {y sparse {0 3}}}
> (mul (make-polynomial 'x (adjoin-term '(0 3) (make-sparse-term-list)))
       (make-polynomial 'y (adjoin-term '(0 3) (make-sparse-term-list))))
coercing to x from (y sparse (0 3))
. . Polys not in same var -- MUL-POLY {{x sparse {0 3}} {y sparse {0 3}}}
> (add (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
       (make-polynomial 'y (adjoin-term (list 2 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (make-sparse-term-list)))))
                                        (adjoin-term '(1 5)
                                                     (adjoin-term (list 0 (make-polynomial 'z (adjoin-term '(1 1) (adjoin-term '(0 1)
                                                                                                                               (make-sparse-term-list)))))
                                                                  (make-sparse-term-list))))))
coercing to x from (y sparse (2 (polynomial x sparse (2 1) (1 2))) (1 5) (0 (polynomial z sparse (1 1) (0 1))))
. . Polys not in same var -- ADD-POLY {{x sparse {2 1} {1 2} {0 3}} {y sparse {2 {polynomial x sparse {2 1} {1 2}}} {1 5} {0 {polynomial z sparse {1 1} {0 1}}}}}
>

変換しようとする所まで来ている事が確認出来る。
coerce-polyを実装する。考え方は

元の多項式の各項に対して
係数を取り出す(c)
変数部分を多項式の形に変換する(v)
cとvを掛け合わせる
これ変換しようとしている変数の多項式でなければ、それを定数項にした多項式を作る
最後に変換された各項を足し合わせる

  (define (coerce-poly new-v p)
    (let ((empty (the-empty-termlist (term-list p))))
      (define (coerce-term t)  ; returns polynomial (tagged)
        (let ((v (tag (make-poly (variable p) (adjoin-term (list (order t) 1) empty))))
              (c (coeff t)))
          (let ((new-c (mul c v)))
            (if (and (eq? (type-tag new-c) 'polynomial) (same-variable? (variable (contents new-c)) new-v))
                new-c  ; this is already polynomial of new-v
                (make-polynomial new-v (adjoin-term (make-term 0 new-c) empty))))))
      (define (coerce-term-list tl)  ; collects coerced terms
        (if (empty-termlist? tl)
            '()
            (cons (coerce-term (first-term tl)) (coerce-term-list (rest-terms tl)))))
      (let ((r (coerce-term-list (term-list p))))
        (contents (if (< 1 (length r))
                      (apply add r)
                      (car r))))))

色々動作確認

> (add 2 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (make-sparse-term-list)))))
{polynomial x sparse {2 1} {1 2} {0 {integer . 2}}}
> (add 2 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list))))))
{polynomial x sparse {2 1} {1 2} {0 {integer . 5}}}
> (mul 2 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list))))))
{polynomial x sparse {2 {integer . 2}} {1 {integer . 4}} {0 {integer . 6}}}
> (add (make-polynomial 'x (adjoin-term '(0 3) (make-sparse-term-list)))
       (make-polynomial 'y (adjoin-term '(0 3) (make-sparse-term-list))))
{integer . 6}
> (add (make-polynomial 'x (adjoin-term '(1 2) (make-sparse-term-list)))
       (make-polynomial 'x (adjoin-term '(2 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
{polynomial x sparse {2 2} {1 2} {0 3}}
> (mul (make-polynomial 'x (adjoin-term '(1 2) (make-sparse-term-list)))
       (make-polynomial 'x (adjoin-term '(2 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
{polynomial x sparse {3 {integer . 4}} {1 {integer . 6}}}
> (mul (make-polynomial 'x (adjoin-term '(0 2) (make-sparse-term-list)))
       (make-polynomial 'x (adjoin-term '(2 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
{polynomial x sparse {2 {integer . 4}} {0 {integer . 6}}}
> (add (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
       (make-polynomial 'y (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list))))))
{polynomial x sparse {2 1} {1 2} {0 {polynomial y sparse {2 {integer . 1}} {1 {integer . 2}} {0 {integer . 6}}}}}
> (mul (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
       (make-polynomial 'y (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list))))))
{polynomial
 x
 sparse
 {2 {polynomial y sparse {2 {integer . 1}} {1 {integer . 2}} {0 {integer . 3}}}}
 {1 {polynomial y sparse {2 {integer . 2}} {1 {integer . 4}} {0 {integer . 6}}}}
 {0 {polynomial y sparse {2 {integer . 3}} {1 {integer . 6}} {0 {integer . 9}}}}}

> (add (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
       (make-polynomial 'y (adjoin-term (list 2 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (make-sparse-term-list)))))
                                        (adjoin-term '(1 5)
                                                     (adjoin-term (list 0 (make-polynomial 'z (adjoin-term '(1 1) (adjoin-term '(0 1)
                                                                                                                               (make-sparse-term-list)))))
                                                                  (make-sparse-term-list))))))
{polynomial
 x
 sparse
 {2 {polynomial y sparse {2 {integer . 1}} {0 {integer . 1}}}}
 {1 {polynomial y sparse {2 {integer . 2}} {0 {integer . 2}}}}
 {0 {polynomial y sparse {1 {integer . 5}} {0 {polynomial z sparse {1 {integer . 1}} {0 {integer . 4}}}}}}}
> (mul (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (adjoin-term '(0 3) (make-sparse-term-list)))))
       (make-polynomial 'y (adjoin-term (list 2 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 2) (make-sparse-term-list)))))
                                        (adjoin-term '(1 5)
                                                     (adjoin-term (list 0 (make-polynomial 'z (adjoin-term '(1 1) (adjoin-term '(0 1)
                                                                                                                               (make-sparse-term-list)))))
                                                                  (make-sparse-term-list))))))
     
{polynomial
 x
 sparse
 {4 {polynomial y sparse {2 {integer . 1}}}}
 {3 {polynomial y sparse {2 {integer . 4}}}}
 {2 {polynomial y sparse {2 {integer . 7}} {1 {integer . 5}} {0 {polynomial z sparse {1 {integer . 1}} {0 {integer . 1}}}}}}
 {1 {polynomial y sparse {2 {integer . 6}} {1 {integer . 10}} {0 {polynomial z sparse {1 {integer . 2}} {0 {integer . 2}}}}}}
 {0 {polynomial y sparse {1 {integer . 15}} {0 {polynomial z sparse {1 {integer . 3}} {0 {integer . 3}}}}}}}
>

Extended exercise: Rational functions

有理関数。多項式を分子、分母に持つ。

Exercise 2.93

add-rat等はネイティブの演算子を使っているが、これをこれまで作って来た総称関数に変更。make-ratは取り敢えずgcdを呼び出さない様にする。

(define (install-rational-package)
…
  (define (make-rat n d)
    (cons n d))
  (define (add-rat x y)
    (make-rat (add (mul (numer x) (denom y))
                   (mul (numer y) (denom x)))
              (mul (denom x) (denom y))))
  (define (sub-rat x y)
    (make-rat (sub (mul (numer x) (denom y))
                   (mul (numer y) (denom x)))
              (mul (denom x) (denom y))))
  (define (mul-rat x y)
    (make-rat (mul (numer x) (numer y))
              (mul (denom x) (denom y))))
  (define (div-rat x y)
    (make-rat (mul (numer x) (denom y))
              (mul (denom x) (numer y))))
…
  'done)

次に型変換が問題となる。今までinteger→rational→realとしていたが、integer→realとし、poynomial→rationalとする。rationalをpolynomialにprojectする時には取り敢えず分子のみとする。div演算をして商だけを取る方が自然だが、Exercise 2.92の結果、div-termsの結果もdropしようとしてしまい非常にややこしい事になるので、どの道drop出来ないのであれば簡単な方法を取る事にする。

(define (install-polynomial-package)
…
  (put 'raise '(polynomial)
       (lambda (p) ((get 'raise 'polynomial) p)))
…
  'done)

(define (install-coercion)
  (put 'raise 'integer
       (lambda (x) (make-real x)))
;  deleted raise method for rational
  (put 'raise 'real
       (lambda (x) (make-complex-from-real-imag x 0)))
…
  (put 'raise 'polynomial
       (lambda (x) (make-rational (make-polynomial (variable x) (term-list x)) (make-integer 1))))
  (put 'project 'rational
       (lambda (n d) n))
  (put 'project 'real
       (lambda (r) (make-integer (round r))))
…
  'done)

(define (compare-types a b)
  (define (comp-aux a b l)
    (cond ((null? l) #f)
          ((eq? a b) '=)
          ((eq? a (car l)) '<)
          ((eq? b (car l)) '>)
          (else (comp-aux a b (cdr l)))))
  (comp-aux a b '(rational polynomial complex real integer)))

これに伴いvariableとterm-listはパッケージの外に出す。

(define (variable p) (car p))
(define (term-list p) (cdr p))

動作確認。p1、p2の定義方法が教科書とは少し変わる。

> (define p1 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(0 1) (make-sparse-term-list)))))
> (define p2 (make-polynomial 'x (adjoin-term '(3 1) (adjoin-term '(0 1) (make-sparse-term-list)))))
> (define rf (make-rational p2 p1))
> (add rf rf)
{rational
 {polynomial x sparse {5 {integer . 2}} {3 {integer . 2}} {2 {integer . 2}} {0 {integer . 2}}}
 polynomial
 x
 sparse
 {4 {integer . 1}}
 {2 {integer . 2}}
 {0 {integer . 1}}}
>

Exercise 2.94

gcdを総称関数化する。取り敢えずintegerとpolynomialに実装する。

(define (greatest-common-divisor a b) (apply-operation 'gcd a b))

(define (install-integer-package)
…
  (put 'gcd '(integer integer)
       (lambda (a b) (gcd a b)))
  'done)

(define (install-polynomial-package)
…
  (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 a b)
    (if (empty-termlist? b)
        a
        (gcd-terms b (remainder-terms a b))))
  (define (remainder-terms a b)
    (cadr (div-terms a b)))
…
  (put 'gcd '(polynomial polynomial)
       (lambda (p1 p2) (tag (gcd-poly p1 p2))))
…

動作確認

> (greatest-common-divisor 20 16)
4
> (define p1 (make-polynomial 'x (adjoin-term '(4 1) (adjoin-term '(3 -1) (adjoin-term '(2 -2) (adjoin-term '(1 2) (make-sparse-term-list)))))))
> (define p2 (make-polynomial 'x (adjoin-term '(3 1) (adjoin-term '(1 -1) (make-sparse-term-list)))))
> (greatest-common-divisor p1 p2)
{polynomial x sparse {2 {integer . -1}} {1 {integer . 1}}}
>

検算。実際に割り算をすると： $\frac{x^4-x^3-2x^2+2x}{x^2-x}=(x-1)(x^2-x)+(-x^2+x)$ 。
$\frac{x^4-x^3-2x^2+2x}{-x^2+x}=-x^2+2$ で余り無し。また $\frac{x^3-x}{-x^2+x}=-x-1$ で余り無し。

Exercise 2.95

$Q1=P1 \times P2=(x^2-2x+1)(11x^2+7)=11x^4-22x^3+18x^2-14x+7$
$Q2=P1 \times P3=(x^2-2x+1)(13x+5)=13x^3-19x^2+3x+5$
div-termsにdisplayを挿入して実行してみる。

  (define (div-terms L1 L2)
    (display "div-terms ")(newline)
    (display " L1:")(display L1)(newline)
    (display " L2:")(display L2)(newline)
    (if (empty-termlist? L1)
        (list (the-empty-termlist L1) (the-empty-termlist L1)) ; need to add argument to determine sparse or dense.
        (let ((t1 (first-term L1))
              (t2 (first-term L2)))
          (if (> (order t2) (order t1))
              (list (the-empty-termlist L1) L1) ; need to add argument to determine sparse or dense.
              (let ((quotient (make-term (- (order t1) (order t2))
                                         (div (coeff t1) (coeff t2)))))
                (display " quotioent:")(display quotient)(newline)
                (let ((rest-of-result (div-terms
                                       (add-terms L1 (neg-terms (mul-terms (adjoin-term quotient (the-empty-termlist L1))
                                                                           L2)))
                                       L2)))
                  (list (adjoin-term quotient (car rest-of-result)) (cadr rest-of-result))
                  ))))))

実行結果

> (greatest-common-divisor q1 q2)
div-terms 
 L1:(sparse (4 (integer . 11)) (3 (integer . -22)) (2 (integer . 18)) (1 (integer . -14)) (0 (integer . 7)))
 L2:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
 quotioent:(1 (integer . 11/13))
div-terms 
 L1:(sparse (3 (real . -55/13)) (2 (real . 201/13)) (1 (real . -237/13)) (0 (integer . 7)))
 L2:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
 quotioent:(0 (real . -55/169))
div-terms 
 L1:(sparse (2 (real . 1458/169)) (1 (real . -2916/169)) (0 (real . 1458/169)))
 L2:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
div-terms 
 L1:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
 L2:(sparse (2 (real . 1458/169)) (1 (real . -2916/169)) (0 (real . 1458/169)))
 quotioent:(1 (real . 2197/1458))
div-terms 
 L1:(sparse (2 (integer . 5)) (1 (integer . -10)) (0 (integer . 5)))
 L2:(sparse (2 (real . 1458/169)) (1 (real . -2916/169)) (0 (real . 1458/169)))
 quotioent:(0 (real . 845/1458))
div-terms 
 L1:(sparse)
 L2:(sparse (2 (real . 1458/169)) (1 (real . -2916/169)) (0 (real . 1458/169)))
{polynomial x sparse {2 {real . 8 106/169}} {1 {real . -17 43/169}} {0 {real . 8 106/169}}}
>

多項式の割り算をする時に、このケースではまずQ1の $x^4$ の項を消す商を計算する。 $13x^3 \times p = 11x^4$ となるpを最初に決めるので $\frac{11}{13}x$ と言う様な半端な係数の式になってしまう。

この問題を解決する為に割り算で係数が分数にならない様に予め多項式に整数を掛けておく。
割られる多項式の次数をO₁、割る多項式の次数をO₂、割る数の最初の項の係数をcとする。先の例ではO₁=4、O₂=3、c=13なので $c^{1+O_1-O_2}=13^{1+4-3}=13^2=169$ を最初にQ1に掛けておく事になる。
試してみる。

> (greatest-common-divisor (mul 169 q1) q2)
div-terms 
 L1:(sparse (4 (integer . 1859)) (3 (integer . -3718)) (2 (integer . 3042)) (1 (integer . -2366)) (0 (integer . 1183)))
 L2:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
 quotioent:(1 (integer . 143))
div-terms 
 L1:(sparse (3 (integer . -715)) (2 (integer . 2613)) (1 (integer . -3081)) (0 (integer . 1183)))
 L2:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
 quotioent:(0 (integer . -55))
div-terms 
 L1:(sparse (2 (integer . 1458)) (1 (integer . -2916)) (0 (integer . 1458)))
 L2:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
div-terms 
 L1:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
 L2:(sparse (2 (integer . 1458)) (1 (integer . -2916)) (0 (integer . 1458)))
 quotioent:(1 (integer . 13/1458))
div-terms 
 L1:(sparse (2 (integer . 5)) (1 (integer . -10)) (0 (integer . 5)))
 L2:(sparse (2 (integer . 1458)) (1 (integer . -2916)) (0 (integer . 1458)))
 quotioent:(0 (integer . 5/1458))
div-terms 
 L1:(sparse)
 L2:(sparse (2 (integer . 1458)) (1 (integer . -2916)) (0 (integer . 1458)))
{polynomial x sparse {2 {integer . 1458}} {1 {integer . -2916}} {0 {integer . 1458}}}
>

結果は確かに分数にはならないがまだP1にはほど遠い。

Exercise 2.96

a)
paeudoremainder-termsは

  (define (pseudoremainder-terms a b)
    (let ((p (first-term a))
          (q (first-term b)))
      (cadr (div-terms (mul-terms a (adjoin-term (make-term 0 (exponent (coeff q) (- (+ 1 (order p)) (order q)))) (the-empty-termlist a))) b))))

引数のaとbはterm-listのレベルで既に変数とかは剥がされた後なので総称関数のmulはこの時点では使えない。また係数にはタグ付きの整数が入って来てしまうため組み込みのexptが使えない。総称関数exponentを取り敢えずintegerにだけ実装。

(define (exponent a b) (apply-generic 'expt a b))

(define (install-integer-package)
…
  (put 'expt '(integer integer)
       (lambda (a b) (expt a b)))
  'done)

実行結果

> (greatest-common-divisor q1 q2)
div-terms 
 L1:(sparse (4 (integer . 1859)) (3 (integer . -3718)) (2 (integer . 3042)) (1 (integer . -2366)) (0 (integer . 1183)))
 L2:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
 quotient:(1 (integer . 143))
div-terms 
 L1:(sparse (3 (integer . -715)) (2 (integer . 2613)) (1 (integer . -3081)) (0 (integer . 1183)))
 L2:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
 quotient:(0 (integer . -55))
div-terms 
 L1:(sparse (2 (integer . 1458)) (1 (integer . -2916)) (0 (integer . 1458)))
 L2:(sparse (3 (integer . 13)) (2 (integer . -21)) (1 (integer . 3)) (0 (integer . 5)))
div-terms 
 L1:(sparse (3 (integer . 27634932)) (2 (integer . -44641044)) (1 (integer . 6377292)) (0 (integer . 10628820)))
 L2:(sparse (2 (integer . 1458)) (1 (integer . -2916)) (0 (integer . 1458)))
 quotient:(1 (integer . 18954))
div-terms 
 L1:(sparse (2 (integer . 10628820)) (1 (integer . -21257640)) (0 (integer . 10628820)))
 L2:(sparse (2 (integer . 1458)) (1 (integer . -2916)) (0 (integer . 1458)))
 quotient:(0 (integer . 7290))
div-terms 
 L1:(sparse)
 L2:(sparse (2 (integer . 1458)) (1 (integer . -2916)) (0 (integer . 1458)))
{polynomial x sparse {2 {integer . 1458}} {1 {integer . -2916}} {0 {integer . 1458}}}
>

b.
aの結果の各係数を最大公約数で割るとの事。bが空となりaをgcdの多項式として返す直前で各係数をそれらの最大公約数で割る。折角gcdを総称関数にしたが、ここでは全ての係数を一遍にgcdに渡す必要がある為使えない。また各係数はタグが付いているためcontentsでデータを取り出す。

  (define (gcd-terms a b)
    (if (empty-termlist? b)
        (let ((c (apply gcd (map (lambda (t) (contents (coeff t))) (term-list a))))) ; gcd of all coefficients
          (car (div-terms a (adjoin-term (make-term 0 c) (the-empty-termlist a)))))  ; return only quotient
        (gcd-terms b (pseudoremainder-terms a b))))

動作確認
> (greatest-common-divisor q1 q2)
{polynomial x sparse {2 {integer . 1}} {1 {integer . -2}} {0 {integer . 1}}}
>
|

これで漸くP₁に戻った。

このGCDを使って多項式の約分が出来る。GCDを求めた時と同様に分子、分母をそれぞれGCDで割るときには係数が分数にならない様に係数を掛ける。結果として分子、分母の多項式の係数は大きな整数になる。最後に分子、分母の全ての係数の最大公約数でそれぞれを割る。

Exercise 2.97

a.
reduce-termsとreduce-polyを実装しろとの事だが、これだけではテストのしようがないのでb.のreduceを先取りして実装する。
pseudoremainder-termsの実装の一部が今回また必要で、今度は余りではなく商の方なので共通のpseudo-div-termsを実装し、その余りを取るpseudoremainder-termsと商を取るpesudoquotient-termsを実装。但しintegerizing factorの計算方法が異なるので、これは引数で与える事にする。

(define (reduce a b) (apply-generic 'reduce a b))

(define (install-polynomial-package)
…
  (define (gcd-terms a b)
    (if (empty-termlist? b)
        (let ((c (apply gcd (map (lambda (t) (contents (coeff t))) (term-list a))))) ; gcd of all coefficients
          (car (div-terms a (adjoin-term (make-term 0 c) (the-empty-termlist a)))))
        (gcd-terms b (pseudoremainder-terms a b (integerizing-factor a b)))))
  (define (remainder-terms a b)
    (cadr (div-terms a b)))

  (define (pseudoremainder-terms a b i)
    (cadr (pseudo-div-terms a b i)))
  (define (pseudoquotient-terms a b i)
    (car (pseudo-div-terms a b i)))
  (define (pseudo-div-terms a b i)
    (div-terms (mul-terms a (adjoin-term (make-term 0 i) (the-empty-termlist a))) b))
  (define (integerizing-factor a b)
    (let ((p (first-term a))
          (q (first-term b)))
      (exponent (coeff q) (- (+ 1 (order p)) (order q)))))
  (define (integerizing-factor-3 a b g)
    (let ((p (first-term a))
          (q (first-term b)))
      (integerizing-factor (if (> (order p) (order q))
                              a
                              b)
                          g)))
  
  (define (reduce-poly p1 p2)
    (if (same-variable? (variable p1) (variable p2))
        (let ((l (reduce-terms (term-list p1)
                               (term-list p2))))
              (list (make-poly (variable p1) (car l))
                    (make-poly (variable p2) (cadr l))))
        (error "Polys not in same var -- REDUCE-POLY"
               (list p1 p2))))   
  (define (reduce-terms a b)
    (let ((c (gcd-terms a b)))
      (let ((i (integerizing-factor-3 a b c)))
        (list (pseudoquotient-terms a c i) (pseudoquotient-terms b c i)))))

…
  (put 'reduce '(polynomial polynomial)
       (lambda (p1 p2)
         (let ((l (reduce-poly p1 p2)))
           (list (tag (car l)) (tag (cadr l))))))
…
  'done)

Exercise 2.95の例で実行すると

> (define p1 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 -2) (adjoin-term '(0 1) (make-sparse-term-list))))))
> (define p2 (make-polynomial 'x (adjoin-term '(2 11) (adjoin-term '(0 7) (make-sparse-term-list)))))
> (define p3 (make-polynomial 'x (adjoin-term '(1 13) (adjoin-term '(0 5) (make-sparse-term-list)))))
> (define q1 (mul p1 p2))
> (define q2 (mul p1 p3))
> q1
{polynomial x sparse {4 {integer . 11}} {3 {integer . -22}} {2 {integer . 18}} {1 {integer . -14}} {0 {integer . 7}}}
> q2
{polynomial x sparse {3 {integer . 13}} {2 {integer . -21}} {1 {integer . 3}} {0 {integer . 5}}}
> (greatest-common-divisor q1 q2)
{polynomial x sparse {2 {integer . 1}} {1 {integer . -2}} {0 {integer . 1}}}
> (reduce q1 q2)
{{polynomial x sparse {2 {integer . 11}} {0 {integer . 7}}} {polynomial x sparse {1 {integer . 13}} {0 {integer . 5}}}}
>

GCDとしてp1が抽出されていて、reduceによってp2とp3が復元されている。

b.
reduceをintegerにも実装してmake-rationalが整数と多項式と両方で機能する様にする。また、元々分子と分母をconsで繋いでいたが、listで繋ぐ様にしたので、denomの変更が必要。

(define (install-integer-package)
…
(define (reduce-integers n d)
(let *1
(define (make-rat n d)
(reduce n d))
…
|

整数とExercise 2.95の例で動作確認。

> (make-rational 6 4)
{rational 3 2}
> (define p1 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(1 -2) (adjoin-term '(0 1) (make-sparse-term-list))))))
> (define p2 (make-polynomial 'x (adjoin-term '(2 11) (adjoin-term '(0 7) (make-sparse-term-list)))))
> (define p3 (make-polynomial 'x (adjoin-term '(1 13) (adjoin-term '(0 5) (make-sparse-term-list)))))
> (define q1 (mul p1 p2))
> (define q2 (mul p1 p3))
> (make-rational q1 q2)
{rational {polynomial x sparse {2 {integer . 11}} {0 {integer . 7}}} {polynomial x sparse {1 {integer . 13}} {0 {integer . 5}}}}
>

整数については型変換がややこしい事になっているので、これ以上の事はここではしない。q1/q2がちゃんとp2/p3に約分されている。同じ変数を使うのでもう一度走らせて：

> (define p1 (make-polynomial 'x (adjoin-term '(1 1) (adjoin-term '(0 1) (make-sparse-term-list)))))
> (define p2 (make-polynomial 'x (adjoin-term '(3 1) (adjoin-term '(0 -1) (make-sparse-term-list)))))
> (define p3 (make-polynomial 'x (adjoin-term '(1 1) (make-sparse-term-list))))
> (define p4 (make-polynomial 'x (adjoin-term '(2 1) (adjoin-term '(0 -1) (make-sparse-term-list)))))
> (define rf1 (make-rational p1 p2))
> (define rf2 (make-rational p3 p4))
> (add rf1 rf2)
{rational
 {polynomial x sparse {3 {integer . -1}} {2 {integer . -2}} {1 {integer . -3}} {0 {integer . -1}}}
 {polynomial x sparse {4 {integer . -1}} {3 {integer . -1}} {1 {integer . 1}} {0 {integer . 1}}}}
>

確認
$\frac{x+1}{x^3-1}+\frac{x}{x^2-1}=\frac{x^4+x^3+x^2-2x-1}{x^5-x^3-x^2+1}$
GCDは

> (greatest-common-divisor (add (mul p1 p4) (mul p2 p3)) (mul p2 p4))
{polynomial x sparse {1 {integer . -1}} {0 {integer . 1}}}
>

と言う事なので（これを信用するとして）、割り算をすると
$\frac{x^4+x^3+x^2-2x-1}{-x+1}=-x^3-2x^2-3x-1$
$\frac{x^5-x^3-x^2+1}{-x+1}=-x^4-x^3+x+1$
と割り切れるので約分としては合っていそう。
GCDの計算は信用するとして、これ以上約分出来ない事を確認する。

> (define a (make-polynomial 'x (adjoin-term '(3 -1) (adjoin-term '(2 -2) (adjoin-term '(1 -3) (adjoin-term '(0 -1) (make-sparse-term-list)))))))
> a
{polynomial x sparse {3 -1} {2 -2} {1 -3} {0 -1}}
> (define b (make-polynomial 'x (adjoin-term '(4 -1) (adjoin-term '(3 -1) (adjoin-term '(1 1) (adjoin-term '(0 1) (make-sparse-term-list)))))))
> b
{polynomial x sparse {4 -1} {3 -1} {1 1} {0 1}}
> (greatest-common-divisor a b)
{polynomial x sparse {0 {integer . 1}}}
>

GCDは1なのでこれ以上は無理との事。

多項式のGCDの計算はこの手の計算の肝となる部分で、効率的な計算方法が研究されているとの事。

*1:g (gcd n d))) (list (/ n g) (/ d g)))) … (put 'reduce '(integer integer) (lambda (a b) (reduce-integers a b))) 'done) (define (install-rational-package) … (define (denom x) (cadr x