9. 再帰のバリエーション (生成的再帰) プログラミング論 I 本日の内容 生成的な再帰による繰り返し計算 1. 最大公約数：ユークリッドの互除法 2. クイックソート – 分割統治法の考え方 3. 中間値定理による方程式の求解 4. ニュートン法による非線形方程式の求解 生成的再帰 一般に複雑な問題を解きたい場合には… → 解決可能になるまで小さな問題に分割 • (今までは)データ定義に従った分割：構造的再帰(structural) • 例：リストを先頭と残りに分けて，残りを再帰処理 • 異なる観点による問題の分割： 生成的再帰(generative) • そもそも再帰的データを処理しない場合 • 特別な解き方を使用する場合 – 単純な構造的再帰より効率向上を目指す場合，など – 生成的再帰は ad hoc(設計でなくアルゴリズムの発明?) • 単なる数の列挙，数学的定理を利用，など様々 • 複雑なアルゴリズムの発明は専門家任せ? しかし – 理解や利用は必要 (単純なものなら自分で開発) 例題1. 最大公約数の計算 2つの自然数m, n から最大公約数を求める関数． 構造的再帰と生成的再帰で作成し比較 入力：2つの自然数 例：180 と 32 gcd 出力：1つの正整数 4 設計レシピに従い，規約，目的記述文，ヘッダを作成． 規約は正確な入力を指定： (0でなく) 1 以上の自然数 ;; gcd : N[>= 1] N[>= 1] -> N[>= 1] ;; to find the greatest common divisior of n and m (define (gcd n m) ...) 関数本体の具体化は? 構造的 vs. 生成的再帰 • 構造的再帰：自然数の再帰定義を利用しても解ける • 生成的再帰：しかしユークリッドの互除法が有用 最大公約数：構造的再帰 最大公約数：自然数 n, mの小さい方かそれ以下の正整数 [自然数の再帰的定義を利用した構造的再帰] 単純な発想(小さい方の数から始め，順に調べていく) iに，小さい方の自然数を代入して開始．n と m の両方を割 り切る最大の i を探す．割切れない場合 i を1だけ小さくし， 割切れるiが見つかるまで，以下の手順を繰返す． – i = 1 の場合: 最大公約数は 1 – i ≠ 1 の場合: • n を i で割った余りと m を i で割った余りが共に0? – そうなら最大公約数は i – そうでなければ i を 1 だけ小さくして繰り返す 最大公約数：構造的再帰 ;; gcd-structural : N[>= 1] N[>= 1] -> N ;; to find the greatest common divisior of n and m ;; structural recursion using data definition of N[>= 1] (define (gcd-structural n m) (local ((define (first-divisior i) (cond [(= i 1) 1] [else (cond [(and (= (remainder n i) 0) (= (remainder m i) 0)) i] [else (first-divisior (- i 1))])]))) (first-divisior (min m n)))) ;; 注： min は最小の数を返す組み込み関数 実は非効率的．小さい自然数なら良いが… • 例： (gcd-structural 1888819 528921) 結果は17 – DrRacketの (time 式) で時間測定→結構時間がかかっている → 「ユークリッドの互除法」による生成的再帰で時間短縮 gcd-structural without local ;;gcd-structural-wo-local:N N ->N (define (first-divisor n m i) (cond [(= i 1) i] [(and (= (remainder n i) 0) (= (remainder m i) 0)) i] [else (first-divisor n m (- i 1))])) ;; (define (gcd-structural n m) (first-divisor n m (min n m))) 生成的再帰：ユークリッドの互除法 2つの自然数 larger と smaller の最大公約数は，smaller と larger を smaller で割った剰余rとの最大公約数gに等しい larger = g*a, smaller=g*b，r=larger – smaller*c = g(a-bc) larger – smaller => smaller – (larger/smallerの余り) (gcd larger smaller) = (gcd smaller (remainder larger smaller)) 一種の再帰．0 はすべての数で割り切れるので smaller が 0 なら答は larger (再帰が停止する時)．つまり – smaller = 0 なら (基底(再帰しない)ケース) 最大公約数は larger – smaller ≠ 0 なら (再帰するケース) 最大公約数は， 「larger を smaller で割った余 り」 と smaller の最大公約数に等しい 最大公約数：生成的再帰 ;; gcd-generative : N[>= 1] N[>=1] -> N[>=1] ;; to find the greatest common divisior of m and n ;; generative recursion: ;; (gcd m n) = (gcd n (remainder m n)) where (>= m n) (define (gcd-generative n m) (local ((define (clever-gcd larger smaller) (cond [(= smaller 0) larger] [else (clever-gcd smaller (remainder larger smaller))]))) (clever-gcd (max m n) (min m n)))) ;; 注： max は最大、min は最小の数を返す組み込み関数 生成的再帰：ユークリッドの互除法 簡単のため，前ページの clever-gcdに焦点をあて説明 (前ページプログラムの larger を m，smaller を n とする) ;; clever-gcd: N N -> N ;; to find the greatest common divisor of n and m ;; Example: (clever-gcd 180 32) = 4 (define (clever-gcd m n) (cond 終了条件 [(= n 0) m] 解(再帰しない基底ケース) [else (clever-gcd n (remainder m n))])) (ユークリッド互除法に固有な)再帰 • clever-gcd の内部に clever-gcd の呼び出し – 「再帰」で clever-gcd の実行が繰り返される 例： (clever-gcd 180 32) = (clever-gcd 32 20) ⇒ 「データ構造(自然数)」でなく，問題固有の「特別なアル ゴリズム」による生成的再帰 (clever-gcd 180 32) から 4 を得る過程の概略 (clever-gcd 180 32) 最初の式 ① = … = (clever-gcd 32 20) ② = … = (clever-gcd 20 12) ③ = … = (clever-gcd 12 8) ④ = … = (clever-gcd 8 4) ⑤ = … = (clever-gcd 4 0) ⑥ 途中の計算過程 = … = 4 実行結果 m=180, n=32 のとき， clever-gcd は 6回の呼出。 練習1. 公約数の計算コスト 例題1 の最大公約数について 1888819 と 528921 の最大公約数 17 を求める 場合，構造的再帰と生成的再帰のプログラム で，それぞれの再帰の回数は何回か？ 参考：実際の実行時間を DrScheme の time 関 数を使って測定してみる．使用法： (time 式) 例題2. クイックソート • Hoare のクイックソートのアルゴリズムを使ったプ ログラム quick-sort を作る．(分割統治法) 入力はソート対 象となるリスト 例：(list 8 10 6 3 5) quick-sort 出力はソー ト済のリスト (list 3 5 6 8 10) 規約 contract (入力、出力)は第6回例題5のインサーション ソートと同じだが，その名の通り高速なソートを目指すもの。 (インサーションソートは再帰の回数がリストの長さ n に対 し n2 に比例するものだったが，クイックソートは平均で n log2n ) 分割統治法(divide and conquer) 基本方針：小さな問題のほうが解きやすい 問題を，複数の部分問題に分割して解き，部分問 題の解を統合して元の問題の解とする手法． ※ 各部分問題は元の問題と同類の問題であるが，入力デー タの大きさが小さい問題． – インサーションソートでは(データ構造に沿って)先頭要 素と残りの部分に分割 – クイックソートでは分割を工夫 • 「 pivot (分割基準)より大きい要素のリスト」に分割しソート • 「 pivot より小さい要素のリスト」に分割しソート クイックソートの考え方 リスト 8 10 6 3 5 先頭を pivot に 8 pivotより小 635 pivot より大 empty ス(リストが空)に なるまで pivot の 選択と pivot によ るリストの分割を 繰返す(再帰) empty 統合： (再帰呼出 10 pivot pivot 6 35 10 empty empty 分割：基底ケー 中略 35 6 356 empty empty 8 3 5 6 8 10 10 10 による) 中間結果 と pivot のリストを 順序を保存し結合 クイックソートの処理手順 1. 基底ケースかどうかを調べる – 基底ケース(empty)の場合， empty を返して終了 – そうでなければ 手順 2 へ 2. 基準値(pivot)より大きな要素のリストと小さな 要素のリストを生成し，それらを再帰的にソート (より小さな部分問題の生成) – pivotはソート対象の中から1つ選ぶ(今回は先頭) 3. pivotより大きな(および小さな)要素のリストに 対するソート結果(部分問題の解)と pivot を，順 序関係を保存して， 1つのリストへと結合 例題2. クイックソート • 分割統治法によるアルゴリズムを使った生成的再帰 1. 部分問題への分割が必要か調べ，必要なら分割． 2. 第7回例題2の関数 filter1 を活用し， pivotより大きな要素 のリストと小さなリストに分割し，それらに再帰的に実行 3. 分割の pivot は，リストの先頭要素を使う． 4. 分割した2つの小問題の結果リストと pivot (のリスト)を結合する には組込み関数の append を使う． 補足：append は Scheme が備えているリストを結合する関数 (append リストの並び) 例： (append (list 1 2) (list 3 4)) = (list 1 2 3 4) (append (list 1 2) (list 3 4) (list 5)) = (list 1 2 3 4 5 6) (append (list 1 2) 3 (list 4 5)) 単純値の3はリストでないのでエラー クイックソートのプログラム ;; quick-sort : (listof number) -> (listof number) ;; to create a list of numbers with the same numbers as ;; alon sorted in ascending order 小問題に分 (define (quick-sort alon) (cond 終了条件 pivotはリスト 割して再帰 [(empty? alon) empty] の先頭要素 基底(再帰しない)ケース [else (append (quick-sort (filter1 < alon (first alon))) append で結合 (list (first alon)) (quick-sort (filter1 > alon (first alon))))])) • quick-sort の内部で quick-sort を呼出して繰り返し ⇒ まさに「再帰」．ただし小問題への分割はデータ構造 でなく固有のアルゴリズムによる「生成的再帰」 (quick-sort (list 6 2 4)) からの過程 (quick-sort (list 6 2 4)) = … =(append (quick-sort (filter1 < (list 6 2 4) (first (list 6 2 4)))) (list (first (list 6 2 4))) (quick-sort (filter1 > (list 6 2 4) (first (list 6 2 4)))) ) 3つのリストを結合したものが答 (quick-sort (filter1 < (list 6 2 4) 6) → 6 より小さい部分 (list 6) → (list 6) (quick-sort (filter1 > (list 6 2 4) 6) → 6 より大きい部分 に分割して問題を解き進めている 練習2. クイックソート 例題2のクイックソートプログラムについて 1. 例題のプログラムは要素を昇順に並べるが，これ を降順になるように作り変えよ． 2. 例題のままのプログラムだと ，pivot と同じ要素が リスト中に複数あってもソート結果にはひとつしか 残らない(要素数が減ってしまう)．これを改善し， 要素数を変えないプログラムを作成せよ． 3. リストが空の場合を基底ケースとしたが，リストの 要素数が1の場合でも基底ケースになりえる． これを反映したプログラムを作成せよ． 例題3. 二分探索法(区間二分法) 区間二分法：分割統治法の一種(分割後の一部に着目) 中間値定理による生成的再帰で，数値関数 f のルート値 (f(r) = 0 を満たす数) r を求める関数find-rootを作成 f(x) 中間値定理：連続関数 f は f(a) と f(b) の符号が異なれ ば閉区間 [a,b] で根を持つ 入力は1つの関 数と2つの数値 (区間の両端) f(b) a 0 f(a) find-root b x ルート値 (f (a) < 0 < f(b) の例) 出力は1つ の数値 1.4142131805419921875 例：f(x)=x2-2, 0, 2 f(x) = x2 - 2 の解√２と－√２のうち区間[0,2]にある√２に相当 区間二分法の手順 • 関数 f とルートの存在を期待する区間の両端：left, right • アルゴリズムが機能するための前提： 関数が left と right に対し異なる符号を持つ •区間を二分しルートが存在する側の選択を繰返す (以下は f (left) < 0 < f(right) の場合) 1. 区間 [left, right] を2等分。中点は m = (left+right) / 2 2. f(m) の値を計算し 0 と比較： f(right) f(m) 1. f(m) > 0 → 解は[left, m]にある 2. f(m) < 0 → 解は[m, right]にある left f(left) m right 3. ルートが存在する区間を選択し上記を繰り返す．「区 間」の幅を小さくし，ルートが存在する非常に小さな範 囲を決定→ 近似解 区間二分法の処理手順 1. 区間の両端をleft, right, m を (left+right)/2とする • f(left)が0なら終了．結果はleft (m,rightも同様に) • 区間が非常に小さければ終了．結果は m (収束) • そうでなければ以下を実行 2. 区間を半分にして以下の判断で区間を選択 f(left) と f(m) が互いに異符号： ルートは左側と期待→左の区間を選択し再帰 f(m) と f(right) が互いに異符号： ルートは右側と期待→右の区間を選択し再帰 区間二分法の処理手順 f(x) = x2-2, left = 0, right = 2 の場合 left = 0, right = 2 このとき、m=1, f(left)=-2, f(m)= -1, f(right) = 2 f(m) < 0 < f(2) なので ⇒ m を次の再帰の left に left = 1, right = 2 このとき、m=1.5, f(left)=-1, f(m)= 0.25, f(right) = 2 f(left) < 0 < f(m) なので ⇒ m を次の再帰の right に left = 1, right = 1.5 m=1.25, f(left)=-1, f(m)= -0.4375, f(right) = 0.25 f(m) < 0 < f(right) なので ⇒ m を次の再帰の left に 以下略 区間二分法 find-root 中間値mをlocal式で変数定義し何度も同一計算するのを回避 ;;find-root : (number -> number) number number -> number ;;determine R such that f has a root between R±TOLERANCE/2 ;;ASSUMPTION:f is continuous and monotonic, and satisfies ;; (or (<= (f left) 0 (f right)) (<= (f right) 0 (f left))) (define (find-root f left right) (local ((define m (/ (+ left right) 2))) (cond [(= (f right) 0) right] 基底：たまたまルートになった [(= (f left) 0) left] [(= (f m) 0) m] それぞ [(<= (- right left) TOLERANCE) m ] 基底：区間幅が十分小 れ、右、[else 左の区 (cond 間を選 [(<= (* (f m) (f right)) 0) (find-root f m right)] 択して [(<= (* (f left) (f m)) 0) (find-root f left m)] [ else (list "can't find root" left m right)])]))) 再帰 ;； 前提を満たさない区 (define TOLERANCE 0.000001) 間では計算できない 例題4. ニュートン法 ニュートン法による非線形方程式 f(x) = 0 の求解 newton 入力：関数 f と， 解の推測値 guess 出力：f(x) = 0 の 近似解の1つ newton は，関数を入力とする関数(高階関数) – 近似解 x0, x1, x2 ... の収束の様子を理解 – 繰返しによる近似計算で解が求まらない（停止しな い）場合がありえることを理解 ニュートン法 • 初期近似値 x0 から出発 • 反復公式 f(xi)=f’(xi)(xi – xi+1) xi 1 f ( xi ) xi f ' ( xi ) を繰り返す。これは，点（xi, f(xi) ）における 接線 と，x軸 (y=0) との交点 (xi+1, 0) を求めている。 • x1, x2, x3 ... は徐々に解に近づくと期待できる。 8 x3 – 6x2+11x –6 6 4 解 2 解 解 x 0 0 -2 -4 -6 -8 0.5 1 1.5 2 2.5 3 3.5 4 4.5 8 6 4 2 x 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2 -4 -6 -8 関数上の点 関数上の点を1つ選ぶ 8 6 接線 4 2 x 0 0 0.5 1 1.5 2 2.5 3 -2 -4 -6 -8 関数上の点 接線を引く 3.5 4 4.5 8 6 接線 4 接線とx軸 2 の交点 x 解の1つ 0 0 0.5 1 1.5 2 2.5 3 3.5 4 -2 -4 -6 -8 関数上の点 接線と x 軸の交点は解 の1つに近づいたと期待 4.5 8 接線と x 軸の交点 6 (x0 - f(x0)/f '(x0), 0) 4 接線 y = f '(x0)(x-x0) + f(x0) 2 x 0 -2 -4 -6 -8 0 1 2 3 4 5 関数上の点 xi , f ( xi ) から 関数上の点 (x0, f(x0)) 接線とx軸の交点のx座標 xi 1 xi f ( xi ) / f ' ( xi ) を求めることを繰返す ニュートン法の例 f(x) = x3 – 6x2+11x –6, x0 = 0 では： x0 = 0 f(x0) = -6, f '(x0) = 11 ⇒ x1 = x0 - f(x0)/f '(x0) = 0.54545454... x1 = 0.54545454... f(x1) = -1.62283996..., f '(x1) = 5.3471074... ⇒ x2 = x1 - f(x1)/f '(x1) = 0.84895321... x2 = 0.84895321... f(x2) = 0.37398512..., f '(x2) = 2.9747261... ⇒ x3 = x2 - f(x2)/f '(x2) = 0.97460407... 8 6 接線と x 軸の交点 (0.54545454..., 0) x14 = 0.54545454... 2 x 0 0 0.5 1 1.5 2 2.5 -2 -4 -6 -8 関数上の点 (0, -6) x0 = 0 3 3.5 4 4.5 8 6 接線と x 軸の交点 (0.84895321..., 0) x4 2 = 0.84895321... 2 x 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -2 -4 -6 -8 関数上の点 (0.54545454..., -1.62283996...) x1 = 0.54545454... 収束するまで xiの計算を繰り返す ニュートン法の繰り返し処理 1. 「収束」したなら⇒再帰終了：現在の xn が解 2. そうで無ければ： 次の近似値xn+1の値「 xnでの 接線と，x軸の交点のx 座標値」を使い再帰 ここで「収束」とは? ニュートン法では現在の xn がどれだけ 真の値に近いかは正確には分からない⇒今回は次の判定 ある小さな正数δに対し次が成り立つか? 成立した時点で再帰を終了し xn を解とする f ( xn ) f ( xn ) なら true，さもなくば false は次のように書く (abs は「絶対値」を求める関数，TOLERANCE はδ) (<= (abs (f xn)) TOLERANCE) 接線の傾きを求める ;; d/dx: (number->number) number number -> number ;; inclination of the tangent (x+h,f(x+h)) (define (d/dx f x h) h f(x+h)-f(x-h) (/ (- (f (+ x h)) (f (- x h))) (* 2 h))) (x-h,f(x-h)) x+2h 接線と x 軸の交点のx座標を求める ;; find-root-tangent: (number -> number) number -> number ;; to find the root of the tagent of f at r0 (define (find-root-tangent f r0) (local ((define H 0.00001)) (- r0 (/ (f r0) (d/dx f r0 H))))) ニュートン法による再帰関数 ;; newton : (number -> number) number -> number ;; to find a number guess such that ;; (<= (abs (f guess)) TOLERANCE) (define (newton f guess) (cond [(<= (abs (f guess)) TOLERANCE) guess] [else (newton f (find-root-tangent f guess))])) ;; (define TOLERANCE 0.0001) (define (f2 x) (+ (* x x x) (* -6 x x) (* 11 x) -6)) 補足. 接線の傾き • 接線の傾きを求める関数 d/dx – 関数 f と座標値 x, 近似幅 hから，x における f の傾 き (つまり f ' (x) ) の近似値を求める 接線 f(x) f ( x h) f ( x h) 0 xh xh 傾きは f ( x h) f ( x h) 2 h 0 に近い h の値で、 接線の傾きf '(x) の 「近似値」を求める x ;;d/dx:(number->number) number number -> number ;; inclination of the tangent (define (d/dx f x h) (/ (- (f (+ x h)) (f (- x h))) (* 2 h))) (newton f2 0) から結果が得られる過程 (newton f2 0) 最初の式 = ... = (newton f2 0.5454545449586776864012021...) = ... = (newton f2 0.8489532098092286727580919...) = ... = (newton f2 0.9746740703633078705978850...) = ... = (newton f2 0.9990915478975544584656251...) = ... = (newton f2 0.9999987646860979280445576...) コンピュータ内部での計算 = ... = 0.9999987646910051384499455... 実行結果 ニュートン法の能力と限界 • 虚数解は求まらない f(x) • 求まる解はあくまで「近似解」 – 近似解→どの程度の精度 で計算するか(やめるか)に より、収束時間、丸め誤差の 影響など変化 • 初期近似値の選び方によって – 求まる解が変わり得る – 収束に問題：(解と離れて いると)時間がかかる、収束 しない場合がある、など ２ ３ 0 １ x 関数ｆ(x)が単調でなく変曲点を持つ （つまりｆ’(x)の符号が変わる）とき 例）上の図の1から開始 2 (f'(x) = 0となる点)が選ばれる 3 ： y 軸との交点が求まらない （負の無限大に発散） → 収束しない 収束対策 繰り返しの上限回数 を設定→課題 課題 課題1. クイックソート 1. 例題2のクイックソートプログラムについて，pivot と比 較してリストを分割する際の filter1 の比較演算に < ではなく <= を使ってしまった場合(他の部分はそ のまま)，停止性に問題が生じる．(list 5) をソート する場合を例に，その理由を述べよ． 2. 次のような構造体 AddressNote のデータを要素とす るリストを考える．そのリスト要素のAddressNote データを名前 nameの昇順にソートするクイックソート のプログラムを作れ． (define-struct AddressNote (name age address)) ここで name と address は文字列，age は自然数 課題2. Newton 法 例題4 Newton 法プログラムを次のように変更せよ 1. 精度を変更できるよう， TOLERANCE を変数定義で はなく関数 newton への引数 t として与える． 2. 収束しない場合でも停止するように変更する． 最大の繰り返し回数を引数として指定でき，指定 回数に達した場合は値が収束しなくても再帰を停 止するようなプログラムを作成せよ．
© Copyright 2024 ExpyDoc
