SICP問題2.96

a. 擬除算(pseudodivision), 擬剰余(pseudoremainder)を用いた修正。擬除算, 擬剰余は

PとQを多項式とし、O1をPの次数(Pの最高項の次数),O2をQの次数とする。cをQの先頭の係数とする。Pに整数化因子(integerizing factor) c^(1+O1-O2)を掛けると, 結果の多項式はdiv-termsアルゴリズムを使うと分数を用いることなくQで割り切ることが示せる.被除数にこの定数を掛け次に割る計算はPのQによる擬除算(pseudodivision)ということがあり、この除算の剰余を擬剰余(pseudoremainder)という

(計算機プログラムの構造と解釈 第二版 P124より)
ということらしいです。良く分かりませんが。

expの汎用演算子

(define (exp x y) (apply-generic 'exp x y))

scheme-number-packageへのexpの組み込み

(define (install-scheme-number-pacakge)
  (define (tag x)
    (attach-tag 'scheme-number x))
  (put 'add '(scheme-number scheme-number)
       (lambda (x y) (tag (+ x y))))
  (put 'sub '(scheme-number scheme-number)
       (lambda (x y) (tag (- x y))))
  (put 'mul '(scheme-number scheme-number)
       (lambda (x y) (tag (* x y))))
  (put 'div '(scheme-number scheme-number)
       (lambda (x y) (tag (/ x y))))
  (put 'make 'scheme-number
       (lambda (x) (tag x)))
  (put '=zero? '(scheme-number) 
       (lambda (x) (= x 0)))
  (put 'negate '(scheme-number)
       (lambda (x) (tag (- x))))
  (put 'greatest-common-divisor '(scheme-number scheme-number)
       (lambda (x y) (tag (gcd x y))))
  (put 'exp '(scheme-number scheme-number)
       (lambda (x y) (tag (expt x y))))
  'done)
(define (make-scheme-number n)
  ((get 'make 'scheme-number) n))
(install-scheme-number-pacakge)

; pseudoremainder-terms実装を追加
(define (install-polynomial-package)
  (define (make-poly variable termlist)
    (cons variable termlist))
  (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-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)
        (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 (add (order t1) (order t2))
                       (mul (coeff t1) (coeff t2)))
            (mul-term-by-all-terms t1 (rest-terms L))))))
  (define (adjoin-term term termlist)
    (if (=zero? (coeff term))
        termlist
        (cons term termlist)))
  (define (the-empty-termlist) '())
  (define (first-term termlist) (car termlist))
  (define (rest-terms termlist) (cdr termlist))
  (define (empty-termlist? termlist) (null? termlist))
  (define (make-term order coeff) (list order coeff))
  (define (order term) (car term))
  (define (coeff term) (cadr term))
  (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 L1 
                                   (negate-termlist (mul-terms (list (make-term new-o new-c)) L2)))
                        L2)
                       ))
                  (list (add-terms (list (make-term new-o new-c))
                                   (car rest-of-result))
                        (cadr rest-of-result))
                  ))))))
  (define (gcd-terms a b)
    ; トレースするために以下を追加
    ; (display (list 'parameters a b))
    ; (newline)
    (if (empty-termlist? b)
        a
;         (gcd-terms b (remainder-terms a b))))
        (gcd-terms b (pseudoremainder-terms a b))))
  (define (remainder-terms a b)
    (cadr (div-terms a b)))
  (define (pseudoremainder-terms a b)
    (let ((O1 (order (first-term a)))
          (O2 (order (first-term b)))
          (C (coeff (first-term b))))
      (let ((PI (mul-terms a (list (make-term 0 (exp C (add (sub O1 O2) 1)))))))
        (cadr (div-terms PI b)))))
  (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 (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 (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 (=zero-poly? L)
    (define (=zero-term? x)
      (or (empty-termlist? x)
          (if (=zero? (coeff (first-term x)))
              (=zero-term? (rest-terms x))
              #f)))
    (=zero-term? (term-list L)))
  (define (negate-termlist L)
    (if (empty-termlist? L)
        L
        (let ((f (first-term L))
              (r (rest-terms L)))
          (adjoin-term
           (make-term (order f)
                      (negate (coeff f)))
           (negate-termlist r)))))
  (define (sub-poly x y)
    (add-poly x 
              (make-poly (variable y)
                         (negate-termlist (term-list y)))))
  ;; システムの他の部分とのインタフェース
  (define (tag p) (attach-tag 'polynomial 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))))
  (put '=zero? '(polynomial)
       (lambda (x) (=zero-poly? x)))
  (put 'negate '(polynomial)
       (lambda (x) (tag (make-poly (variable x)
                                   (negate-termlist (term-list x))))))
  (put 'sub '(polynomial polynomial)
       (lambda (p1 p2) (tag (sub-poly p1 p2))))
  (put 'div '(polynomial polynomial)
       (lambda (p1 p2) (tag (div-poly p1 p2))))
  (put 'greatest-common-divisor '(polynomial polynomial)
       (lambda (x y) (tag (gcd-poly x y))))
  'done)
(install-polynomial-package)
(define (make-polynomial var terms)
  ((get 'make 'polynomial) var terms))

9.95と同じ内容で確認

(define p1 (make-polynomial 'x '((2 1) (1 -2) (0 1))))
; p1
(define p2 (make-polynomial 'x '((2 11) (0 7))))
; p2
(define p3 (make-polynomial 'x '((1 13) (0 5))))
; p3
(define q1 (mul p1 p2))
; q1
(define q2 (mul p1 p3))
; q2
(greatest-common-divisor q1 q2)
; (polynomial x (2 1458) (1 -2916) (0 1458))
; ※ 2.95ではは (polynomial x (2 1458/169) (1 -2916/169) (0 1458/169))

OK

b. gcd-termsを変更し、全係数を(整数の)最大公約数で割ることで、答えの整数から共通因子を除く修正

(define (install-polynomial-package)
  (define (make-poly variable termlist)
    (cons variable termlist))
  (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-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)
        (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 (add (order t1) (order t2))
                       (mul (coeff t1) (coeff t2)))
            (mul-term-by-all-terms t1 (rest-terms L))))))
  (define (adjoin-term term termlist)
    (if (=zero? (coeff term))
        termlist
        (cons term termlist)))
  (define (the-empty-termlist) '())
  (define (first-term termlist) (car termlist))
  (define (rest-terms termlist) (cdr termlist))
  (define (empty-termlist? termlist) (null? termlist))
  (define (make-term order coeff) (list order coeff))
  (define (order term) (car term))
  (define (coeff term) (cadr term))
  (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 L1 
                                   (negate-termlist (mul-terms (list (make-term new-o new-c)) L2)))
                        L2)
                       ))
                  (list (add-terms (list (make-term new-o new-c))
                                   (car rest-of-result))
                        (cadr rest-of-result))
                  ))))))
  ; 修正
  (define (gcd-terms a b)
    (if (empty-termlist? b)
        (reduce-coeff a)
        (gcd-terms b (pseudoremainder-terms a b))))
  (define (remainder-terms a b)
    (cadr (div-terms a b)))
  (define (pseudoremainder-terms a b)
    (let ((O1 (order (first-term a)))
          (O2 (order (first-term b)))
          (C (coeff (first-term b))))
      (let ((PI (mul-terms a (list (make-term 0 (exp C (add (sub O1 O2) 1)))))))
        (cadr (div-terms PI b)))))
  ; 追加
  (define (gcd-list L)
    (define (iter-gcd tmp-gcd rest-list)
      (if (null? rest-list)
          tmp-gcd
          (iter-gcd (gcd tmp-gcd (car rest-list)) (cdr rest-list))))
    (iter-gcd (car L) L))
  ; 追加
  (define (reduce-coeff a)
    (let ((coeff-gcd (gcd-list (map coeff a))))
      (define (iter-reduce reduced rest)
        (if (null? rest)
            reduced
            (iter-reduce
              (adjoin-term (make-term (order (first-term rest))
                                      (div (coeff (first-term rest)) coeff-gcd))
                           reduced)
              (cdr rest))))
      ; 2009/12/01 修正 reverse を追加
      (reverse (iter-reduce '() a))))
  (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 (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 (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 (=zero-poly? L)
    (define (=zero-term? x)
      (or (empty-termlist? x)
          (if (=zero? (coeff (first-term x)))
              (=zero-term? (rest-terms x))
              #f)))
    (=zero-term? (term-list L)))
  (define (negate-termlist L)
    (if (empty-termlist? L)
        L
        (let ((f (first-term L))
              (r (rest-terms L)))
          (adjoin-term
           (make-term (order f)
                      (negate (coeff f)))
           (negate-termlist r)))))
  (define (sub-poly x y)
    (add-poly x 
              (make-poly (variable y)
                         (negate-termlist (term-list y)))))
  ;; システムの他の部分とのインタフェース
  (define (tag p) (attach-tag 'polynomial 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))))
  (put '=zero? '(polynomial)
       (lambda (x) (=zero-poly? x)))
  (put 'negate '(polynomial)
       (lambda (x) (tag (make-poly (variable x)
                                   (negate-termlist (term-list x))))))
  (put 'sub '(polynomial polynomial)
       (lambda (p1 p2) (tag (sub-poly p1 p2))))
  (put 'div '(polynomial polynomial)
       (lambda (p1 p2) (tag (div-poly p1 p2))))
  (put 'greatest-common-divisor '(polynomial polynomial)
       (lambda (x y) (tag (gcd-poly x y))))
  'done)
(install-polynomial-package)
(define (make-polynomial var terms)
  ((get 'make 'polynomial) var terms))

テスト

(greatest-common-divisor q1 q2)
; (polynomial x (0 1) (1 -2) (2 1))
; ※ a. では(polynomial x (2 1458) (1 -2916) (0 1458))

OK
2009/12/01追加
OKじゃない。並びが逆になっているのでreduce-coeffにreverseを追加。

(greatest-common-divisor q1 q2)
; (polynomial x (2 1) (1 -2) (0 1))
; ※ a. では(polynomial x (2 1458) (1 -2916) (0 1458))

これでOK