カテゴリー

最新の記事

最近のコメント

最近のトラックバック

月別アーカイブ

ブログ検索

RSSフィード

ブロとも申請フォーム

この人とブロともになる

スポンサーサイト

スポンサー広告
--.--.--
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

ベンチャーブログのランキングに参加しています。
下のバナーをクリックして応援していただけると嬉しいです。
にほんブログ村 ベンチャーブログへ

SICP Chapter 1を読了

Scheme/Gauche
2008.12.13
今週の火曜日(2008年12月9日)から「Structure and Interpretation of Computer Programs (SICP)」 を読み始めている。あらかじめ日本語訳の書籍もAmazonで購入しておいたのだが、英文の原著がネット上に全文掲載されているのを見つけたため、英語で読み進んでいる。英語のリーディングの訓練もできるので、一石二鳥だと思ったからだ。

で、昨日、5章構成のうちのChapter 1を読了した。実に面白い。Chapter 1のテーマはプロシジャの抽象化。すでに他のLISPチュートリアルや「プログラミング Gauche」、東京大学の「Scheme 演習」に掲載されていた内容もあったが、私にとって新鮮なネタも数多くあった。まだ足を踏み出したばかりの段階だが、この内容の濃さは素晴らしいの一言である。

例題としては様々な数学ネタが取り上げられていたが、私が特に面白いと思ったのが、フィボナッチ数列の計算法だった。「またフィボナッチ数列かい?」などとおっしゃるなかれ。アルゴリズムによって必要なリソースが大幅に変わってくるということが、非常に明確な形で示されているのだ。

フィボナッチ数列とは何かについては前回触れているので省略。これをごく普通にコード化すると、以下のようになる。(これも前回掲載したが、こちらは再掲させていただく)

(define (fib n)
  (cond ((= n 1) 1)
     ((= n 2) 1)
     (else (+ (fib (- n 1)) (fib (- n 2))))))


さて、このような「再帰型」アルゴリズムの場合、フィボナッチ数の計算量はnの二乗に比例する形で大きくなるらしい。そのためPentium 4/2.4GHzのマシンでは、n=35で2.47秒、n=40で28秒かかるようになり、もはやそれ以上の値では計算したくなくなるのである。

しかしこれを「再帰型」ではなく「反復型」に置き換えると、処理量は一気に減り、nに比例するようになる。

(define (fib n)
 (define (fib-iter a b count)
  (if (< count 2)
     a
     (fib-iter (+ b a) a (- count 1))))
 (fib-iter 1 0 n))


やっていることは、fibの中に反復用の関数を定義し、以下のように処理を進めるように変更しただけである。

Step.1
まず b = 0, a = 1 とする。これは n = 0, n = 1 の時のフィボナッチ数である。

Step.2
与えられたnをカウンター(count)にセット。

Step.3
カウンターの値が2より小さければ、aの値を答えとする。

Step.4
もしカウンターが2以上であれば、以下の処理を行う。

a ← b + a
b ← a
count ← count - 1

Step.5
上記のStep.3~4を繰り返す。

これで n = 35 の場合を計算すると、計算時間は1ミリ秒未満となり、time関数では計測不能になる。ちなみに n = 100000(十万)で、ようやく計算時間が3.585秒となる。つまり計算方法をちょこっと変えるだけで、莫大な効率化を達成できるわけである。

とはいえ、このレベルのリファクタリングは難しいものではない。他のLISP本や「プログラミング Gauche」でも紹介されている。本当に面白いのはここからだ。

実はさらに処理をスケールできる方法がある。それが次の方法だ。基本的な考え方はChapter 1のExercise 1.19に掲載されているのだが、その内容をここで簡単に紹介しておく。

Step.1
まず次のような変換を考える。

a ← bq + aq + ap
b ← bp + aq

なぜこのような式が出てくるのかは、私に聞かないで欲しい。ひょっとすると群論の基本的な知識なのかもしれないが、大学で群論に挫折した私には正直いってよくわからない。それはともかくこの「飛躍」を行うことで、フィボナッチ数の計算をより抽象化できるのである。

フィボナッチ数の計算を行う場合には、p = 0、q = 1 とすればいい。そうすれば上記の式は

a ← b + a
b ← a

となる。

Step.2
上記の変換を「Tpq」と定義しよう。この変換を2回行ったのと等価な変換を「Tp'q'」とすると(つまり「Tp'q'=(Tpq)^2」ということ)、p、q、p'、q'の間には次の関係が成立する。(この関係は計算してみればすぐにわかる)

p' = p^2 + q^2
q' = 2pq + q^2

Step.3
ここまで導けば、ある条件ではTpqではなく、Tp'q'を使うことで、処理量を一気に半分にできることがわかる。つまり「Tpqを実行する回数nが偶数回である」場合には、Tpqをn回実行する代わりに、Tp'q'をn/2回実行することで、目的の値が得られるのだ。n/2 が偶数なら、(Tp'q')^2を求めることで、さらに処理量を半分にできる。後はこれを、nが奇数になるまで繰り返せばいい。

つまり、nが偶数の場合

p ← p^2 + q^2
q ← 2pq + q^2
n ← n/2

を実行するのである。

Step.4
それではnが奇数になったら?その場合は「現時点の」p、qを使い、Tpqを適用した上で、n=n-1とすればいい。なおnが奇数の場合、n=n-1を行った時点で偶数になる。

Step.5
上記のStep.3~4を、nが0になるまで適宜実行する。

これらの処理を組み合わせると、次のような関数が書ける。

(define (fib n)
 (define fib-iter a b p q count)
  (cond ((< count 1)
      b)
     ((even? count)
      (fib-iter
       a
       b
       (+ (* p p) (* q q))
       (+ (* 2 p q) (* q q))
       (/ count 2)))
     (else
       (fib-iter
        (+ (* b q) (* a q) (* a p))
        (+ (* a q) (* b p))
        p
        q
        (- count 1)))))
 (fib-iter 1 0 0 1 n))


これを使ってn=100000(十万)の値を求めると、0.237秒で処理が終わる。つまり単なる反復型より、1桁速くなるわけだ。

ちなみにこのアルゴリズムでは、処理量は log(n)に比例するようになるという。ということは、nが大きいほど効果が大きいというわけである。

n=百万にすると、第3のコードは27秒で答えを出した。反復型では460秒だった。n=十万の時よりは若干差が広がっている。再帰型は気が遠くなりそうなので試していない。

処理内容によってはアルゴリズムをちょっと工夫するだけで、これだけスケーラビリティが変わってしまうのである。再帰型と比べれば、3つ目のコードは2万5000倍もスケールしている。これは本当に面白いと思う。
スポンサーサイト

ベンチャーブログのランキングに参加しています。
下のバナーをクリックして応援していただけると嬉しいです。
にほんブログ村 ベンチャーブログへ

GaucheとPHPの速度比較

Scheme/Gauche
2008.12.10
GahceとPHPの処理速度比較をやってみた。Gaucheのスタディを始めてからずっと気になっていたのが処理速度だったからだ。いくらコーディングしやすくても、処理時間がやたらかかるようでは実用的なアプリケーションなんて作れない。場合によってはC言語も併用しながら、チューニングを行うはめになってしまう。とはいえ、まあ半分は興味本位なんだけどね。

とりあえずやらせてみたのは、フィボナッチ数列の計算だ。階乗計算でもいいかなあと思ったのだが、PHPで階乗計算をやってみると、ある程度以上の桁数になった時点で答えが「INF」になってしまうのだ。しかも「INF」にならないようにすると、あまり処理時間が稼げないので比較が難しい。というわけで、答えの桁数があまり大きくなく、しかもやたらと処理に時間がかかるものをということで、定番のフィボナッチ数列を選んだわけである。

ご存じのとおり、フィボナッチ数列とは、以下のようなものである。

fib(1) = 1
fib(2) = 1
fib(3) = 2
...
fib(n) = fib(n-1) + fib(n-2)


処理時間を稼ぎたいので、処理方式は再帰型とする。これならnの二乗のペースで処理プロセスが大きくなっていくからだ。

Gaucheのコードは以下の通り。n=30の値を計算している。

(define (fib n)
  (cond ((= n 1) 1)
     ((= n 2) 1)
     (else (+ (fib (- n 1)) (fib (- n 2))))))

(time (fib 30))



PHPのコードは以下の通り。計算内容は上と同じだ。なお now_time()は現在時刻をマイクロ秒まで取得するためのもの。単位は秒である。

function now_time (){
  list ($msec, $sec) = explode(' ', microtime());
  return( (float)$msec + (float)$sec );
}

function fib( $n ){
  if( $n == 1 ){
    return( 1 );
  }else if( $n == 2){
    return( 1 );
  }else{
    return( fib( $n-1 ) + fib( $n-2 ) );
  }
}

$n = 30;
$start = now_time();
echo "Ans=" . fib( $n ) . "<br />";
$end = now_time();
echo "Time=" . ($end - $start);



使用したマシンは、昨年の秋に中古屋で買った日立のフローラ。プロセッサはPentium 4 の2.4GHz、今年の春にCentOSを導入して、そのまま使っているものだ。

計測結果は以下の通り。

Gauche
real:0.226 user:0.230 sys:0.000

PHP
Time=2.026


たったひとつのテストで結論を出すべきではないが、この結果を見る限り、Gaucheの方がPHPよりも約9倍速いことになる。ちなみに n=35 にすると、Gaucheでは 2.47秒、PHPでは22.86秒かかる。PHPからGaucheへの移行は、パフォーマンス面のメリットもあるのかもしれない。もちろんGCIによる呼び出しやDBアクセスのオーバーヘッドも考える必要があるので、内部処理だけで比較するのは問題があるのだが。

でも個人的には、いくぶんホッとしたのである。

ベンチャーブログのランキングに参加しています。
下のバナーをクリックして応援していただけると嬉しいです。
にほんブログ村 ベンチャーブログへ

Scheme/GaucheスタディのためのEmacs設定

Scheme/Gauche
2008.12.09
Scheme/Gaucheをスタディする環境として、LinuxとCywin上でEmacsを使っている。メインはLinux上のEmacsなんだけど、モバイル環境ではノートPCにインストールしたCygwinも使えるようにしているのだ。ここで重要なのが.emacsの設定である。基本的には「Think IT」に掲載されていた吉田 裕美さんの「Gaucheを始めてみませんか?」にあった設定を、ほぼそのままパクっているのだが、そのままではどうもしっくりこない。そこで若干カスタマイズして使っている。

ここでは自分の備忘録の意味も含めて、現在の.emacsの内容を掲載しておきたいと思う。

;;; http://www.thinkit.co.jp/article/74/1/3.html から引用。
(setq process-coding-system-alist
(cons '("gosh" utf-8 . utf-8)
process-coding-system-alist))
(setq scheme-program-name
"/usr/local/bin/gosh -i")
(autoload 'scheme-mode "cmuscheme"
"Major mode for Scheme." t)
(autoload 'run-scheme "cmuscheme"
"Run an inferior Scheme process." t)
(defun scheme-other-window ()
"Run scheme on other window"
(interactive)
(switch-to-buffer-other-window
(get-buffer-create "*scheme*"))
(run-scheme scheme-program-name))
(define-key global-map
"\C-cs" 'scheme-other-window)
(show-paren-mode)
;;; 引用はここまで。

;;; ここからは独自の設定。

;; バックスペースキーの定義。
;; やっぱりバックスペースはBSキーでやりたいから。
(define-key global-map
"\C-h" 'delete-backward-char)

;; 画面幅をはみ出した場合、次の行に表示する。
(setq truncate-lines nil)
(setq truncate-partial-width-windows nil)

;; 自動的に画面を縦に分割し、左側でgoshを起動。
;; いちいちC-x3、C-cs、C-xoをタイプするのが面倒だから。
(split-window-horizontally)
(other-window 1)
(scheme-other-window)

;; 制御を右側(エディタ)画面に移しておく。
(other-window 1)

;; エディタ側の画面を少し大きく。
;; 幅140文字のターミナルで、goshを60文字、
;; エディタを80文字の幅にしておく。
(enlarge-window-horizontally 10)


独自設定の部分は、他の人の参考にもなると思う。私自身、こういう設定を行っている人がいてもおかしくないと思うのだが(だってEmacsを起動するたびにC-x3 C-cs C-xoをタイプするなんて面倒でしょう?)、ググった範囲では見つけられなかった。もしよろしければご参考までに。

なおEmacsの設定に関するページはけっこう存在するが、やはり「GNU Emacs Lispリファレンスマニュアル」が最も情報量が多く、使い物になるということもわかった。こちらもご参考までに。

ベンチャーブログのランキングに参加しています。
下のバナーをクリックして応援していただけると嬉しいです。
にほんブログ村 ベンチャーブログへ

向かい風とギア付き自転車の速度

Scheme/Gauche
2008.12.07
Newton法が簡単に解けたので、「向かい風はなぜ強く感じるか」で導き出した方程式を解いてみたいと思う。

この方程式は、連続的に変化するギヤを持つ自転車を想定し、向かい風の風速を w 、無風状態の巡航速度を v0 、ケイデンス一定、かける力一定の場合、自転車の速度 v がどうなるのかを考えたものだ。この方程式が正しいかどうかはよくわからないが、とりあえず解いてみよう。方程式は以下の通りである。

 f(v) = v*v*v + 2*w*v*v + w*w*v - v0*v0*v0 = 0


これを前回作成したnewton関数で解く。ついでにリファクタリングも行っておく。前回のnewton関数は、再帰ループに入るたびに微分関数のdefineを評価していたので、このdefine部分と再帰ループ部分を分離するのだ。以下のようになる。

 (define (newton func a dx dy)
  (define (deriv x) (/ (- (func (+ x dx)) (func x)) dx))
  (define (newton-rec a)
   (cond ((< (func a) dy) a)
    (else (newton-rec (- a (/ (func a) (deriv a)))))))
  (newton-rec a))


パラメーターとしてwとv0を与えて解く必要があるので、これらのパラメーターを受け取ってnewton関数を実行するspeed関数を定義する。

 (define (speed w v0)
  (newton (lambda(x)
       (+ (* x x x) (* 2 w x x) (* w w x) (* -1 v0 v0 v0)))
       1000 0.001 0.001))


実数解がひとつしかない3次方程式の場合、初期値が中途半端だと無限ループに入る可能性があるため、ある程度大きな値を初期値にしておく。計算精度はそれほど必要ないので、dx、dy共に0.001程度にした。

このspeed関数を、wの値を変えながら計算し、wとvの関係を求めていく。またv0も複数のパターンを用意する。たとえばwを0~30の間で2ずつ変化させて、v0=15とv0=30のvの値を求めるには、次のようなプログラムを実行すればいい。

 (use srfi-1)

 (map (lambda(w)
  (print "w=" w
      " v(15)=" (speed w 15)
      " v(30)=" (speed w 30)))
  (iota 16 0 2))


上記プログラムをちょっと手直しし、v0=5、15、25、35、45、w=-30~+30で、vがどのようになるのかをグラフ化してみた。

グラフ

方程式が正しいかどうかはともかくとして、手軽に解ける手段を持つのは悪くない。グラフを見ると、v0が変化してもグラフの傾きはそれほど変わらない。これは無風状態の巡航速度との比率で考えれば、v0が大きいほど風の影響が小さく感じられることを意味する。

もちろんwではなくw/v0の影響という観点で見れば影響はすべての場合でほぼ同じになる。しかしw/v0ではなくwの影響を考えた方が実感には合うだろう。つまり同じ風が吹いていた場合、無風状態の巡航速度が速いライダーほど、風の影響を感じなくなるのだ。

私自身の経験から言えば、この結論は妥当性があると思う。以前に比べて巡航速度が上がった結果、あまり向かい風が気にならなくなっているからである。

みなさんはどう思われるだろうか。

ベンチャーブログのランキングに参加しています。
下のバナーをクリックして応援していただけると嬉しいです。
にほんブログ村 ベンチャーブログへ

Newton法で方程式を解く

Scheme/Gauche
2008.12.06
新しい言語のスタディを始めると、やはり簡単なプログラムを作って、これまで使っていた言語との比較を行いたくなる。前回ちょっと言及した東京大学理学部の「Scheme演習」の中に「Newton法で方程式を解け」という問題があったので、やってみることにした。

解を求めたい方程式を

 f(x)=0


とすると、Newton法による方程式の解は以下のように求められる。

Step1.
まず任意のxにおけるf(x)の微分値を求める関数を作成する。この関数を deriv(x) とする。deriv(x)は以下のように定義できる。

 deriv(x) = ( f(x+dx) - f(x) ) / dx



dxが十分に小さい値であれば、deriv(x)の値は十分に微分値に近づくと考えられる。

Step2.
適当なxの値(これを a とする)におけるf(x)の値を求める。つまりf(a)を求めるわけである。

Step3.
x=a における deriv(x) 、つまり deriv(a) を求める。

Step4.
( a - f(a) / deriv(a) ) を求め、これを新しい a とする。

Step5.
Step2~4を繰り返し、f(a)の値を0に近づける。0との差が十分に小さい値になったときに(これを dy とする)、その時の a の値を方程式の解とする。

方程式によっては、最初のaの値によって求められる解が変わってしまったり、解が得られないこともあるが、とりあえず上記の方法で解を求めてみよう。

Scheme/Gaucheでプログラムを作成した場合、コードは次のようになる。

 (define (newton func a dx dy)
  (define (deriv x)
    (/ (- (func (+ x dx))(func x)) dx))
  (cond ((< (func a) dy) a)
    (else (newton func (- a (/ (func a) (deriv a))) dx dy))))


ここで以下のようにf(x)を定義し

 (define (f x)(- (* x x) 2))



以下のようにnewtonを呼び出せば、

 (newton f 5 0.0000000001 0.0000000001)



方程式のひとつの解が得られる。

 gosh> 1.4142135623730887



さてこれを、PHPで書いたらどうなるだろう?
再帰処理を使う場合には、次のようになるはずだ。

 function deriv($f,$x,$dx){
  return((call_user_func($f,$x+$dx)-call_user_func($f,$x))/$dx);
 }

 function newton($f,$a,$dx,$dy ){
  if(call_user_func($f,$a)<$dy){
   return($a);
  }else{
   return(newton($f,$a-(call_user_func($f,$a)/deriv($f,$a,$dx)),$dx,$dy));
  }
 }


さらに f(x) を定義し、newtonを呼び出せばいい。

 function f($x){
  return($x*$x-2);
 }

 echo newton('f',5,0.0000000001,0.0000000001);


Webブラウザからこれを実行すれば、以下のように表示される。

 1.4142135623731



Scheme/Gaucheのように関数内で関数を定義することもできるが、再帰処理を行う場合には内部で定義した関数が重複定義となるため、実行時にエラーになる。そのため微分値を求める deriv 関数は newton 関数の外で定義している。

再帰処理をループに置き換えれば、deriv 関数を newton 関数の中に定義することもできる。

 function newton( $f, $a, $dx, $dy ){
  function deriv( $f, $x, $dx){
   return( ( call_user_func( $f, $x + $dx ) - call_user_func( $f, $x ) ) / $dx );
  }
  while( call_user_func( $f, $a ) >= $dy ){
   $a = $a - ( call_user_func( $f, $a ) / deriv( $f, $a, $dx ) );
  }
  return( $a );
 }


使用法は同じである。

さて、この程度のプログラムであれば、どちらでも簡単に実現できる。ここで私が注目したいのは、どちらの方がシンプルで美しく見えるか、ということである。個人的にはScheme/Gaucheの方がまとまりがいいような気がする。例えばPHPの call_user_func は、ちょっとしたハックの匂いがする。関数名ではなくて、関数そのものを渡せればいいのになあ、と思っちゃうわけだ。

まあこれは、好みの問題なのかもしれませんけどね。

ベンチャーブログのランキングに参加しています。
下のバナーをクリックして応援していただけると嬉しいです。
にほんブログ村 ベンチャーブログへ

FC2Ad

相続 会社設立

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。