9. 再帰のバリエーション (生成的再帰) プログラミング論 I

本日の内容 生成的な再帰による繰り返し計算 最大公約数：ユークリッドの互除法 クイックソート 中間値定理による方程式の求解 分割統治法の考え方 中間値定理による方程式の求解 ニュートン法による非線形方程式の求解

生成的再帰 一般に複雑な問題を解きたい場合には… → 解決可能になるまで小さな問題に分割 (今までは)データ定義に従った分割：構造的再帰(structural) 例：リストを先頭と残りに分けて，残りを再帰処理 異なる観点による問題の分割： 生成的再帰(generative) そもそも再帰的データを処理しない場合 特別な解き方を使用する場合 単純な構造的再帰より効率向上を目指す場合，など 生成的再帰は ad hoc(設計でなくアルゴリズムの発明?) 単なる数の列挙，数学的定理を利用，など様々 複雑なアルゴリズムの発明は専門家任せ? しかし 理解や利用は必要 (単純なものなら自分で開発)

例題1. 最大公約数の計算 2つの自然数m, n から最大公約数を求める関数． 構造的再帰と生成的再帰で作成し比較 gcd 入力：2つの自然数 出力：1つの正整数 gcd 例：180 と 32 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)))

生成的再帰：ユークリッドの互除法 smaller = 0 なら (基底(再帰しない)ケース) 最大公約数は larger 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 の最大公約数に等しい

最大公約数：生成的再帰 ;; 注： max は最大、min は最小の数を返す組み込み関数 ;; 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 180 32) = (clever-gcd 32 20) 簡単のため，前ページの 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 を作る．(分割統治法) 入力はソート対 象となるリスト 出力はソー ト済のリスト quick-sort 例：(list 8 10 6 3 5) (list 3 5 6 8 10) 規約 contract (入力、出力)は第6回例題5のインサーションソートと同じだが，その名の通り高速なソートを目指すもの。(インサーションソートは再帰の回数がリストの長さ n に対し n2 に比例するものだったが，クイックソートは平均で n log2n )

分割統治法(divide and conquer) 基本方針：小さな問題のほうが解きやすい 問題を，複数の部分問題に分割して解き，部分問題の解を統合して元の問題の解とする手法．　※各部分問題は元の問題と同類の問題であるが，入力データの大きさが小さい問題． インサーションソートでは(データ構造に沿って)先頭要素と残りの部分に分割 クイックソートでは分割を工夫 「 pivot (分割基準)より大きい要素のリスト」に分割しソート 「 pivot より小さい要素のリスト」に分割しソート

クイックソートの考え方 分割：基底ケー 統合： (再帰呼出 リスト 先頭を pivot に ス(リストが空)に なるまで pivot の 8 10 6 3 5 分割：基底ケー ス(リストが空)に なるまで pivot の 選択と pivot によ るリストの分割を 繰返す(再帰) 先頭を pivot に pivot より大 8 pivotより小 10 6 3 5 pivot pivot 6 10 3 5 empty empty empty 中略 統合： (再帰呼出 による) 中間結果 と pivot のリストを 順序を保存し結合 3 5 6 empty empty 10 empty 3 5 6 8 10 3 5 6 8 10

クイックソートの処理手順 基底ケースかどうかを調べる 基底ケース(empty)の場合， empty を返して終了 そうでなければ 手順 2 へ 基準値(pivot)より大きな要素のリストと小さな要素のリストを生成し，それらを再帰的にソート(より小さな部分問題の生成) pivotはソート対象の中から1つ選ぶ(今回は先頭) pivotより大きな(および小さな)要素のリストに対するソート結果(部分問題の解)と pivot を，順序関係を保存して， 1つのリストへと結合

例題2. クイックソート 分割統治法によるアルゴリズムを使った生成的再帰 部分問題への分割が必要か調べ，必要なら分割． 第7回例題2の関数 filter1 を活用し， pivotより大きな要素のリストと小さなリストに分割し，それらに再帰的に実行 分割の pivot は，リストの先頭要素を使う． 分割した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 の内部で quick-sort を呼出して繰り返し ;; 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 [(empty? alon) empty] [else (append (quick-sort (filter1 < alon (first alon))) (list (first alon)) (quick-sort (filter1 > alon (first alon))))])) 小問題に分 割して再帰 終了条件 pivotはリスト の先頭要素 基底(再帰しない)ケース append で結合 quick-sort の内部で quick-sort を呼出して繰り返し ⇒ まさに「再帰」．ただし小問題への分割はデータ構造でなく固有のアルゴリズムによる「生成的再帰」

(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のクイックソートプログラムについて 例題のプログラムは要素を昇順に並べるが，これを降順になるように作り変えよ． 例題のままのプログラムだと ，pivot と同じ要素がリスト中に複数あってもソート結果にはひとつしか残らない(要素数が減ってしまう)．これを改善し，要素数を変えないプログラムを作成せよ． リストが空の場合を基底ケースとしたが，リストの要素数が1の場合でも基底ケースになりえる． これを反映したプログラムを作成せよ．

例題3. 二分探索法(区間二分法) 区間二分法：分割統治法の一種(分割後の一部に着目) 中間値定理による生成的再帰で，数値関数 f のルート値 (f(r) = 0 を満たす数) r を求める関数find-rootを作成 f(x) f(b) 中間値定理：連続関数 f は f(a) と f(b) の符号が異なれ ば閉区間 [a,b] で根を持つ a x b f(a) ルート値 (f (a) < 0 < f(b) の例) f(x) = x2 - 2 の解√２と－√２のうち区間[0,2]にある√２に相当 find-root 1.4142131805419921875 入力は1つの関 数と2つの数値 (区間の両端) 出力は1つ の数値 例：f(x)=x2-2, 0, 2

区間二分法の手順 関数 f とルートの存在を期待する区間の両端：left, right アルゴリズムが機能するための前提： 区間を二分しルートが存在する側の選択を繰返す (以下は f (left) < 0 < f(right) の場合) 区間 [left, right] を2等分。中点は m = (left+right) / 2 f(m) の値を計算し 0 と比較： f(m) > 0 → 解は[left, m]にある f(m) < 0 → 解は[m, right]にある ルートが存在する区間を選択し上記を繰り返す．「区間」の幅を小さくし，ルートが存在する非常に小さな範囲を決定→ 近似解 f(right) left right f(left) f(m) m

区間二分法の処理手順 区間の両端を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 [(<= (* (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 近似解 x0, x1, x2 ...　の収束の様子を理解 繰返しによる近似計算で解が求まらない（停止しない）場合がありえることを理解 newton 出力：f(x) = 0 の 近似解の1つ 入力：関数 f と， 解の推測値 guess newton は，関数を入力とする関数(高階関数)

f(xi)=f’(xi)(xi – xi+1) ニュートン法 初期近似値 x0 から出発 反復公式 を繰り返す。これは，点（xi, f(xi) ）における 接線と，x軸 (y=0) との交点 (xi+1, 0) を求めている。 x1, x2, x3 ... は徐々に解に近づくと期待できる。 f(xi)=f’(xi)(xi – xi+1)

x3 – 6x2+11x –6 x 解 解 解

x 関数上の点 関数上の点を1つ選ぶ

接線 x 接線を引く 関数上の点

接線 接線とx軸 の交点 x 解の1つ 接線と x 軸の交点は解の1つに近づいたと期待 関数上の点

x 接線と x 軸の交点 接線 y = f '(x0)(x-x0) + f(x0) 関数上の点 から 接線とx軸の交点のx座標 (x0 - f(x0)/f '(x0), 0) 接線 y = f '(x0)(x-x0) + f(x0) x 関数上の点 　　　　　　から 接線とx軸の交点のx座標 を求めることを繰返す 関数上の点 (x0, f(x0))

ニュートン法の例 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...

x 接線と x 軸の交点 (0.54545454..., 0) x1 = 0.54545454... 関数上の点 (0, -6)

接線と x 軸の交点 (0.84895321..., 0) x2 = 0.84895321... x 関数上の点 (0.54545454..., -1.62283996...) x1 = 0.54545454... 収束するまで xiの計算を繰り返す

ニュートン法の繰り返し処理 「収束」したなら⇒再帰終了：現在の xn が解 そうで無ければ： 次の近似値xn+1の値「 xnでの接線と，x軸の交点のx 座標値」を使い再帰 ここで「収束」とは? ニュートン法では現在の xn がどれだけ真の値に近いかは正確には分からない⇒今回は次の判定 ある小さな正数δに対し次が成り立つか?　 成立した時点で再帰を終了し xn を解とする なら true，さもなくば false は次のように書く (abs は「絶対値」を求める関数，TOLERANCE はδ) (<= (abs (f xn)) TOLERANCE)

h 接線の傾きを求める ;; 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))) ;; 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)) (x+h,f(x+h)) h f(x+h)-f(x-h) (x-h,f(x-h)) x+2h 接線と x 軸の交点のx座標を求める ニュートン法による再帰関数

補足. 接線の傾き 接線の傾きを求める関数 d/dx 関数 f と座標値 x, 近似幅 hから，x における f の傾き (つまり f ' (x) ) の近似値を求める x f(x) 傾きは 接線 0 に近い h の値で、 接線の傾きf '(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) x 収束対策 繰り返しの上限回数 を設定→課題 虚数解は求まらない 求まる解はあくまで「近似解」 近似解→どの程度の精度で計算するか(やめるか)により、収束時間、丸め誤差の影響など変化 初期近似値の選び方によって 求まる解が変わり得る 収束に問題：(解と離れていると)時間がかかる、収束しない場合がある、など １ ２ ３ x 関数ｆ(x)が単調でなく変曲点を持つ（つまりｆ’(x)の符号が変わる）とき 例）上の図の1から開始 　　2 (f'(x) = 0となる点)が選ばれる 　　3 ： y 軸との交点が求まらない （負の無限大に発散） →　収束しない 収束対策 繰り返しの上限回数 を設定→課題

課題

課題1. クイックソート 例題2のクイックソートプログラムについて，pivot と比較してリストを分割する際の filter1 の比較演算に < ではなく <= を使ってしまった場合(他の部分はそのまま)，停止性に問題が生じる．(list 5) をソートする場合を例に，その理由を述べよ． 次のような構造体 AddressNote のデータを要素とするリストを考える．そのリスト要素のAddressNote データを名前 nameの昇順にソートするクイックソートのプログラムを作れ． (define-struct AddressNote (name age address)) ここで name と address は文字列，age は自然数

課題2. Newton 法 例題4 Newton 法プログラムを次のように変更せよ 精度を変更できるよう， TOLERANCE を変数定義ではなく関数 newton への引数 t として与える． 収束しない場合でも停止するように変更する． 最大の繰り返し回数を引数として指定でき，指定回数に達した場合は値が収束しなくても再帰を停止するようなプログラムを作成せよ．