«前の日記(2008年10月28日) 最新 次の日記(2008年10月30日)» 編集

ema log


2008年10月29日 [長年日記]

_ [Gauche]SICP 1.29の説明

シンプソンの公式」によって,x^3 を 0〜1 の区間で数値積分するプログラムを作れというもの.一章の練習問題ってのが恐ろしい.実装は単純だし,行数は短いのだけど,Cだと大騒ぎか.「シンプソンの公式」がなにかはWikipediaを読んでもらうとして,

0 が 0 の時のエラー処理はどうなんだ?とかはさておき.

という数式をプログラム化すればいい.まずは,「Σ」に相当する関数 sum.

(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a)
         (sum term (next a) next b))))

これは,本文中から得られる.term, next に関数を渡せるのが Scheme の良いところ.

(cube 0) + (cube 0.2) + (cube 0.4) + (cube 0.6) + (cube 0.8) + (cube 1.0) + 0

みたいに,演算される様を想像できればOK.というか,Σを思い出せればOK.

後は,Wikipedia の数式の定義通りに実装している.y, h の定義は,Wikipedia の本文にある(Wikipediaではy_kではなくx_i).なお let は,その関数スコープでのみ参照できる変数を作る(実は,この演習問題時点では let が使えないらしいのですがw)

なお,非Gaucheユーザは1/3 は 1÷3 ではなく,1/3 という有理数リテラルであることに注意.

() のネスト具合がわかりにくいのは,インデントで判断するしかない.正直,まともなエディタの補助がないとかけない.

(define (cube x) (* x x x))

(define (simpson term a b n)
  (let ((h (/ (- b a) n)))
    (define (y k) (term (+ a (* k h))))
    (define (inc2 n) (+ n 2))
    (* h 1/3 (+ (y 0)
                (* 4 (sum y 1 inc2 (- n 1)))
                (* 2 (sum y 2 inc2 (- n 1)))
                (y n)))))

一番下の部分が,数式ほぼそのまま.Schemeすげえ!と思う瞬間.Simpson 自体が再帰ではないあたりが面白いというか,スマートというか.ただ,()の数もすごいw

C はおろか,Ruby ではここまでスマートにかけない.結局関数がファーストクラスオブジェクトじゃないからなぁ.同じコンセプトを実装できはするんだけど.

Ruby でまねてみると(next は succ に置き換えている.というのも,識別子のため,使用できなかった).n=10ではなぜか精度が出ず,n=1000で,再帰が深すぎると怒られた.読み取りやすくはあるかな.proc の呼び出しが [] なのがきもいなー.Gaucheの有理数を尊重してRational使うと

def sum(term, a, succ, b)
  if (a > b)
    0
  else
    term[a] + sum(term, succ[a], succ, b)
  end
end

cube = proc { |x| x*x*x }

def simpson(term, a, b, n)
  h = Rational((b - a), n)
  y = proc { |k| a + k * h }
  inc2 = proc { |n| n + 2 }
  h/3 * (y[0] +
         4 * sum(y, 1, inc2, n-1) +
         2 * sum(y, 2, inc2, n-1) +
         y[n])
end

puts simpson(cube, 0, 1, 100)

JavaScriptで似せると,以下の通り.return を素で忘れます.有理数がないので,普通に小数点演算.ちなみに,JavaScriptでは,1/2 は 0.5 になる.

function sum(term, a, next, b) {
  if (a > b) {
    return 0;
  } else {
    return term(a) + sum(term, next(a), next, b);
  }
}

var cube = function (x) { return x*x*x; };

function simpson(term, a, b, n) {
  var h = (b - a) / n;
  var y = function(k) { return a + k * h; };
  var inc2 = function(n) { return n + 2; };

  return h / 3 * (y(0) +
                  4 * sum(y, 1, inc2, n-1) +
                  2 * sum(y, 2, inc2, n-1) +
                  y(n));
}

alert(simpson(cube, 0, 1, 100));

sum は,Ruby だと後置ifの例外処理と通常処理

def sum(term, a, next, b)
  return 0 if (a > b)
  return term[a] + sum(term, next[a], next, b)
end

JavaScriptだと 三項演算子使った方が,好きかなぁとか思ったり.

function sum(term, a, next, b) {
  return (a > b) ? 0
                 : term(a) + sum(term, next(a), next, b);
}

やっぱり,シンプルにかけるかどうかは考え方次第だなー.

以上,気分転換.現実逃避.