テレスコーピングについて

テレスコーピングとは

『アルゴリズム辞典』（共立出版）の「最良近似式計算法」より引用

多項式の形の最良近似式があって，これより低次の最良に近い近似式を次のようにして簡便に求めることができる．xⁿが n 次以下のチェビシェフ多項式の１次結合で表わされるので，もとの多項式をチェビシェフ多項式の１次結合に展開し，最高次から数値を削除して，再びベキ級数に戻す．これをテレスコーピング（telescoping）という．

この頁では、三角関数を例に、具体的に多項式の係数を求める方法を解説する。計算には Scheme を使用する（実際に使用した処理系は Gauche である）。

多項式近似

各三角関数とその微分の sin'(x) = cos(x), cos'(x) = - sin(x) という関係から、以下のように三角関数はテイラー展開される。（ ! == 階乗）

cos(x) = 1 - x²/2! + x⁴/4! - x⁶/6! + x⁸/8! ...

sin(x) = x - x³/3! + x⁵/5! - x⁷/7! + x⁹/9! ...

この多項式展開から、低次の数項を取り出したものは、0 の前後で三角関数をよく近似する近似式として利用できる。

チェビシェフ多項式

チェビシェフ多項式は T_n = cos(n*cos^-1(x)) を満たす多項式で、オンラインでは Wolfram の Mathworld に詳しい解説がある（ http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html ）。

プログラムで具体的に係数を求めたい場合は以下の漸化式が利用できる。

T₀(x) = 1

T₁(x) = x

T_n+1(x) = 2*x*T_n(x) - T_n-1(x)

関数近似の他には、信号処理のフィルタなどでも使われる（wikipedia:チェビシェフフィルタ等を参照）。

テレスコーピングの実際

まず、簡単なところで、sin の 5 次の近似 x - x³/3! + x⁵/5! のテレスコーピングをやってみる。

階乗を計算し、次数順に並べると x⁵/120 - x³/6 + x

T₅ = 16*x⁵ - 20*x³ + 5*x, T₃ = 4*x³ - 3*x, T₁ = x

T₅ の係数は 1/1920

x⁵/120 - x³/6 + x - T₅/1920 = -5*x³/32 + 383*x/384

T₃ の係数は -5/128

-5*x³/32 + 383*x/384 - (-5)*T₃/128 = 169*x/192

T₁ の係数は 169/192

T₅ を削除して

-5*T₃/128 + 169*T₁/192

= -5*x³/32 + 15*x/128 + 169*x/192 = -5*x³/32 + 383*x/384 が、テレスコーピングで得られた多項式となる。

次に、sin 関数を、11 次の項までふつうに近似したものと、13 次の近似式からテレスコーピングをおこない 11 次にしたものを比べてみる。

処理系の sin 関数が返す値（PCを使用しているので、Intel の数値演算プロセッサの sin 命令が返す値で、ほぼ真値と思われる）との誤差を、0～1 の範囲で比較するとこうなっている。

　　グラフ 1

上のグラフからわかるように、ふつうに近似した場合、端での誤差が大きくなっている。

ここで、0～0.6の範囲を拡大してみると以下のようになっている。

　　グラフ 2

テレスコーピングにより、誤差を（0～1の間に）うまく分散させた関数がえられたことがわかる。

次に、13 次まで近似したものと、それをテレスコーピングにより 11 次で近似したものを比べてみる。

　　グラフ 3

テレスコーピングがどのような操作であるのか、直感的にわかると思う。

関数の定義域

To be written.

プログラム

以上の計算のために作成したプログラムを示す。

telescope.scm

簡単な多項式処理系がコードの大半である

; telescope.scm

(use srfi-1)

; polynomial library

(define (make-polynomial factors)
  (let ((len (length factors)))
    (let ((degree (- len 1)))
      (cons degree factors) )))

(define (polynomial-degree polynomial) (car polynomial))

(define (polynomial-factors polynomial) (cdr polynomial))

(define (negate-polynomial polynomial)
  (let ((degree (polynomial-degree polynomial)))
    (cons degree (map - (polynomial-factors polynomial))) ))

(define (delete-zero-term polynomial)
  (let ((d (polynomial-degree polynomial))
        (fs (polynomial-factors polynomial)) )
    (if (zero? d)
        polynomial
        (if (zero? (car fs))
            (delete-zero-term (cons (- d 1) (cdr fs)))
            polynomial ))))

(define (add-polynomial poly0 poly1)
  (let ((d0 (polynomial-degree poly0))
        (d1 (polynomial-degree poly1))
        (fs0 (polynomial-factors poly0))
        (fs1 (polynomial-factors poly1)) )
    (delete-zero-term
      (if (eq? d0 d1)
          (cons d0 (map + fs0 fs1))
          (if (< d0 d1)
              (add-polynomial poly1 poly0)
              (let ((diff (- d0 d1)))
                (let ((lead (take fs0 diff))
                      (rest (drop fs0 diff)) )
                  (cons d0 (append lead (map + rest fs1))) )))))))

(define (sub-polynomial poly0 poly1)
  (add-polynomial poly0 (negate-polynomial poly1)) )

(define (factor-polynomial factor poly)
    (cons (polynomial-degree poly) (map (lambda (x) (* factor x)) (polynomial-factors poly))) )

(define (divmod-polynomial poly0 poly1)
  (define (division fs)
    (if (null? fs)
        (make-polynomial '(0))
        (make-polynomial (reverse fs)) ))
  (define (sub-factors list0 list1 tmp)
    (if (null? list1)
        (delete-zero-term (make-polynomial (append (reverse tmp) list0)))
        (sub-factors (cdr list0) (cdr list1) (cons (- (car list0) (car list1)) tmp)) ))
  (define (divmod-polynomial-it poly0 poly1 fs)
    (let ((d0 (polynomial-degree poly0))
          (d1 (polynomial-degree poly1)) )
      (if (< d0 d1)
          (list (division fs) poly0)
          (if (zero? d1)
              (list (factor-polynomial (/ 1 (car (polynomial-factors poly1))) poly0) (make-polynomial '(0)))
              (let ((fs0 (polynomial-factors poly0))
                    (fs1 (polynomial-factors poly1)) )
                (let ((f (/ (car fs0) (car fs1))))
                  (divmod-polynomial-it
                    (sub-factors fs0 (map (lambda  (x) (* f x)) fs1) '())
                    poly1
                    (cons f fs) )))))))
  (divmod-polynomial-it poly0 poly1 '()) )

(define (lift-polynomial poly) ; XXX
  (let ((d (polynomial-degree poly)))
    (cons (+ d 1) (append (polynomial-factors poly) '(0))) ))

; chebyshev library

; Tn
(define (chebyshev n)
  (cond ((eq? n 0) (make-polynomial '(1)))
        ((eq? n 1) (make-polynomial '(1 0)))
        (else (sub-polynomial
                (factor-polynomial 2 (lift-polynomial (chebyshev (- n 1))))
                (chebyshev (- n 2)) ))))

(define (chebyshev-list n)
  (if (zero? n)
      (list (chebyshev 0))
      (cons (chebyshev n) (chebyshev-list (- n 1))) ))

; chebyshev-ex

(define (chebyshev-ex poly)
  (define (chebyshev-ex-loop poly chebyshevs tmp)
    (if (null? chebyshevs)
        (reverse tmp)
        (let ((c (car chebyshevs)))
          (let ((divmod (divmod-polynomial poly c)))
            (chebyshev-ex-loop
              (cadr divmod)
              (cdr chebyshevs)
              (cons (car (polynomial-factors (car divmod))) tmp) )))))
  (chebyshev-ex-loop poly (chebyshev-list (polynomial-degree poly)) '()) )

; chebyshev-mix

(define (chebyshev-mix factors)
  (define (chebyshev-mix-loop factors chebyshevs poly)
    (if (null? factors)
        poly
        (chebyshev-mix-loop
          (cdr factors)
          (cdr chebyshevs)
          (add-polynomial (factor-polynomial (car factors) (car chebyshevs)) poly) )))
  (chebyshev-mix-loop factors (chebyshev-list (- (length factors) 1)) (make-polynomial '(0))) )

; factrial

(define (fact n)
  (define (fact-loop n a)
    (if (zero? n)
        a
        (fact-loop (- n 1) (* n a)) ))
  (fact-loop n 1) )

; trigon-ex 

(define (trigon-ex n)
  (define (trigon-ex-loop n tmp)
    (cond ((zero? n) (reverse (cons 1 tmp)))
          ((eq? n 1) (reverse (cons 0 (cons 1 tmp))))
          (else
            (let ((sign (if (< (modulo n 4) 2) 1 -1)))
              (trigon-ex-loop (- n 2) (cons 0 (cons (* sign (/ 1 (fact n))) tmp))) ))))
  (trigon-ex-loop n '()) )

telescope-test.scm

ユニットテスト

; telescope-test.scm

(use gauche.test)

(add-load-path ".")
(test-start "telescope")
(load "telescope")

(test-section "polynomial")
(define poly (make-polynomial '(3 2 1)))
(test* "make-polynomial" '(2 3 2 1) poly)
(test* "negate-polynomial" '(2 -3 -2 -1) (negate-polynomial poly))
(test* "delete-zero-term 1" '(0 0) (delete-zero-term '(0 0)))
(test* "delete-zero-term 2" '(2 3 2 1) (delete-zero-term '(3 0 3 2 1)))
(test* "add-polynomial 1" '(2 6 4 2) (add-polynomial poly poly))
(define poly2 (make-polynomial '(5 4 3 2 1)))
(test* "add-polynomial 2" '(4 5 4 6 4 2) (add-polynomial poly poly2))
(test* "sub-polynomial 1" '(0 0) (sub-polynomial poly poly))
(test* "sub-polynomial 2" '(4 -5 -4 0 0 0) (sub-polynomial poly poly2))
(test* "sub-polynomial 3" '(4 5 4 0 0 0) (sub-polynomial poly2 poly))
(test* "factor-polynomial" '(2 3/2 1 1/2) (factor-polynomial 1/2 poly))
(test* "divmod-polynomial 1" '((0 0) (2 3 2 1)) (divmod-polynomial poly poly2))
(test* "divmod-polynomial 2" '((0 1) (0 0)) (divmod-polynomial poly poly))
(test* "divmod-polynomial 3" '((0 1) (0 0)) (divmod-polynomial poly2 poly2))
(test* "divmod-polynomial 4" '((2 5/3 2/9 8/27) (1 32/27 19/27)) (divmod-polynomial poly2 poly))
(define poly3 (make-polynomial '(3)))
(test* "divmod-polynomial 5" '((2 1 2/3 1/3) (0 0)) (divmod-polynomial poly poly3))
(test* "divmod-polynomial 6" '((4 5/3 4/3 1 2/3 1/3) (0 0)) (divmod-polynomial poly2 poly3))


(test-section "chebyshev")
(test* "chebyshev 0" (make-polynomial '(1)) (chebyshev 0))
(test* "chebyshev 1" (make-polynomial '(1 0)) (chebyshev 1))
(test* "chebyshev 2" (make-polynomial '(2 0 -1)) (chebyshev 2))
(test* "chebyshev 3" (make-polynomial '(4 0 -3 0)) (chebyshev 3))
(test* "chebyshev 6" (make-polynomial '(32 0 -48 0 18 0 -1)) (chebyshev 6))
(test* "chebyshev-list" '((3 4 0 -3 0) (2 2 0 -1) (1 1 0) (0 1)) (chebyshev-list 3))
(test* "chebyshev-ex" '(1/2 1 3/2) (chebyshev-ex (make-polynomial '(1 1 1))))
(test* "chebyshev-mix" (make-polynomial '(1 1 1)) (chebyshev-mix '(1/2 1 3/2)))

(test-section "misc")
(test* "fact" 24 (fact 4))
(test* "trigon-ex 0" '(1) (trigon-ex 0))
(test* "trigon-ex 1" '(1 0) (trigon-ex 1))
(test* "trigon-ex 2" '(-1/2 0 1) (trigon-ex 2))
(test* "trigon-ex 3" (list (/ -1 (fact 3)) 0 1 0) (trigon-ex 3))
(test* "trigon-ex 4" (list (/ 1 (fact 4)) 0 -1/2 0 1) (trigon-ex 4))
(test* "trigon-ex 10" (list (/ -1 (fact 10)) 0 (/ 1 (fact 8)) 0 (/ -1 (fact 6)) 0 (/ 1 (fact 4)) 0 -1/2 0 1) (trigon-ex 10))
(test* "trigon-ex 11" (list (/ -1 (fact 11)) 0 (/ 1 (fact 9)) 0 (/ -1 (fact 7)) 0 (/ 1 (fact 5)) 0 (/ -1 (fact 3)) 0 1 0) (trigon-ex 11))

telescope-main.scm

多項式を求め、それを利用してユーザー定義版三角関数を定義している。main では gnuplot で plot できるデータを出力する。

簡単な解説

(trigon-ex n) 三角関数の n 次まで展開した近似式を返す。n が偶数か奇数かに応じて cos か sin のそれを返す
(chebyshev-ex poly) n 次の多項式 poly を n 次以下のチェビシェフ多項式の１次結合で表した時の係数のリストを返す
(chebyshev-mix (cdr cheb-facts-cos)) 最高次を削除してから、チェビシェフ多項式の係数をもとにベキ級数に戻している
実験の結果、16 次・17 次までの近似からテレスコーピングした結果で、倍精度浮動小数点として十分な誤差になっているようであった
級数の合計を計算する際、数値計算の精度を上げるために、負になる項と正になる項をまとめて計算しているが、ホーナー法で計算したほうがよいのかもしれない

; telescope-main.scm

(load "./telescope.scm")

;(define poly-cos (make-polynomial (trigon-ex 14)))  ; cos
;(define poly-sin (make-polynomial (trigon-ex 15)))  ; sin
(define poly-cos (make-polynomial (trigon-ex 16)))  ; cos
(define poly-sin (make-polynomial (trigon-ex 17)))  ; sin

(define cheb-facts-cos (chebyshev-ex poly-cos))
(define cheb-facts-sin (chebyshev-ex poly-sin))

;(define factors-cos
;  (map exact->inexact (polynomial-factors (chebyshev-mix cheb-facts-cos))) )
;(define factors-sin
;  (map exact->inexact (polynomial-factors (chebyshev-mix cheb-facts-sin))) )
(define factors-cos
  (map exact->inexact (polynomial-factors (chebyshev-mix (cdr cheb-facts-cos)))) )
(define factors-sin
  (map exact->inexact (polynomial-factors (chebyshev-mix (cdr cheb-facts-sin)))) )

(define (list-powers x n)
  (define (list-powers-loop n tmp)
    (if (zero? n)
        tmp
        (list-powers-loop (- n 1)  (cons (* x (car tmp)) tmp)) ))
  (list-powers-loop n '(1.0)) )

(define (sum xs)
  (define (sum-loop xs)
    (cond ((null? xs) '(0))
          ((null? (cdr xs)) xs)
          (else (sum-loop (sort (cons (+ (car xs) (cadr xs)) (cddr xs))))) ))
  (let ((negs (sort (map - (filter negative? xs))))
        (poss (sort (filter positive? xs))) )
    (- (car (sum-loop poss)) (car (sum-loop negs))) ))

(define (my-cos x)
  (let ((powers (list-powers x (- (length factors-cos) 1))))
;    (sum (cons (- (cos x)) (map * factors-cos powers))) ))
    (sum (map * factors-cos powers)) ))

(define (my-sin x)
  (let ((powers (list-powers x (- (length factors-sin) 1))))
;    (sum (cons (- (sin x)) (map * factors-sin powers))) ))
    (sum (map * factors-sin powers)) ))

(define (main args)
  (define (main-loop i n)
    (if (> i n)
        (exit)
        (let ((x (exact->inexact (/ i n))))
          (display x) (display " ") (display (my-cos x)) (newline)
          (display x) (display " ") (display (my-sin x)) (newline)
          (main-loop (+ i 1) n) )))
  (main-loop 0 1000) )