タグ「一人読書会」が付けられているもの

 

1.2.1 線形再帰と反復

問題 1.9

次の2つの手続きはどちらも引数を 1増やす手続きと、引数を 1減らす手続きdec をつかって2つの正の数を足す方法を定義している。

(define (plus a b)
  (if (= a 0)
  b
  (inc (plus (dec a) b))))
(define (plus a b)
  (if (= a 0)
  b
  (plus (dec a) (inc b))))

置き換えモデルを使い、それぞれの手続きが (plus 4 5) を評価するときに生成するプロセスを示し、反復的か再帰的かを述べよ。

 

まz、「置き換えモデル」 1.1.5 手続き作用の置き換えモデルで出てきた、各仮パラメータを対応する引数で取り替え、手続きの本体を評価するプロセス。手続きを考える時の補助であり、解釈系の実際の働きを述べるものではない。

やってみよう。

(define (plus 4 5)                  (if (= 4 0) 5 (inc (plus  3 5))))
    (define (plus 3 5)              (if (= 3 0) 5 (inc (plus  2 5))))
        (define (plus 2 5)          (if (= 2 0) 5 (inc (plus  1 5))))
            (define (plus 1 5)      (if (= 1 0) 5 (inc (plus  0 5))))
                (define (plus 0 5)   5) 
            (define (plus 1 5)       6)
        (define (plus 2 5)           7)
    (define (plus 3 5)               8)
(define (plus 4 5)                   9)
9 

(define (plus 4 5)                 (if (= 4 0) 5 (plus  3 6)))
    (define (plus 3 6)             (if (= 3 0) 6 (plus  2 7)))
        (define (plus 2 7)         (if (= 2 0) 7 (plus  1 8)))
            (define (plus 1 8)     (if (= 1 0) 8 (plus  0 9)))
                (define (plus 0 9) 9)))
            (define (plus 1 8)     9)))
        (define (plus 2 7)         9)))  
    (define (plus 3 6)             9)))
(define (plus 4 5)                 9)))
9

本書で、紹介されている書き方で表現するのはちょっと困難なので、再帰をインデントで表現してみた。

前回確認したのと同様、前者が再帰的プロセスにより結果が膨張と収縮にて生成されており、後者は反復的プロセスにより結果が生成されている。

後者のプロセスは末尾再帰的というらしいのだが、何故か?と思ったら、関数の最後の(末尾の)ステップが再帰呼び出しになるような形になっていることによるらしい。上記の置き換えモデルがどのように処理されるかを見れば分かるが、再帰的プロセスと異なり、反復的プロセスにでは、計算結果が出た後は、値を呼び出しの逆順に返すだけだ。このため、末尾再帰的なプロセスは、呼び出しをジャンプに最適化できるということらしい。

Wikipedia 末尾再帰 を参照。

今日はここまで。

 

1.2.1 線形再帰と反復

反復的プロセス(iterative process)

再帰的プロセス(recursive process) の確認を前回行った、次は、同じく階乗を求める処理を反復的プロセス(iterative process) で行う。

(define (factorial n)
  (fact-iter 1 1 n))

(define (fact-iter product counter max-count)
  (if (> counter max-count)
      product
      (fact-iter (* counter product)
                 (+ counter 1)
                 max-count)))

・・・反復という響きから、再帰ではなく、繰り返しを行うイメージだったが、例示されているコードは、思いっきり再帰している様に見える。

そもそもSchemeでは、ループが書けないのか?

反復的プロセスが『再帰』してしまっていたら、再帰的プロセスと、反復的プロセスの違いは何なのだ?

 

前回の Java での書き換えはメソッド呼び出しが見えにくかったので、今回はもう少しわかりやすく、メソッドの入り口と出口で、内容をトレースする様にして、再帰的プロセス、反復的プロセス、ループ処理のそれぞれで、途中経過を出力しながら階乗を計算させてみる。

package cap1_2_1;

public class Factorial {

    public static void main(String[] args) {
        Factorial me = new Factorial();
        long num = Long.parseLong(args[0]);
        
        // 再帰処理
        System.out.printf("%d! = %d\n", num, me.recursiveFactorial(num));
        // 反復処理(再帰)
        System.out.printf("%d! = %d\n", num, me.iterativeFactorial(1L,1L,num));
        // 反復処理(ループ)
        System.out.printf("%d! = %d\n", num, me.loopFactorial(num));
    }
    
    /* 再帰処理により階乗を計算 */
    public long recursiveFactorial(long num) {
        trace(num,"recursive", true);    
        if (num == 1L) {
            return trace(1L, "recursive" , false);
        }
        return trace(num * recursiveFactorial(num - 1L),"recursive" , false);
    }

    /* 反復処理(再帰)により階乗を計算 */
    public long iterativeFactorial(long product, long counter, long num) {
        trace(product,"iterative", true);
        if (counter > num) {
            return trace(product, "iterative", false);
        }
        return trace(iterativeFactorial(product * counter, counter + 1, num), "iterative", false);
    }
    
    /* 反復処理(ループ)により階乗を計算 */
    public long loopFactorial(long num) {
        trace(num,"loop", true);

        long product = 1L; 
        long counter = 1L;
        
        while(true) {
            if (counter > num) {
                return trace(product, "loop", false);
            }
            product *= counter;
            counter++;
        }
    }

    /* 
     * メソッドの呼び出しをトレース 
     * メソッドに入った時、リターン時に呼ぶ
     */
    private long trace(long num, String method, boolean isEnter) {
        /*
        // スタックトレースを吐く
        StackTraceElement[] elms = Thread.currentThread().getStackTrace();
        for (StackTraceElement e : elms) {
            System.out.printf("\t%s ", e.getMethodName());
        }
        */
        // (E) Enter (R) Return 
        System.out.printf("\t(%s) %s %6d\n", (isEnter)?"E":"R", method, num);
        return num;
    }
}

8 を引数に処理を行ったところ。(E) は、メソッドの入り口、(R)は出口でのトレース。

recursive が再帰プロセスで、iterative が反復プロセス。

	(E) recursive      8
	(E) recursive      7
	(E) recursive      6
	(E) recursive      5
	(E) recursive      4
	(E) recursive      3
	(E) recursive      2
	(E) recursive      1
	(R) recursive      1
	(R) recursive      2
	(R) recursive      6
	(R) recursive     24
	(R) recursive    120
	(R) recursive    720
	(R) recursive   5040
	(R) recursive  40320
8! = 40320
	(E) iterative      1
	(E) iterative      1
	(E) iterative      2
	(E) iterative      6
	(E) iterative     24
	(E) iterative    120
	(E) iterative    720
	(E) iterative   5040
	(E) iterative  40320
	(R) iterative  40320
	(R) iterative  40320
	(R) iterative  40320
	(R) iterative  40320
	(R) iterative  40320
	(R) iterative  40320
	(R) iterative  40320
	(R) iterative  40320
	(R) iterative  40320
8! = 40320
	(E) loop      8
	(R) loop  40320
8! = 40320

おー 再帰、反復の両プロセスでも、『再帰』を利用しているが、行われている処理内容が全く異なることが分かった。

再帰プロセスでは、行きに計算処理をスタックに積んでいき(膨張)、帰りに計算を行っている(収縮)が、反復処理では徐々に成長していっている。

本書に「再帰的プロセス(process)と再帰的手続き(procedure)を混同しないように注意しなければならない」とあるのは、この違いのことか。

「手続きが再帰的」 というときには、手続の構文に再帰を含むことをいう。上記Java に書き換えた iterativeFactorial() メソッドは、再帰的手続きだが、反復的プロセスということ。なるほど。

「(Ada,Pascal,Cを含む)通常の言語の実装が、再帰手続きの実行で消費する記憶量が、プロセスは原理的に反復的であっても、手続き呼び出しの数とともに増加するように設計してある」 要するに、通常の言語では、本質的にプロセスが反復的でも、構文が再帰だと、メソッド呼び出しがスタックされる。このために、これらの区別がまぎらわしいとのこと。

また、『これらの言語では、反復プロセスは、do、repeat、until、for、while のような特殊目的の「ループ構造」に頼ってしか記述できない。』

確かに、上記 Java 例で iterativeFactorial はメソッド呼び出しが、一番底まで到達すると結果は計算されていて、あとはメソッド呼び出しを戻りながら結果を呼び出し元に返すだけにもかかわらず、メソッド呼び出しはスタックに積まれていくのだから、実際のプロセスが反復的だからといってメモリの使用量など変わらないだろう。trace() メソッドのスタックトレースを吐くコメントを外すと同じようにメソッドがスタックにつまれていって取り出されていく様がよく分かる。

繰り返し構文で書き直すと、loopFactorial() メソッドとなる。

じゃあ、Scheme ではどうなのよ?と思ったら、5章で検討すると。また、実装上この欠点はないと。

へー

構文上の再帰は、シンタックスシュガーで、実際は固定スペースで実行できるらしい。この性質の実装を 末尾再帰的(tail recursive) と言うらしい。

へー

今日はここまで。

 

1.2.1 線形再帰と反復

再帰的プロセス(recursive process)

任意の正の整数 n に対する階乗関数 n! は、 n ・ (n-1)! で得られる。つまり、n! は、 (n-1)!  を計算し、結果に n をかけて計算できる。1! は 1に等しいので、以下の手続きに変換できる。

(define (factorial n)
  (if (= n 1)
      1
      (* n (factorial (- n 1)))))

この形のプロセスを再帰的プロセス(recursive process) といい、遅延演算(deferred operations) の列で特徴づけられる。

遅延演算とは、プロセスが膨張と収縮の形をとるとき、膨張時に遅延演算の列として作成される。収縮は演算が実際に実行される時に起きる。

解釈系は、後に実行する演算を覚えておく必要があり、覚えておく必要のある情報量は n に比例して(線形に)成長するため、こういうプロセスを線形再帰的プロセス(liner recuresive porceee) という。

例によって、Schemeでは、その様子を確認するだけのスキルがまだ無いので、たまには Java で確認してみる。

やっていることの、本質的には、これだけのこと。

private long factorial(long num){
package cap1_2_1;

public class Capter1_2_1 {
	public long factorial(long num){
		if (num == 1L){
			return 1L;
		}
		return num * this.factorial(num - 1L);
	}
}

処理の経過をトレースするために、いろいろと組み込み、引数に 8 を与えて実行してみる。

package cap1_2_1;

import java.util.Stack;

public class Capter1_2_1 {
	// トレース用
	private Stack stk = new Stack();
	
	public static void main(String[] args) {
		Capter1_2_1 me = new Capter1_2_1();
		long num = Long.parseLong(args[0]);
		me.testFactorial(num);
	}
	public void testFactorial(long num){
		stk.push(String.format("factional(%d)", num));
		long result = this.factorial(num);
		System.out.printf("*************\n");
		System.out.printf("%d! = %d\n", num, result);
	}
	public long factorial(long num){
		System.out.println(stk.toString());
		stk.pop();
		stk.push(String.format("%d", num));
		
		if (num == 1L){
			System.out.println(stk.toString());
			return 1L;
		}
		
		stk.push(String.format("factorial(%d)", num - 1L));
		long num_minus_1 = this.factorial(num - 1L);
		
		stk.pop(); // 自身
		stk.pop(); // ひとつ手前を計算結果に置き換え
		stk.push(String.format("%d", num * num_minus_1));
		
		System.out.println(stk.toString());
		return num * num_minus_1;
	}
}

膨張と、収縮の様子がよく分かる。

[factional(8)]
[8, factorial(7)]
[8, 7, factorial(6)]
[8, 7, 6, factorial(5)]
[8, 7, 6, 5, factorial(4)]
[8, 7, 6, 5, 4, factorial(3)]
[8, 7, 6, 5, 4, 3, factorial(2)]
[8, 7, 6, 5, 4, 3, 2, factorial(1)]
[8, 7, 6, 5, 4, 3, 2, 1]
[8, 7, 6, 5, 4, 3, 2]
[8, 7, 6, 5, 4, 6]
[8, 7, 6, 5, 24]
[8, 7, 6, 120]
[8, 7, 720]
[8, 5040]
[40320]
*************
8! = 40320

Java 好きだし、おもしろいけど、やはりPython と比べると、思考から実装までのラグがちょっと大きいなぁ~

今日はここまで。

 

1.1.8 ブラックボックス抽象としての手続き

 

新しい言語を覚えようとすると、まず変数の宣言の仕方や、制御構文を覚えるが、1つなにか知っていればそれらを覚えるくらいは何のことはない。

しかし、パラダイムが異なるとこれは大変。多分プログラム言語で習得が難解とされていることがらは、パラダイムの相違に起因しているのではないかと、この章を読んで思った。

例えば、C言語のポインターなんて、高級言語であるC言語に、アセンブラのパラダイムが持ち込まれているのが難解といわれる原因だろうし、構造化言語に浸かっていればオブジェクト指向はちんぷんかんぷんだろう。ただし、逆からみれば、あまりにも自明過ぎて、何が難解か分からないだろう。

なぜ、そんなことを思ったかというと、ちょっと前に、クロージャー関数がはやった(?)が、当初何のことやら分からなかった。また、レキシカル変数とか言う用語の意味もよく分からなかった。そこそこプログラムを書いていたにもかかわらず。

それらの概念が、自分が主にやってきた言語には無かったから、理解するのに手間取ったが、この章で述べられている用語をきちんと理解していれば、そのどちらも悩むことなく理解できただろうなぁと。

自分の興味の変遷は Java → JavaScript → Python → Scheme  な感じだが、今思うと何というか自然に思える。

手続き抽象(procedural abstraction)

 分割統治法を用いるとき、部分の手続きがどう結果を計算するのかには関心を持たない。Wikipediaでは、「アルゴリズムは自由に選択できる」とされている。それは「手続き」ではなく、「手続き抽象(procedural abstraction)」である。

局所名

手続きの仮パラメータは、どんな名前でもかまわない。手続き本体に対して局所的で無ければならない。

束縛変数(bound variable)

手続き定義のなかで、仮パラメータはどんな名前でもかまわないし、名前をすべて変更しても手続きの意味は変わらないという意味で、手続き定義は、仮パラメータを束縛している(bind)。そういう名前を束縛変数(bound variable)という。変数が束縛されていなければ、自由である(free)。

有効範囲(Scope)

名前が束縛されている式の範囲

ブロック構造(block structure)

定義の入れ子。単純な名前保護の機構。ここで作成した、sqrt2 関数から呼び出される関数、average、square、good-enough? ・・・ を sqrt2 関数のブロックの中に定義できる。sqrt2 ブロックの外で、同名の手続きを定義しても、それぞれには影響を与えない。

(define (sqrt2 x) 
  (define (average x y) (/ (+ x y) 2))
  (define (square x) (* x x))
  (define (good-enough? guess x)
    (< (abs (- (square guess) x)) 0.001))
  (define (improve guess x)
    (average guess (/ x guess)))
  (define (sqrt-itr guess x) 
    (if (good-enough? guess x)
        guess
        (sqrt-itr (improve guess x)
                  x)))
  (sqrt-itr 1.0 x))

静的有効範囲(lexicacl scoping)

上記例で、x は、sqrt2 の定義に束縛されている。なので、その内側の関数は、x の有効範囲内にあるため明示的に渡さなくてもよいので、自由変数にできる。こうしたやり方を静的有効範囲(lexical scoping) という。

(define (sqrt2 x) 
  (define (average a b) (/ (+ a b) 2))
  (define (good-enough? guess)
    (< (abs (- (* guess guess) x)) 0.001))
  (define (improve guess)
    (average guess (/ x guess)))
  (define (sqrt-itr guess) 
    (if (good-enough? guess)
        guess
        (sqrt-itr (improve guess))))
  (sqrt-itr 1.0))

今日はここまで。

 

1.1.7 Newton 法による平方根

 

問題 1.8

立方根をとるNewton法は y が x の立方根の近似値なら、よりよい近似は以下の値で与えられるという事実によっている。この手続きを実装せよ。

newton06_01

(define (good-enough? guess pre-guess)
  (< (abs (- guess pre-guess)) 0.001))
(define (improve guess x)
  (/ (+ (/ x (* guess guess)) (* 2 guess)) 3))
(define (cuberoot-itr guess x) 
  (if (good-enough? (improve guess x) guess)
      (improve guess x)
      (cuberoot-itr (improve guess x)
                x)))
(define (cuberoot x) (cuberoot-itr 1.0 x))

newton07_01

こんな感じかな。

 

1.1.7 Newton 法による平方根

 

問題 1.8

立方根をとるNewton法は y が x の立方根の近似値なら、よりよい近似は以下の値で与えられるという事実によっている。この手続きを実装せよ。

newton06_01

まず、平方根同様に、式を確認

newton06_02

納得いった。

いきなりSchemeでかけるほど訳がわかっていないので、まずは、イメージをつかむために、Python でグラフを書いてみる。5 の立方根を計算。

# -*- encoding:utf-8 -*-
'''
Newton 近似により、平方根を求めるサンプル
'''

from pylab import *
from numpy import *
from numpy.ma.core import arange
from decimal import *

def d(x):
    return Decimal(x)

cube = lambda x, b : d(x) ** d('3') - d(b)         
newton = lambda target, guess : (d(target) / (d(guess) ** d('2')) + (d('2') * d(guess))) / d('3')

def f1(x_range, target, interval):
    x = arange(d(x_range[0]), d(x_range[1]), d(interval))
    y = [cube(y1, target) for y1 in x]

    return (x, y)

if __name__ == '__main__':
    grid(True)
    getcontext().prec = 8
    target = d('5')
    
    # 3次元のグラフ
    x, y = f1(('-6','8'),target, '0.1')
    plot(x,y)

    guess = d('6.0')
    repeat = 8
    for i in range(repeat):
        x = (guess, guess)
        y = (d('0'), cube(guess, target))
        plot(x, y, 'k--')

        new_guess = newton(target, guess)
        
        x = (guess, new_guess)
        y = (cube(guess, target), d('0'))
        plot(x, y)
    
        guess = new_guess
        print guess

    show()

newton06_03

今日はここまで。

次は、これを Scheme に移植してみる。

 

1.1.7 Newton 法による平方根

 

問題 1.7

1.ある繰返しから guess の変化に注目し、変化が予測値に比べ非常に小さくなった時点で止める手続きを設計する。これは小さい数、大きい数に対してうまく動くか?

(define (good-enough? guess pre-guess)
  (< (abs (- guess pre-guess)) 0.001))
(define (average x y) (/ (+ x y) 2))
(define (improve guess x)
  (average guess (/ x guess)))
(define (sqrt-itr guess x) 
  (if (good-enough? (improve guess x) guess)
      (improve guess x)
      (sqrt-itr (improve guess x)
                x)))
(define (sqrt2 x) (sqrt-itr 1.0 x))

まず、単純に予測値の差が一定の値(0.001) 未満になったら、処理を終了するように good-enough? を書き換え、sqrt-itr の実装をそれにあわせて変更してみる。

問題があった値に対して、確認してみる。

newton05_01

小さい値で、誤差が非常に大きくなってしまったり、大きい値で無限ループに陥ることはなくなった。

また、good-enough? での判定用定数を小さくすることで、精度もあがるようだ。

sqrt-itr で、improve が 2回 呼ばれるのが気になるが、まだ、Scheme で変数の使い方を知らないので、よしとしておく。

一般的には、終了判定は、以下とするようだ

newton05_02

「これ以外にも計算の終了を決めることは可能なので、状況に応じて、計算の打ち切り方法を決めればよい。」 とのこと。

大きい数、小さい数に対しても、うまく動く。・・・ と思うんだけどなぁ。「それなりにうまく動く」か?

 

1.1.7 Newton 法による平方根

 

問題 1.7

さて、ようやく問題の内容自体を理解できる準備が整ったので、考えてみる。問題の概要。

(define (square x) (* x x))
(define (good-enough? guess x)
  (< (abs (- (square guess) x)) 0.001))
(define (average x y) (/ (+ x y) 2))
(define (improve guess x)
  (average guess (/ x guess)))
(define (sqrt-itr guess x) 
  (if (good-enough? guess x)
      guess
      (sqrt-itr (improve guess x)
                x)))
(define (sqrt2 x) (sqrt-itr 1.0 x))

まず、動きを検証するのに、デバッグ文を仕込んだりしたいのだが、Scheme でどうするのかさっぱりわからないので、Python に上記を移植してみる。

# -*-encoding:utf-8-*-
'''
http://ja.wikipedia.org/wiki/IEEE_754
'''
from decimal import *

def init_env():
    '''環境の設定'''
    global use_decimal   # Decimal を利用するか
    global judge         # 計算打ち切り判定
    global initial_guess # 推測初期値
    
    use_decimal = True  
    getcontext().prec  = 6          
    judge              = dcm('0.00001') 
    initial_guess      = '1.0'         
    
def dcm(x):
    '''数値を扱う型を切り替える float, decimal'''
    if use_decimal:
        return Decimal(x)
    else:
        return float(x)
def square(x):
    return dcm(x) * dcm(x)
def good_enough(guess ,x):
    g2 = square(guess)
    v = g2 - dcm(x)
    if v < dcm('0'):
        v = v * dcm('-1')
    print '*** %s -  %s = %s ***' % (g2, dcm(x), v)
    return v < judge
def avarage(x, y):
    return (dcm(x) + dcm(y)) / dcm('2')
def improve(guess, x):
    av = avarage((dcm(x) / dcm(guess)), guess) 
    print av
    return av
def sqrt_itr(guess, x):
    if good_enough(guess, x):
        return dcm(guess)
    else:
        return sqrt_itr(improve(guess, x), x)
def sqrt2(x):
    init_env()
    return sqrt_itr(initial_guess, x)

if __name__ == '__main__':
    target = '0.00003'
    print sqrt2(target)

init_env 関数に、設定を行うようにしている。

1.good-enough? テストは、非常に小さい数の平方根をとるときには効果的ではない。

good-enough? テストは、誤差が 0.001 という定数 以下になったら、計算を打ち切るようになっている。

まず、定数になっていることから、直感的に非常に小さい値の場合、誤差の割合が相対的に大きくなると思われる。

また、判定基準が 0.001 となっているが、この値を比較するのは、平方根をとる前(とろうとしている)値のため、求める値が1以下の値の場合平方根もさらに小さな値になり、判定基準の値、0.001 が、求める値に対して相対的にかなり大きくなってしまう。

newton04_01

    use_decimal = True  
def init_env():
    (略)
    use_decimal = True                       # Decimal を使う  
    getcontext().prec  = 6                   # Decimal の有効桁数を 6桁に
    judge              = dcm('0.00001')      # 0.001 より精度をあげる
    initial_guess      = '1.0'               

if __name__ == '__main__':
    target = '0.0003'
    print sqrt2(target)

上記 Python の init_env の設定を上記にして検証してみる。

*** 予測値 - 平方根をとる値 = good-enough? の判定値

*** 1.00 -  0.0003 = 0.9997 ***
0.50015
*** 0.250150 -  0.0003 = 0.249850 ***
0.250375
*** 0.0626876 -  0.0003 = 0.0623876 ***
0.125786
*** 0.0158221 -  0.0003 = 0.0155221 ***
0.0640855
*** 0.00410695 -  0.0003 = 0.00380695 ***
0.0343834
*** 0.00118222 -  0.0003 = 0.00088222 ***     ← ここで判定OKとなってしまう。
0.0215542
*** 0.000464584 -  0.0003 = 0.000164584 ***
0.0177363
*** 0.000314576 -  0.0003 = 0.000014576 ***
0.0173254
*** 0.000300169 -  0.0003 = 1.69E-7 ***
0.0173254                                             

かなり(0.001より)精度が悪い結果になってしまっている。

 

2.限られた精度で実行される場合、非常に大きい数にも不適切。

翻って、非常に大きい値の場合。おそらく、内部数値表現の精度によって、対象が大きい値だと 0.001 という精度自体を保てないためだと思われるが、実際どうなのか確認してみる。

newton04_02

どんどん値を大きくしていくと、結果がかえってこなくなった。

Python で確認してみる。

    use_decimal = True  
def init_env():
    (略)
    use_decimal = False                      # Decimal は使わない(floatとなる)
    judge              = dcm('0.001')  
    initial_guess      = '1.0'               

if __name__ == '__main__':
    target = '33333333333333333333333333333333333'
    print sqrt2(target)

 

やはり、ある程度以上は、予測値が近づかず、変化しないため、無限ループになってしまっているようだ。

    (略)

*** 3.33333333333e+34 -  3.33333333333e+34 = 4.61168601843e+18 ***
1.82574185835e+17
*** 3.33333333333e+34 -  3.33333333333e+34 = 4.61168601843e+18 ***
1.82574185835e+17

    (同じ処理の繰り返し・・・)

Traceback (most recent call last):
  File "C:\test.py", line 50, in <module>
    print sqrt2(target)
  File "C:\test.py", line 46, in sqrt2
    return sqrt_itr(initial_guess, x)
  File "C:\test.py", line 43, in sqrt_itr
   (略)
RuntimeError: maximum recursion depth exceeded while calling a Python object

Python で、 計算に利用する数値データ型を、float から Decimal に変えて、十分に桁数をとれば、処理は完了するので、あくまで 内部数値データ型の問題だろう。

 

1.1.7 Newton 法による平方根

 

問題 1.7

さて、Newton 法 も何となくわかった気になったが、本書に書かれているロジック とは何か違う。

前回は なんや、微分して、x軸との交点を求めて・・・ みたいなイメージだったのだが、本書では、予測値で、平方根を出したい値を割って、さらにその平均をとることを繰り返している。

√2 を計算する

予測値 平均
1 2/1 = 2 (1 + 2) / 2 = 1.5
1.5 2/1.5 = 1.3333 (1.5 + 1.3333) /2 = 1.4167
1.4167 ・・・ ・・・

うーん 違和感を感じる。何故平均でいいんだ?

最近、「数学ガール」 を読了して、数式をみて思考停止しなくてすむようになったのと、ちょっと数式を展開するのがおもしろく感じたこともあるし、今年は、中学程度も怪しい数学レベルをせめて高校程度までもっていこうという野望を抱いているので、ここはちょっとこだわっておこう。そのうち息子に算数を教えることができるようにしとくためにも・・・

まず、上記で行っている作業は、平方根を求める対象を b として、予測値を a とし、素直に以下の式とする。

newton02_01

で、一方、Newton 近似の 一般式は、

exp  でした。

本書で言う、予測値の部分が、an1 箇所なので、やはり平方根を求める対象を b とすると、

newton02_02 となるので、

newton02_03 2a で通分して、

newton02_04 分子を計算して、

newton02_05 分母のa を消すと・・・

newton02_06 お、一致しましたね。

これでようやく、本書で実装されたところの意味を理解することができました。

 

一人読書会 「計算機プログラムの構造と解釈」 (2)

1.1.7 Newton 法による平方根

Newton法 により平方根を求める手続きを実装する。

問題 1.7

平方根の計算で使った good-enough? テストは非常に小さい値の平方根をとるときには効果的ではない。また限られた精度において非常に大きい数にも不適切である。それぞれどのようにテストが失敗するか?また、変化が予測値にくらべ非常に小さくなった時に終了する手続きを設計せよ。

 

まず、残念ながら、Newton 法 を知らないので、今年の目標は数学の勉強をすることなので、読み流さずにWikipediaで調べることとする。

重ねて残念ながら、Wikipedia の説明では、知識不足で、わかったような、わからないような。

偶然、所有している本「微分・積分のしくみ」に、ニュートン近似として、丁寧な説明があったので、まずそちらを参考にして、Newton法による平方根の近似とはどんなものかの理解を試みる。

newton01

たとえば、√5 の値を調べるとする。

2次関数のグラフを書く

まず、上の図の2次関数のグラフは、y = x^2 - 5 のグラフ。

√5  を 満たすx、 x = √5  は、 x^2 = 5 であり、 x^2 - 5 = 0 となる。

Python で、

exp  = lambda x : x**2 - 5

と定義してみる。

適当な開始位置を決める

求めたいのは、上記の2次関数のグラフで、y=0 となるときの x の値。グラフを見ると、0 ~ 5 の間にあるので、とりあえず開始位置を x = 10 としてみる。すると、x = 10 の場合、y = 95 に値が定まる。(上図①)

接線を求め、x軸との交点を求める

y = x^2 - 5 のグラフ上の、上記で決めた (10, 95) での接線が、x1 で、y = 0 となる接線の方程式を求める

y = f(x) = x^2 - 5  を 微分すると  f'(x) = 2x で、f'(10) = 20 が 接線の傾きになる。

とすると、

接線の方程式 : y - f(a) = f'(a)(x - a)

y - 95 = 20x - 200

y = 20x - 105

( y = 20x + b 、 95 = 200 + b 、 b = -105 )

この接線が、y = 0 となるのは、x = 5.25 となる (上図②)

これで、みてわかるように、最初の値、x = 10 より、y = x^2 - 5 の x軸との交点により近づいた。

ここで取得した値、5.25 を初期値として、上記の処理を繰り返すことで、√5 の値に近似していく。

というのが ニュートン近似だと。

上図の ① を an1 ② を an 、として、一般化すると、exp となる。

計算結果を出力してみる。

ようするに、n をどんどん増やしていくと、徐々に近づいていっているのがわかる。とりあえず、富士山麓に鸚鵡が鳴く条件は満たした。

result01

result02

グラフのコードなど

Python で グラフを書くための準備については、以下を参照

Python で グラフをかけるようにしとくと、こういった勉強の時に役立つなぁ。と。簡単にかけるし。Python イイ!

/tips/wiki.cgi?page=Python+NumPy

/tips/wiki.cgi?page=Python+matplotlib

上記のコード

# -*- encoding:utf-8 -*-
from pylab import *
from numpy import *
from numpy.ma.core import arange

# 二次関数 f(x) y = x^2 - 5
exp  = lambda x : x**2 - 5
# 二次関数の微分 f'(x) y = x^2 - 5 
expd = lambda x : 2 * x
        
def f1(x_range, intarval):
    '''
        関数expに対して、xの値 x_range[0] ~ x_range[1] 間隔 intarval に対応する y を取得
    '''
    x = arange(x_range[0], x_range[1], intarval)
    y = [exp(tx) for tx in x]
    return (x, y)

if __name__ == '__main__':
    grid(True)
    
    # 二次関数のグラフを書く
    (x, y) = f1((-11.0, 11.0),0.1)
    plot(x, y)

    # 適当な値を設定
    tmp_x2 = 10.0
    
    repeat = 5
    for i in range(repeat):
        print '%d : %s' % (i+1, tmp_x2)
        tmp_x1 = tmp_x2
        
        x1 = (tmp_x2, tmp_x1) 
        y1 = (0, exp(tmp_x1))
        plot(x1, y1, 'k--')
    
        tmp_x2 = tmp_x1 - (exp(tmp_x1) / expd(tmp_x1))
        x1 = (tmp_x2, tmp_x1) 
        y1 = (0, exp(tmp_x1))
        plot(x1, y1)
      
    show()

ようやく、問題に取りかかる準備ができた。

このペースだと読了するのに1年ではきかんな~。

一人読書会 「計算機プログラムの構造と解釈」 (2)

1.1.7 Newton 法による平方根

Newton法 により平方根を求める手続きを実装する。

問題 1.6

ifは特殊形式であるが、cond を利用して普通の手続きとして定義し、平方根の計算にこれを使おうとすると何が起きるか。

 

以下、特殊形式の if を 通常の手続きである、new-if を定義し、平方根を Newton法 で求める処理を書き直した。

(define (square x) (* x x))
(define (sqrt-itr guess x)
  (new-if (good-enugh? guess x)
      guess
      (sqrt-itr (improve guess x)
                x)))
(define (improve guess x)
  (avarage guess(/ x guess)))
(define (avarage x y)
  (/ (+ x y) 2))
(define (good-enugh? guess x)
  (< (abs (- (square guess) x)) 0.001))
(define (sqrt2 x)
  (sqrt-itr 1.0 x))

(define (new-if predicate then-clause else-clause)
  (cond (predicate then-clause)
        (else else-clause)))

実際に動かしてみると、スタックオーバーフローだろう、Out Of Memory が発生した。

特殊形式 if を使った場合には、想定通りの動きをしていた。ぱっと見、等価に思えるが。

sicp_q_1_6

うーーん 。。。 なぜだろう 。。。

・・・

特殊形式・・・ define は 特殊形式で、作用は行わない。。。 if も特殊形式・・・

通常評価は、部分式を評価し、作用(演算子を適用)させる・・・

おお!そうか、特殊形式ではないから、被演算子がすべて評価されてしまうからか!、sqrt-itr の else 節 で、再帰しているのが、終了条件に合致しても、else節 がなんと、評価されてしまうのだ!

if は特殊形式で、条件に合致した then または、 else のどちらかしか評価しないのだ。

なるほど。

一人読書会 「計算機プログラムの構造と解釈」 (1)

1.1.3 組合せの評価

組合せは、以下のように評価される.

  1. 組合せの部分式を評価
  2. 最左の演算子を被演算子に作用させる

これを 作用的順序の評価(applicative-order evaluation) というが、もう一つ、組合せを完全に展開してから演算を行う評価方法を、正規順序の評価(normal-order evaluation) という.

問題 1.5

作用的順序の評価を使う解釈系での振る舞い、正規順序の評価を使う解釈系での振る舞いを説明せよ

>(define (p) (p))
>(define (test x y)
  (if (= x 0)
      0
      y))
>(test 0 (p))

if がショートサーキットであるとするならば、解釈系が、作用的順序であれば、(test 0 (p)) で、最初に x = 0 が評価されて、結果が 0 となり、test の第2引数 (p) は評価されない。

正規順序の評価であれば、まず組合せが展開される段階で、test の第2引数 (p) が無限ループ(?) に陥ってしまう。

・・・ のではないか。

ちなみに、DrScheme で実行すると、結果が返ってこなかった。

「計算機プログラムの構造と解釈」 を購入。しっかりと読み込みたいので、1年くらいかけて一人読書会としゃれ込もうと思う。

・・・電車で読むだけではもったいない。

1.手続きによる抽象の構築

まず、全編を通して、プログラムの説明には、C や Java や Python ではなく、Lisp の方言である、Scheme (スキーム) を使用する。

そもそも、Scheme の設計者である、ジェラルド・ジェイ・サスマン が、本書の著者。

amazon のレビューをみると、訳に対する苦情が多いようだが、自分も最初のページで違和感を感じた。デバッグのことを「虫とり」とか。

ただ、訳者が、1931年生まれの先生ということを知ると、訳の雰囲気の違和感も哲学書のような雰囲気を醸し出してきた。自分はこういう訳も嫌いではない。

で、本書では、 cheme がわからないと話にならないので、最初の章は、文法および考え方の説明に費やされている。

こちらに簡単にまとめた。さらにまとめると、

確認用の処理系には、DrScheme を使用することとする。

手続きの作用を表現。括弧で囲んで、組合せと呼ぶ。

> (+ 1 2)
3

この例だと、加算(+) が、手続きを表す演算子(オペレーター)、1、2 が作用される被演算子(オペランド)

演算子を左端に書く記法を、前置記法という。利点として、任意個数の引数をとることができる、組合せを入れ子にできる。

特殊形式(special forms)

引数に作用しない場合、組合せとは呼ばず、特殊形式という。たとえば、define

> (define x 2)

これで、変数 x に 2 をセットする。

合成手続き

define をつかって、演算に名前をつけることができる。自乗を行う演算に名前をつけてみる。

> (define (square x) (* x x))
> (square 2)
4

条件式と述語

cond ・・・ 場合分け

if ・・・ 場合が2つの場合

論理演算子

and、or、not

用例は、こちら を参照。

問題

問題1.1

入力して確認。

問題1.2

これでよいと思うけど・・・

> (/ (+ 5 4 (- 2 (- 3 (+ 6 (/ 4 5))))) 
     (* 3 (- 6 2) (- 2 7)))
-37/150

問題1.3

いままで経験した言語との考え方の違いを思い知る。

今の知識、思考方法では、こんな感じにしかかけない・・・

変数や配列を使わずに、2番目に大きい値を取得するって困難。

> (define (smaller x y) (if (< x y) x y))
> (define (last x y z) (smaller x (smaller y z)))
> (define (square x) (* x x))
> (define (q1_3 x y z) (+ 
                        (square (if (> x (last x y z)) x (last x y z)))
                        (square (if (> y (last x y z)) y (last x y z)))))
> (q1_3 4 3 2)
25

問題1.4

引数 b が 負の値の場合、符号を反転し絶対値を取得し、a と加算する

おーなんと、被演算子に演算子を置くと、その演算子が作用するのか!

> (define (a-plus-abs-b a b)
    ((if (> b 0) + -) a b))
> (a-plus-abs-b 2 -3)
5

おー結構しんどいなー。今日はここまで。

ゆっくりでもいいから何とか続けていきたい。。。