foldl,foldr

Hugs.Base> foldl (/) 2 [1, 2, 4, 8]
0.03125
Hugs.Base> foldr (/) 2 [1, 2, 4, 8]
0.5

foldl (/) 2 [1, 2, 4, 8] は、(((2 / 1) / 2) / 4) / 8
foldr (/) 2 [1, 2, 4, 8] は、1 / (2 / (4 / (8 / 2)))ということですね。

Hugs.Base> foldl (/) 3 [1, 3, 5, 7, 2, 1]
0.0142857142857143
Hugs.Base> foldr (/) 3 [1, 3, 5, 7, 2, 1]
1.42857142857143

Schemeで書いてみる。Schemeのライブラリ、SRFIのfoldは

> (fold / 2 '(1 2 4 8))
8

これは8 / (4 / (2 / (1 / 2)))になってしまう。自分で書くしかない。

(define (foldl proc elem list)
  (define (foldl-iter list ans)
    (cond ((null? list) ans)
          (#t (foldl-iter (cdr list) (proc ans (car list))))))
  (foldl-iter list elem))

(define (foldr proc elem list)
  (define (foldr-rec list)
    (cond ((null? list) elem)
          (#t (proc (car list) (foldr-rec (cdr list))))))
  (foldr-rec list))

> (foldl / 2.0 '(1 2 4 8))
0.03125
> (foldr / 2.0 '(1 2 4 8))
0.5

> (foldl / 3.0 '(1 3 5 7 2 1 ))
0.014285714285714287
> (foldr / 3.0 '(1 3 5 7 2 1 ))
1.4285714285714284

いい感じに頭がこんがらがってきました。

次は Pythonで。

def foldl(proc, elem, lis):
    lis.insert(0, elem)
    return reduce(proc, lis)

def foldr(proc, elem, lis):
    lis.append(elem)
    lis.reverse()
    return reduce(lambda x, y : apply(proc, [y, x]), lis)

fl = foldl(lambda x, y : x / y, 2.0, [1.0, 2.0, 4.0, 8.0])
fr = foldr(lambda x, y : x / y, 2.0, [1.0, 2.0, 4.0, 8.0])

>>> fl
0.03125
>>> fr
0.5
>>> print foldl(lambda x, y : x / y, 3.0, [1.0, 3.0, 5.0, 7.0, 2.0, 1.0])
0.0142857142857
>>> print foldr(lambda x, y : x / y, 3.0, [1.0, 3.0, 5.0, 7.0, 2.0, 1.0])
1.42857142857

Pythonて絶対こんな書き方をする言語ではないと思った。

ちょっと暇ができたので、scheme分を補充しておこう……。Durand-Kerner法(DK法)というのはニュートン法の改良版で、代数方程式の全ての解を一度に求めることができる。

数式の表示にはMathMLを使っています。また僕は数値計算にあまり詳しくないため、間違いなどありましたら指摘してください。

Durand-Kerner法で多項式の近似解を求めるプログラム

(require (lib "1.ss" "srfi"));SRFI-1のロード

(define PI 3.1415)
(define (square x) (* x x))
  
(define (dka input)     ;Durand-Kerner法
  (define formula (if (proper-list? input) 
                      (drop-while zero? input) ;最高次係数が0なら削除
                      (error "INVALID INPUT")))
  (define degree (- (length formula) 1))
  (define degree_list (iota degree 1 1))　
                        ;次数nに対して、'(1 2 ... n)というリストを作成
  (define zcenter 
    (if (>= 0 degree) (error "INVALID DEGREE")
        (- (/ (cadr formula) (* (car formula) degree)))))
  
  (define (get-R list)  ;解を囲む円の半径の候補を求める
    (let ((denom (car list)))
      (map (lambda (x k) 
             (expt (/ (abs x) denom) (/ 1 k))) 
           (cdr list) degree_list)))

  (define trial-list　　　　;Aberthの初期値
    (let* ((angle_base 
            (/ (* 2 PI (make-rectangular 0 1)) degree))
           (R (* degree 
                 (fold max 0　(get-R (kumitate-exp formula zcenter))))))
      (map (lambda (x) 
             (+ zcenter (* R (exp (* angle_base (- x 0.75)))))) 
           degree_list)))

  (define (dka-iter list)　　;反復
    (let ((deriv-list　
           (map (lambda (x)
                  (subst　(get-form (kumitate formula x)) x))
                list)))　　　　;それぞれの近似値における導関数値を求める
      (cond ((< (abs-list 
                 (map (lambda (x) (subst formula x)) list))
                0.0001) list);それぞれの近似値を与式に代入し、誤差を調べる
            (#t (dka-iter 
                 (map (lambda (x y) 
                        (if (= y 0) x (- x (/ (subst formula x) y)))) 
                      list deriv-list))))))
  (cond ((= degree 1) zcenter)
        (#t (dka-iter trial-list))))
             
(define (kumitate-exp list x) ;組み立て除法を用いて式を変形
  (define (kumitate-t-iter list ans)
  (cond ((null? list) ans)
        (#t (let ((form (kumitate list x)))
              (kumitate-t-iter (get-form form) (cons (last form) ans))))))
  (kumitate-t-iter list '()))

(define (kumitate list x);組み立て除法
  (define (kumitate-iter list v)
    (cond ((null? list) '())
          (#t (cons (+ (car list) (* v x)) 
                    (kumitate-iter (cdr list) (+ (car list) (* v x)) )))))
  (cons (car list) (kumitate-iter (cdr list) (car list))))

(define (get-form list)	;組み立て除法後のリストlistから、式の部分を取得
  (drop-right list 1))

(define (abs-list list)
  (fold + 0 (map (lambda (x) (+ (magnitude x))) list)))

(define (subst form v)	;式formに値vを代入した結果を返す
  (last (kumitate form v)))

リスト操作ライブラリ

SRFI-1を使う。proper-list?、 last、drop-while、drop-right、fold、iotaあたりがそう。便利すぎ。

式の表現

多項式は最高次係数を先頭とするリストであらわすことにした。つまり$x^2-5x+6$は'(1 -5 6)。

式にたいして値が代入できなきゃどうしようもないので、それを行う関数substを定義。計算には組み立て除法を用いる。

> (subst '(1 -5 6) 2)
0
> (subst '(1 -5 6) 4)
2

組み立て除法

組み立て除法というのはある多項式を (x-a) といった一次式で割るのに適した方法。それを使って商関数Q(x)と剰余Rを求めると、 $f(x)=(x-a)Q(x)+R$ より、 $f(a)=R$

> (kumitate '(1 -4 -2 34 2 4 3) 2)
(1 -2 -6 22 46 96 195)
> (kumitate '(1 -2 -6 22 46 96) 2)
(1 0 -6 10 66 228)

組み立て除法によって導関数値を求める

組み立て除法を繰り返すと、 $x^6-4x^5-2x^4+34x^3+2x^2+4x+3=(x-2)(x^5-2x^4-6x^3+22x^2+46x+96)+195={(x-2)}^2(x^4-6x^2+10x+66)+(x-2)228+195$ ...と更なる変形が得られる。これを使うとaに対するm次の導関数値が比較的簡単に求められる。

いま、f(x)が $$f(x)={(x-a)}^{m}Q(x)+{(x-a)}^{(m-1)}{R}_m...$$ と変形できたとして、点 x = a における m 次導関数値を求めたい。m 回微分すると次数がm-1以下の項はすべて消えることから、 $$f^{(m)}(x)={({(x-a)}^{m}Q(x))}^{(m)}$$ としてよい。

$${(g(x)h(x))}^{(m)}=\sum _{k=1}^m\left(\genfrac{}{}{0ex}{}{m}{k}\right)g^{(m-k)}(x)h^{(m)}(x)$$ を考える。これを展開したとき、( x-a )が含まれる項はxにaを代入したとき全て0になるので、結局 $$f^{(m)}(a)=m!Q_m(a)$$ が言える。まーここでは一次導関数値しか使わないけど。

Aberthの初期値

ある代数方程式の解が含まれる領域を推定する方法 。直感的には、解の分布の中心zcenterから、その円内にすべての解が含まれるように半径Rの円を描いて、その円周上に初期値を配置する。こうすることで近似値それぞれが真値に近づくスピードをそろえることができる。

解の中心は簡単に求められる。なぜなら、$a(x-b)(x-c)(x-d)=a{x}^3-a(b+c+d)x^2+...$

半径を求めるには、与式を解の中心zcenterで組み立て除法を用いて変形し、r = ( z - zcenter ) と置換することで円の半径 r が解となるような方程式を立てればよい。

結局、n個の初期値を円周上に等間隔に配置するには $$z_j=-\frac{a_1}{n{a}_0}+Rexp(\frac{2\pi i}{n}(j-\frac{3}{4}))$$ を j = 1～n について計算すればよい。

こうして求めた初期値はリストとしてtrial-listに格納される。

  (define (get-R list)  ;解を囲む円の半径の候補を求める
    (let ((denom (car list)))
      (map (lambda (x k) 
             (expt (/ (abs x) denom) (/ 1 k))) 
           (cdr list) degree_list)))

  (define trial-list　　　　;Aberthの初期値
    (let* ((angle_base 
            (/ (* 2 PI (make-rectangular 0 1)) degree))
           (R (* degree 
                 (fold max 0　(get-R (kumitate-exp formula zcenter))))))
      (map (lambda (x) 
             (+ zcenter (* R (exp (* angle_base (- x 0.75)))))) 
           degree_list)))

Durand-Kerner法

上の方法で求めた初期値それぞれについて、ニュートン法による近似を行う。本来は導関数値をさらに近似する。

$P_N(x)=a_0(x-a_1)(x-a_2)...(x-a_n)$とすると、その一次導関数は $$P_N^{(1)}(x)=a_0(x-a_2)\cdots (x-a_n)+a_0(x-a_1)(x-a_3)\cdots (x-a_n)+a_0(x-a_1)(x-a_2)(x-a_4)\cdots (x-a_n)+\cdots =a_0\sum _{i=1}^n\prod _{j=1,j\ne i}^n(x-a_j)$$ ここで、$z_i\approx a_i$ とすると、 $a_i$ の含まれる項は全て消えるので、 $$P_N^{(1)}(z_i)==a_0\prod _{j=1,j\ne i}^n(z_i-a_j)$$ このとき、$z_j\approx a_j$もまた言えるので、 $$P_N^{(1)}(z_i)==a_0\prod _{j=1,j\ne i}^n(z_i-z_j)$$

反復処理

  (define (dka-iter list)　　;反復
    (let ((deriv-list　
           (map (lambda (x)
                  (subst　(get-form (kumitate formula x)) x))
                list)))　　　　;それぞれの近似値における導関数値を求める
      (cond ((< (abs-list 
                 (map (lambda (x) (subst formula x)) list))
                0.0001) list);それぞれの近似値を与式に代入し、誤差を調べる
            (#t (dka-iter 
                 (map (lambda (x y) 
                        (if (= y 0) x (- x (/ (subst formula x) y)))) 
                      list deriv-list))))))

まず組み立て除法を用いて各近似値に対応する導関数値をそれぞれ計算し、deriv-listとする。それぞれの近似値を与式に代入してみて、全てが0に近ければ十分収束しているとして解を出力、終了。まだダメそうならmapを用いて、近似値全てをニュートン法で再近似。

導関数値リストはもともと別の場所でも使うはずだったので一次変数にしてあるが、その部分を消したためdka-iterの反復呼び出しの際に直接計算してもいい。でもそうするとすごいわかりにくくなると思うのでそのまま。上で示したように、点aにおける一階導関数値を求めるにはまずkumitateを用いて一回変形し、変形後のリストからget-formを用いて式の部分のみを抜き出し、それにaを代入すればよい。

結果

$$x^2-2x+1=0$$ $$x^3-3x^2-10x+24=0$$ $$x^4-5x^3-21x^2-22x-16=0$$ $$8x^5+24x^4-25x^3+3x^2+98=0$$ について、

> (dka '(1 -2 1))
(1 1)

> (dka '(1 -3 -10 24))
(3.9999998318967864+2.8387151406426836e-008i
 -3.000000000000001+1.441359404959568e-017i
 2.0000000001780975+1.4040173748418967e-010i)

> (dka '(1 -5 -21 -22 -16))
(8.0+0.0i
 -0.4999999999911601+0.8660254037577851i
 -2.0000000003879568+5.235939426694963e-010i
 -0.500000298422228-0.8660254066797882i)

> (dka '(8 24 -25 3 0 98))
(1.151630579359255+0.817203444821391i
 -0.7099165021391798+1.038229298587666i
 -3.883428154439853-9.183549615799121e-041i
 -0.7099165035189984-1.0382292988192343i
 1.151630579277084-0.8172034448086737i)

それぞれ誤差は、

0
2.3890323896364474e-006
7.669790640016116e-006
4.0979890195062324e-007

もっとシンプルにできそうな予感はするんだけど思いつかないなー。