積分のフーリエ変換

以下の積分
 \int_{-\infty}^{t}s(t')dt'
フーリエ変換を求めます。
(※以下に述べることは、斉藤洋一著『信号とシステム (コロナ社)』に書いてあります。)

方針としては、この積分 s(t)とヘビサイド関数の畳み込みなので、ヘビサイド関数のフーリエ変換を求めることに帰着します。(畳み込みのフーリエ変換フーリエ変換の積なので)

ビサイド関数は符号関数から求められます。 \theta(t) = \frac{sgn(t) + 1}{2}です。
ビサイド関数 \theta(t)

符号関数 sgn(t)

ビサイド関数のフーリエ変換を求めるために符号関数のフーリエ変換を求めます。
符号関数は次の関数 g(t)の極限をとったものです。

g(t) = \left\{\begin{array}{l}e^{-at}\theta(t)\,\,\,(t\geq0) \\-e^{at}\theta(-t)\,\,\,(t<0)\end{array}\right.

 \mathcal{F}[sgn(t)] = \lim_{a\to 0}\mathcal{F}[g(t)] = \frac{1}{j\pi f}

 \mathcal{F}[\theta(t)] = \mathcal{F}[\frac{sgn(t) + 1}{2}]=\frac{1}{2}\mathcal{F}[sgn(t)] + \frac{1}{2}\mathcal{F}[1]=\frac{1}{j 2\pi f} + \frac{1}{2}\delta(f)

さて、ヘビサイド関数のフーリエ変換が求まりました。これを使って積分フーリエ変換を求めます。積分は畳み込みです。
 \int_{-\infty}^{t}s(t')dt'=\int_{-\infty}^{\infty}s(t')\theta(t-t')dt' = s(t)*\theta(t)
畳み込みのフーリエ変換フーリエ変換の積です。
 \mathcal{F}[\int_{-\infty}^{t}s(t')dt'] = \mathcal{F}[s(t)*\theta(t)] = \mathcal{F}[s(t)]\mathcal{F}[\theta(t)]=\frac{1}{j 2\pi f}S(f) + \frac{1}{2}S(0)\delta(f)

プログラミングHaskell 7章問題9

import Data.Char

type Bit = Int

bin2int :: [Bit] -> Int
bin2int = foldr (\x y -> x + 2*y) 0

int2bin :: Int -> [Bit]
int2bin 0 = []
int2bin n = n `mod` 2 : int2bin (n `div` 2)

make8 :: [Bit] -> [Bit]
make8 bits = take 8 (bits ++ repeat 0)

encode :: String -> [Bit]
encode = concat . map (add_parity . make8 . int2bin . ord)

chop9 :: [Bit] -> [[Bit]]
chop9 [] = []
chop9 bits = take 9 bits : chop9 (drop 9 bits)

correct :: [Bit] -> Bool
correct = even . parity

check_parity :: [Bit] -> [Bit]
check_parity bits | correct bits = tail bits
                  | otherwise = error "parity is incorrect -- check_parity"

decode :: [Bit] -> String
decode = map (chr . bin2int . check_parity) . chop9

transmit :: String -> String
transmit = decode . channel . encode

channel :: [Bit] -> [Bit]
channel = tail

parity :: [Bit] -> Bit
parity bits | odd (sum bits) = 1
            | otherwise = 0

add_parity :: [Bit] -> [Bit]
add_parity bits = parity bits : bits

プログラミングHaskell 7章問題7

type Bit = Int

unfold :: (a -> Bool) -> (a -> b) -> (a -> a) -> a -> [b]
unfold p h t x | p x = []
               | otherwise = h x : unfold p h t (t x)

chop8 :: [Bit] -> [[Bit]]
chop8 = unfold (\x -> null x) (take 8) (drop 8)

map' :: (a -> b) -> [a] -> [b]
map' f = unfold (\x -> null x) (f . head) tail

iterate' :: (a -> a) -> a -> [a]
iterate' f = unfold (\x -> False) (\x -> x) f

プログラミングHakell 5章問題8

import Data.Char

-- caeser encription (upper case version)

let2int :: Char -> Int
let2int c = ord c - ord 'a'

int2let :: Int -> Char
int2let n = chr (ord 'a' + n)

shift :: Int -> Char -> Char
shift n c | isLower c = int2let ((let2int c + n) `mod` 26)
          | otherwise = toUpper (int2let ((let2int (toLower c) + n) `mod` 26))

encode :: Int -> String -> String
encode n xs = [shift n x | x <- xs]

table :: [Float]
table = [8.2, 1.5, 2.8, 4.3, 12.7, 2.2, 2.0, 6.1, 7.0, 0.2, 0.8, 4.0, 2.4,
         6.7, 7.5, 1.9, 0.1, 6.0, 6.3, 9.1, 2.8, 1.0, 2.4, 0.2, 2.0, 0.1]

percent :: Int -> Int -> Float
percent n m = (fromIntegral n / fromIntegral m) * 100

count :: Eq a => a -> [a] -> Int
count x xs = length [x'| x' <- xs, x' == x]

lowers :: String -> Int
lowers xs = length [x | x <- xs, isLower x]

freqs :: String -> [Float]
freqs xs = [percent (count x (map toLower xs)) n | x <- ['a' .. 'z']]
  where
    n = length xs

chisqr :: [Float] -> [Float] -> Float
chisqr os es = sum [(o - e)^2 / e | (o, e) <- zip os es]

rotate :: Int -> [a] -> [a]
rotate n xs = drop n xs ++ take n xs

positions :: Eq a => a -> [a] -> [Int]
positions x xs = [i | (i, x') <- zip [0 .. n] xs, x' == x]
  where
    n = length xs - 1

crack :: String -> String
crack xs = encode (-factor) xs
  where
    factor = head (positions (minimum chitab) chitab)
    chitab = [chisqr (rotate n table') table | n <- [0 .. 25]]
    table' = freqs xs

SICP問題2.87, 2.88 多項式どうしの引き算

;; 2.87, 2.88
;;;;;;;;;;;;;
;; 汎用演算 ;;
;;;;;;;;;;;;;
(define (equ? x y) (apply-generic 'equ? x y))
(define (project x) (apply-generic 'project x))
(define (raise x) (apply-generic 'raise x))
(define (add x y) (apply-generic 'add x y))
(define (sub x y) (apply-generic 'sub x y))
(define (mul x y) (apply-generic 'mul x y))
(define (div x y) (apply-generic 'div x y))
(define (real-part z) (apply-generic 'real-part z))
(define (imag-part z) (apply-generic 'imag-part z))
(define (magnitude z) (apply-generic 'magnitude z))
(define (angle z) (apply-generic 'angle z))
(define (my-sqrt x) (apply-generic 'my-sqrt x))
(define (my-atan x y) (apply-generic 'my-atan x y))
(define (cosine x) (apply-generic 'cosine x))
(define (sine x) (apply-generic 'sine x))
(define (=zero? x) (apply-generic '=zero? x))

;;;;;;;;;;;;;
;; 数の生成 ;;
;;;;;;;;;;;;;
(define (make-integer n)
  ((get 'make 'integer) n))

(define (make-real n)
  ((get 'make 'real) n))

(define (make-rational n d)
  ((get 'make 'rational) n d))

(define (make-complex-from-real-imag x y)
  ((get 'make-from-real-imag 'complex) x y))

(define (make-complex-from-mag-ang r a)
  ((get 'make-from-mag-ang 'complex) r a))

(define (make-polynomial var terms)
  ((get 'make 'polynomial) var terms))

;;;;;;;;;;;;;;;;;;;
;; apply-generic ;;
;;;;;;;;;;;;;;;;;;;
(define (apply-generic op . args)
  (let ((type-tags (map type-tag args)))
    (let ((proc (get op type-tags)))
      (if proc
          (if (or (eq? op 'project) (eq? op 'raise) (eq? op 'equ?) (eq? op '=zero?))
              (apply proc (map contents args))
              (drop (apply proc (map contents args))))
          (if (= (length args) 2)
              (let ((a1 (car args))
                    (a2 (cadr args)))
                (cond [(> (type-height a1) (type-height a2))
                       (apply-generic op (raise a1) a2)]
                      [(< (type-height a1) (type-height a2))
                       (apply-generic op a1 (raise a2))]
                      [else (error "No method for these types"
                                   (list op type-tags))]))
              (error "No method for these types"
                     (list op type-tags)))))))

;;;;;;;;;;;;;;
;; table操作 ;;
;;;;;;;;;;;;;;
(define table '())

(define (put func type contents)
  (set! table (cons (list func type contents) table)))

(define (get func type)
  (letrec ((get-sub (lambda (func type table)
                      (cond [(null? table) #f]
                            [(and (equal? func (caar table))
                                  (equal? type (cadar table)))
                             (caddar table)]
                            [else (get-sub func type (cdr table))]))))
    (get-sub func type table)))

;;;;;;;;;;;;;;;;
;; その他の関数 ;;
;;;;;;;;;;;;;;;;
;; drop演算
(define (drop x)
  (let ((projected (project x)))
    (if (equ? x (raise projected))
        (if (eq? (type-tag projected) 'integer)
            projected
            (drop projected))
        x)))

;; 型の塔における高さ
(define (type-height x)
  (if (eq? (type-tag x) 'polynomial)
      0
      (+ 1 (type-height (raise x)))))

;; tagの付与
(define (attach-tag type-tag contents)
  (cons type-tag contents))

;; tagの取得
(define (type-tag datum)
  (if (pair? datum)
      (car datum)
      (error "Not pair -- TYPE-TAG" datum)))

;; contentsの取得
(define (contents datum)
  (if (pair? datum)
      (cdr datum)
      (error "Not pair -- CONTENTS" datum)))

;;;;;;;;;;;;;;
;; パッケージ ;;
;;;;;;;;;;;;;;
;; 整数パッケージ
(define (install-integer-package)
  (let ((tag (lambda (x) (attach-tag 'integer x))))
    (put '=zero? '(integer)
         (lambda (x) (= x 0)))
    (put 'add '(integer integer)
         (lambda (x y) (tag (+ x y))))
    (put 'sub '(integer integer)
         (lambda (x y) (tag (- x y))))
    (put 'mul '(integer integer)
         (lambda (x y) (tag (* x y))))
    (put 'div '(integer integer)
         (lambda (x y) (tag (quotient x y))))
    (put 'my-sqrt '(integer)
         (lambda (x) (make-real (sqrt x))))
    (put 'my-atan '(integer integer)
         (lambda (x y) (make-real (atan x y))))
    (put 'cosine '(integer)
         (lambda (x) (make-real (cos x))))
    (put 'sine '(integer)
         (lambda (x) (make-real (sin x))))
    (put 'equ? '(integer integer)
         (lambda (x y) (= x y)))
    (put 'make 'integer
         (lambda (x) (tag (floor x))))
    (put 'project '(integer)
         (lambda (x) (tag x)))
    (put 'raise '(integer)
         (lambda (x) (make-rational x 1)))
    'done))

;; 有理数パッケージ
(define (install-rational-package)
  (letrec* ((numer (lambda (x) (car x)))
            (denom (lambda (x) (cdr x)))
            (make-rat (lambda (n d)
                        (let ((g (gcd n d)))
                          (cons (/ n g) (/ d g)))))
            (gcd (lambda (a b)
                   (if (zero? b)
                       a
                       (gcd b (remainder a b)))))
            (add-rat (lambda (x y)
                       (make-rat (+ (* (numer x) (denom y))
                                    (* (numer y) (denom x)))
                                 (* (denom x) (denom y)))))
            (sub-rat (lambda (x y)
                       (make-rat (- (* (numer x) (denom y))
                                    (* (numer y) (denom x)))
                                 (* (denom x) (denom y)))))
            (mul-rat (lambda (x y)
                       (make-rat (* (numer x) (numer y))
                                 (* (denom x) (denom y)))))
            (div-rat (lambda (x y)
                       (make-rat (* (numer x) (denom y))
                                 (* (denom x) (numer y)))))
            (tag (lambda (x) (attach-tag 'rational x))))
           (put '=zero? '(rational)
                (lambda (x) (= (numer x) 0)))
           (put 'add '(rational rational)
                (lambda (x y) (tag (add-rat x y))))
           (put 'sub '(rational rational)
                (lambda (x y) (tag (sub-rat x y))))
           (put 'mul '(rational rational)
                (lambda (x y) (tag (mul-rat x y))))
           (put 'div '(rational rational)
                (lambda (x y) (tag (div-rat x y))))
           (put 'my-sqrt '(rational)
                (lambda (x) (make-real (sqrt (/ (numer x) (denom x))))))
           (put 'my-atan '(rational rational)
                (lambda (x y) (make-real
                               (atan (/ (numer x) (denom x))
                                     (/ (numer y) (denom y))))))
           (put 'cosine '(rational)
                (lambda (x) (make-real
                             (cos (/ (numer x) (denom x))))))
           (put 'sine '(rational)
                (lambda (x) (make-real
                             (sin (/ (numer x) (denom x))))))
           (put 'equ? '(rational rational)
                (lambda (x y) (= (/ (numer x) (denom x))
                                 (/ (numer y) (denom y)))))
           (put 'make 'rational
                (lambda (n d) (tag (make-rat n d))))
           (put 'raise '(rational)
                (lambda (x) (make-real (/ (numer x) (denom x)))))
           (put 'project '(rational)
                (lambda (x) (make-integer (numer x))))
           'done))

;; 実数パッケージ
(define (install-real-package)
  (let ((tag (lambda (x) (attach-tag 'real x))))
    (put '=zero? '(real)
         (lambda (x) (= x 0)))
    (put 'add '(real real)
         (lambda (x y) (tag (+ x y))))
    (put 'sub '(real real)
         (lambda (x y) (tag (- x y))))
    (put 'mul '(real real)
         (lambda (x y) (tag (* x y))))
    (put 'div '(real real)
         (lambda (x y) (tag (/ x y))))
    (put 'my-sqrt '(real)
         (lambda (x) (tag (sqrt x))))
    (put 'my-atan '(real real)
         (lambda (x y) (tag (atan x y))))
    (put 'cosine '(real)
         (lambda (x) (tag (cos x))))
    (put 'sine '(real)
         (lambda (x) (tag (sin x))))
    (put 'equ? '(real real)
         (lambda (x y) (= x y)))
    (put 'make 'real
         (lambda (x) (tag x)))
    (put 'raise '(real)
         (lambda (x) (make-complex-from-real-imag
                      (make-real x) (make-real 0))))
    (put 'project '(real)
         (lambda (x)(make-rational (floor x) 1)))
    'done))

;; 複素数パッケージ
(define (install-complex-package)
  (let* ((make-from-real-imag (lambda (x y)
                                ((get 'make-from-real-imag 'rectangular) x y)))
         (make-from-mag-ang (lambda (r a)
                              ((get 'make-from-mag-ang 'polar) r a)))
         (add-complex (lambda (z1 z2)
                        (make-from-real-imag
                         (add (real-part z1) (real-part z2))
                         (add (imag-part z1) (imag-part z2)))))
         (sub-complex (lambda (z1 z2)
                        (make-from-real-imag
                         (sub (real-part z1) (real-part z2))
                         (sub (imag-part z1) (imag-part z2)))))
         (mul-complex (lambda (z1 z2)
                        (make-from-mag-ang
                         (mul (magnitude z1) (magnitude z2))
                         (add (angle z1) (angle z2)))))
         (div-complex (lambda (z1 z2)
                        (make-from-mag-ang
                         (div (magnitude z1) (magnitude z2))
                         (sub (angle z1) (angle z2)))))
         (tag (lambda (z) (attach-tag 'complex z))))
    (put '=zero? '(complex)
         (lambda (z) (and (= (real-part z) 0) (= (imag-part z) 0))))
    (put 'add '(complex complex)
         (lambda (z1 z2) (tag (add-complex z1 z2))))
    (put 'sub '(complex complex)
         (lambda (z1 z2) (tag (sub-complex z1 z2))))
    (put 'mul '(complex complex)
         (lambda (z1 z2) (tag (mul-complex z1 z2))))
    (put 'div '(complex complex)
         (lambda (z1 z2) (tag (div-complex z1 z2))))
    (put 'equ? '(complex complex)
         (lambda (z1 z2) (and (equ? (real-part z1) (real-part z2))
                              (equ? (imag-part z1) (imag-part z2)))))
    (put 'real-part '(complex) real-part)
    (put 'imag-part '(complex) imag-part)
    (put 'magnitude '(complex) magnitude)
    (put 'angle '(complex) angle)
    (put 'project '(complex)
         (lambda (z) (real-part z)))
    (put 'raise '(complex)
         (lambda (z) (make-polynomial 'x (list (list 0 (tag z))))))
    (put 'make-from-real-imag 'complex
         (lambda (x y) (tag (make-from-real-imag x y))))
    (put 'make-from-mag-ang 'complex
         (lambda (r a) (tag (make-from-mag-ang r a))))
    'done))

;; 直交座標パッケージ
(define (install-rectangular-package)
  (let* ((real-part (lambda (z) (car z)))
         (imag-part (lambda (z) (cdr z)))
         (make-from-real-imag (lambda (x y) (cons x y)))
         (magnitude (lambda (z) (my-sqrt (add (mul (real-part z) (real-part z))
                                              (mul (imag-part z) (imag-part z))))))
         (angle (lambda (z) (my-atan (imag-part z) (real-part z))))
         (make-from-mag-ang (lambda (r a) (cons (mul r (cosine a))
                                                (mul r (sine a)))))
         (tag (lambda (x) (attach-tag 'rectangular x))))
    (put 'real-part '(rectangular) real-part)
    (put 'imag-part '(rectangular) imag-part)
    (put 'magnitude '(rectangular) magnitude)
    (put 'angle '(rectangular) angle)
    (put 'make-from-real-imag 'rectangular
         (lambda (x y) (tag (make-from-real-imag x y))))
    (put 'make-from-mag-ang 'rectangular
         (lambda (r a) (tag (make-from-mag-ang r a))))
    'done))

;; 極座標パッケージ
(define (install-polar-package)
  (let* ((magnitude (lambda (z) (car z)))
         (angle (lambda (z) (cdr z)))
         (make-from-mag-ang (lambda (r a) (cons r a)))
         (real-part (lambda (z) (mul (magnitude z) (cosine (angle z)))))
         (imag-part (lambda (z) (mul (magnitude z) (sine (angle z)))))
         (make-from-real-imag (lambda (x y)
                                (cons (my-sqrt (add (mul x x) (mul y y)))
                                      (my-atan y x))))
         (tag (lambda (x) (attach-tag 'polar x))))
    (put 'real-part '(polar) real-part)
    (put 'imag-part '(polar) imag-part)
    (put 'magnitude '(polar) magnitude)
    (put 'angle '(polar) angle)
    (put 'make-from-real-imag 'polar
         (lambda (x y) (tag (make-from-real-imag x y))))
    (put 'make-from-mag-ang 'polar
         (lambda (r a) (tag (make-from-mag-ang r a))))
    'done))

;; 多項式パッケージ
(define (install-polynomial-package)
  (letrec* ((make-poly (lambda (variable term-list)
                         (cons variable term-list)))
            (variable (lambda (p) (car p)))
            (term-list (lambda (p) (cdr p)))
            (variable? (lambda (x) (symbol? x)))
            (same-variable? (lambda (v1 v2)
                              (and (variable? v1) (variable? v2) (eq? v1 v2))))
            (adjoin-term (lambda (term term-list)
                           (if (=zero? (coeff term))
                               term-list
                               (cons term term-list))))
            (first-term (lambda (term-list) (car term-list)))
            (rest-terms (lambda (term-list) (cdr term-list)))
            (the-empty-termlist (lambda () '()))
            (empty-termlist? (lambda (term-list) (null? term-list)))
            (make-term (lambda (order coeff) (list order coeff)))
            (order (lambda (term) (car term)))
            (coeff (lambda (term) (cadr term)))
            (eq-termlist (lambda (tl1 tl2)
                           (cond [(empty-termlist? tl1) (empty-termlist? tl2)]
                                 [(empty-termlist? tl2) (empty-termlist? tl1)]
                                 [(and (= (order (first-term tl1))
                                          (order (first-term tl2)))
                                       (= (coeff (first-term tl1))
                                          (coeff (first-term tl2))))
                                  (eq-termlist (rest-terms tl1)
                                               (rest-terms tl2))]
                                 [else #f])))
            (inverse-termlist (lambda (t)
                                (if (empty-termlist? t)
                                    t
                                    (cons (list (order (first-term t))
                                                (mul (make-integer -1) (coeff (first-term t))))
                                          (inverse-termlist (rest-terms t))))))
            (add-terms (lambda (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)))]))])))
            (mul-terms (lambda (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)))))
            (mul-term-by-all-terms (lambda (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)))))))
            (add-poly (lambda (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)))))
            (sub-poly (lambda (p1 p2)
                        (if (same-variable? (variable p1) (variable p2))
                            (make-poly (variable p1)
                                       (add-terms (term-list p1)
                                                  (inverse-termlist
                                                   (term-list p2))))
                            (error "Polys not in same var -- SUB-POLY"))))
            (mul-poly (lambda (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)))))
            (tag (lambda (p) (attach-tag 'polynomial p))))
           (put '=zero? '(polynomial)
                (lambda (p) (empty-termlist? p)))
           (put 'equ? '(polynomial polynomial)
                (lambda (p1 p2) (and (eq? (variable p1) (variable p2))
                                     (eq-termlist (term-list p1) (term-list p2)))))
           (put 'add '(polynomial polynomial)
                (lambda (p1 p2) (tag (add-poly p1 p2))))
           (put 'sub '(polynomial polynomial)
                (lambda (p1 p2) (tag (sub-poly p1 p2))))
           (put 'mul '(polynomial polynomial)
                (lambda (p1 p2) (tag (mul-poly p1 p2))))
           (put 'project '(polynomial)
                (lambda (p) (make-integer (order (first-term (term-list p))))))
           (put 'make 'polynomial
                (lambda (var terms) (tag (make-poly var terms))))
           'done))

プログラムが複雑になってきて大変になってきた。いちおう動いているみたいだけど、いろいろ変な設計になっている。

SICP問題2.86

;; 2.86
;;;;;;;;;;;;;
;; 汎用演算 ;;
;;;;;;;;;;;;;
(define (equ? x y) (apply-generic 'equ? x y))
(define (project x) (apply-generic 'project x))
(define (raise x) (apply-generic 'raise x))
(define (add x y) (apply-generic 'add x y))
(define (sub x y) (apply-generic 'sub x y))
(define (mul x y) (apply-generic 'mul x y))
(define (div x y) (apply-generic 'div x y))
(define (real-part z) (apply-generic 'real-part z))
(define (imag-part z) (apply-generic 'imag-part z))
(define (magnitude z) (apply-generic 'magnitude z))
(define (angle z) (apply-generic 'angle z))
(define (my-sqrt x) (apply-generic 'my-sqrt x))
(define (my-atan x y) (apply-generic 'my-atan x y))
(define (cosine x) (apply-generic 'cosine x))
(define (sine x) (apply-generic 'sine x))

;;;;;;;;;;;;;
;; 数の生成 ;;
;;;;;;;;;;;;;
(define (make-integer n)
  ((get 'make 'integer) n))

(define (make-real n)
  ((get 'make 'real) n))

(define (make-rational n d)
  ((get 'make 'rational) n d))

(define (make-complex-from-real-imag x y)
  ((get 'make-from-real-imag 'complex) x y))

(define (make-complex-from-mag-ang r a)
  ((get 'make-from-mag-ang 'complex) r a))

;;;;;;;;;;;;;;;;;;;
;; apply-generic ;;
;;;;;;;;;;;;;;;;;;;
(define (apply-generic op . args)
  (let ((type-tags (map type-tag args)))
    (let ((proc (get op type-tags)))
      (if proc
          (if (or (eq? op 'project) (eq? op 'raise) (eq? op 'equ?))
              (apply proc (map contents args))
              (drop (apply proc (map contents args))))
          (if (= (length args) 2)
              (let ((a1 (car args))
                    (a2 (cadr args)))
                (cond [(> (type-height a1) (type-height a2))
                       (apply-generic op (raise a1) a2)]
                      [(< (type-height a1) (type-height a2))
                       (apply-generic op a1 (raise a2))]
                      [else (error "No method for these types"
                                   (list op type-tags))]))
              (error "No method for these types"
                     (list op type-tags)))))))

;;;;;;;;;;;;;;
;; table操作 ;;
;;;;;;;;;;;;;;
(define table '())

(define (put func type contents)
  (set! table (cons (list func type contents) table)))

(define (get func type)
  (letrec ((get-sub (lambda (func type table)
                      (cond [(null? table) #f]
                            [(and (equal? func (caar table))
                                  (equal? type (cadar table)))
                             (caddar table)]
                            [else (get-sub func type (cdr table))]))))
    (get-sub func type table)))

;;;;;;;;;;;;;;;;
;; その他の関数 ;;
;;;;;;;;;;;;;;;;
;; drop演算
(define (drop x)
  (let ((projected (project x)))
    (if (equ? x (raise projected))
        (if (eq? (type-tag projected) 'integer)
            projected
            (drop projected))
        x)))

;; 型の塔における高さ
(define (type-height x)
  (if (eq? (type-tag x) 'complex)
      0
      (+ 1 (type-height (raise x)))))

;; tagの付与
(define (attach-tag type-tag contents)
  (cons type-tag contents))

;; tagの取得
(define (type-tag datum)
  (if (pair? datum)
      (car datum)
      (error "Not pair -- TYPE-TAG" datum)))

;; contentsの取得
(define (contents datum)
  (if (pair? datum)
      (cdr datum)
      (error "Not pair -- CONTENTS" datum)))

;;;;;;;;;;;;;;
;; パッケージ ;;
;;;;;;;;;;;;;;
;; 整数パッケージ
(define (install-integer-package)
  (let ((tag (lambda (x) (attach-tag 'integer x))))
    (put 'add '(integer integer)
         (lambda (x y) (tag (+ x y))))
    (put 'sub '(integer integer)
         (lambda (x y) (tag (- x y))))
    (put 'mul '(integer integer)
         (lambda (x y) (tag (* x y))))
    (put 'div '(integer integer)
         (lambda (x y) (tag (quotient x y))))
    (put 'my-sqrt '(integer)
         (lambda (x) (make-real (sqrt x))))
    (put 'my-atan '(integer integer)
         (lambda (x y) (make-real (atan x y))))
    (put 'cosine '(integer)
         (lambda (x) (make-real (cos x))))
    (put 'sine '(integer)
         (lambda (x) (make-real (sin x))))
    (put 'equ? '(integer integer)
         (lambda (x y) (= x y)))
    (put 'make 'integer
         (lambda (x) (tag (floor x))))
    (put 'project '(integer)
         (lambda (x) (tag x)))
    (put 'raise '(integer)
         (lambda (x) (make-rational x 1)))
    'done))

;; 有理数パッケージ
(define (install-rational-package)
  (letrec* ((numer (lambda (x) (car x)))
            (denom (lambda (x) (cdr x)))
            (make-rat (lambda (n d)
                        (let ((g (gcd n d)))
                          (cons (/ n g) (/ d g)))))
            (gcd (lambda (a b)
                   (if (zero? b)
                       a
                       (gcd b (remainder a b)))))
            (add-rat (lambda (x y)
                       (make-rat (+ (* (numer x) (denom y))
                                    (* (numer y) (denom x)))
                                 (* (denom x) (denom y)))))
            (sub-rat (lambda (x y)
                       (make-rat (- (* (numer x) (denom y))
                                    (* (numer y) (denom x)))
                                 (* (denom x) (denom y)))))
            (mul-rat (lambda (x y)
                       (make-rat (* (numer x) (numer y))
                                 (* (denom x) (denom y)))))
            (div-rat (lambda (x y)
                       (make-rat (* (numer x) (denom y))
                                 (* (denom x) (numer y)))))
            (tag (lambda (x) (attach-tag 'rational x))))
           (put 'add '(rational rational)
                (lambda (x y) (tag (add-rat x y))))
           (put 'sub '(rational rational)
                (lambda (x y) (tag (sub-rat x y))))
           (put 'mul '(rational rational)
                (lambda (x y) (tag (mul-rat x y))))
           (put 'div '(rational rational)
                (lambda (x y) (tag (div-rat x y))))
           (put 'my-sqrt '(rational)
                (lambda (x) (make-real (sqrt (/ (numer x) (denom x))))))
           (put 'my-atan '(rational rational)
                (lambda (x y) (make-real
                               (atan (/ (numer x) (denom x))
                                     (/ (numer y) (denom y))))))
           (put 'cosine '(rational)
                (lambda (x) (make-real
                             (cos (/ (numer x) (denom x))))))
           (put 'sine '(rational)
                (lambda (x) (make-real
                             (sin (/ (numer x) (denom x))))))
           (put 'equ? '(rational rational)
                (lambda (x y) (= (/ (numer x) (denom x))
                                 (/ (numer y) (denom y)))))
           (put 'make 'rational
                (lambda (n d) (tag (make-rat n d))))
           (put 'raise '(rational)
                (lambda (x) (make-real (/ (numer x) (denom x)))))
           (put 'project '(rational)
                (lambda (x) (make-integer (numer x))))
           'done))

;; 実数パッケージ
(define (install-real-package)
  (let ((tag (lambda (x) (attach-tag 'real x))))
    (put 'add '(real real)
         (lambda (x y) (tag (+ x y))))
    (put 'sub '(real real)
         (lambda (x y) (tag (- x y))))
    (put 'mul '(real real)
         (lambda (x y) (tag (* x y))))
    (put 'div '(real real)
         (lambda (x y) (tag (/ x y))))
    (put 'my-sqrt '(real)
         (lambda (x) (tag (sqrt x))))
    (put 'my-atan '(real real)
         (lambda (x y) (tag (atan x y))))
    (put 'cosine '(real)
         (lambda (x) (tag (cos x))))
    (put 'sine '(real)
         (lambda (x) (tag (sin x))))
    (put 'equ? '(real real)
         (lambda (x y) (= x y)))
    (put 'make 'real
         (lambda (x) (tag x)))
    (put 'raise '(real)
         (lambda (x) (make-complex-from-real-imag
                      (make-real x) (make-real 0))))
    (put 'project '(real)
         (lambda (x)(make-rational (floor x) 1)))
    'done))

;; 複素数パッケージ
(define (install-complex-package)
  (let* ((make-from-real-imag (lambda (x y)
                                ((get 'make-from-real-imag 'rectangular) x y)))
         (make-from-mag-ang (lambda (r a)
                              ((get 'make-from-mag-ang 'polar) r a)))
         (add-complex (lambda (z1 z2)
                        (make-from-real-imag
                         (add (real-part z1) (real-part z2))
                         (add (imag-part z1) (imag-part z2)))))
         (sub-complex (lambda (z1 z2)
                        (make-from-real-imag
                         (sub (real-part z1) (real-part z2))
                         (sub (imag-part z1) (imag-part z2)))))
         (mul-complex (lambda (z1 z2)
                        (make-from-mag-ang
                         (mul (magnitude z1) (magnitude z2))
                         (add (angle z1) (angle z2)))))
         (div-complex (lambda (z1 z2)
                        (make-from-mag-ang
                         (div (magnitude z1) (magnitude z2))
                         (sub (angle z1) (angle z2)))))
         (tag (lambda (z) (attach-tag 'complex z))))
    (put 'add '(complex complex)
         (lambda (z1 z2) (tag (add-complex z1 z2))))
    (put 'sub '(complex complex)
         (lambda (z1 z2) (tag (sub-complex z1 z2))))
    (put 'mul '(complex complex)
         (lambda (z1 z2) (tag (mul-complex z1 z2))))
    (put 'div '(complex complex)
         (lambda (z1 z2) (tag (div-complex z1 z2))))
    (put 'equ? '(complex complex)
         (lambda (z1 z2) (and (equ? (real-part z1) (real-part z2))
                              (equ? (imag-part z1) (imag-part z2)))))
    (put 'real-part '(complex) real-part)
    (put 'imag-part '(complex) imag-part)
    (put 'magnitude '(complex) magnitude)
    (put 'angle '(complex) angle)
    (put 'project '(complex)
         (lambda (z) (real-part z)))
    (put 'make-from-real-imag 'complex
         (lambda (x y) (tag (make-from-real-imag x y))))
    (put 'make-from-mag-ang 'complex
         (lambda (r a) (tag (make-from-mag-ang r a))))
    'done))

;; 直交座標パッケージ
(define (install-rectangular-package)
  (let* ((real-part (lambda (z) (car z)))
         (imag-part (lambda (z) (cdr z)))
         (make-from-real-imag (lambda (x y) (cons x y)))
         (magnitude (lambda (z) (my-sqrt (add (mul (real-part z) (real-part z))
                                              (mul (imag-part z) (imag-part z))))))
         (angle (lambda (z) (my-atan (imag-part z) (real-part z))))
         (make-from-mag-ang (lambda (r a) (cons (mul r (cosine a))
                                                (mul r (sine a)))))
         (tag (lambda (x) (attach-tag 'rectangular x))))
    (put 'real-part '(rectangular) real-part)
    (put 'imag-part '(rectangular) imag-part)
    (put 'magnitude '(rectangular) magnitude)
    (put 'angle '(rectangular) angle)
    (put 'make-from-real-imag 'rectangular
         (lambda (x y) (tag (make-from-real-imag x y))))
    (put 'make-from-mag-ang 'rectangular
         (lambda (r a) (tag (make-from-mag-ang r a))))
    'done))

;; 極座標パッケージ
(define (install-polar-package)
  (let* ((magnitude (lambda (z) (car z)))
         (angle (lambda (z) (cdr z)))
         (make-from-mag-ang (lambda (r a) (cons r a)))
         (real-part (lambda (z) (mul (magnitude z) (cosine (angle z)))))
         (imag-part (lambda (z) (mul (magnitude z) (sine (angle z)))))
         (make-from-real-imag (lambda (x y)
                                (cons (my-sqrt (add (mul x x) (mul y y)))
                                      (my-atan y x))))
         (tag (lambda (x) (attach-tag 'polar x))))
    (put 'real-part '(polar) real-part)
    (put 'imag-part '(polar) imag-part)
    (put 'magnitude '(polar) magnitude)
    (put 'angle '(polar) angle)
    (put 'make-from-real-imag 'polar
         (lambda (x y) (tag (make-from-real-imag x y))))
    (put 'make-from-mag-ang 'polar
         (lambda (r a) (tag (make-from-mag-ang r a))))
    'done))