Lisp/基本からさらに一歩進んで/数値/Example 1

問題: 複素数と、函数の唯一の根を含む正方形の面積によって定義された函数 f が与えられる。この場合の根(ルート、√)を求めよ。

この函数 1/f(x) のために留数定理を使います。まずは正方形上の経路積分をできるようになる必要があります。

(defun integrate-square-path (f start length precision)
  "f は積分のための函数。
   Start は複素数 : 正方形の下辺の左の角。
   Length は正方形の辺の長さ。
   Precision は許容可能な二点間の距離。"
  (let* ((sum 0) ;;この変数に結果が合計される。
         (n (ceiling (/ length precision))) ;;それぞれの辺にいくつのいくつの点があるか
         (step (float (/ length n))) ;;二点間の距離
         (j 0) ;;index
         (side 0) ;;辺の番号: 0 から 3 まで
         (d (complex step 0)) ;;2点間の複素数
         (cur start)) ;;現在位置
    (loop (incf sum (* (funcall f cur) d)) ;;合計を増加
          (incf cur d) ;;位置の変更
          (incf j) ;;indexを増加させる
          (when (= j n) ;;辺の変更
            (setf j 0)  
            (incf side)
            (setf d (case side  ;;方向の変更
                      (1 (complex 0 step))
                      (2 (complex (- step) 0))
                      (3 (complex 0 (- step)))
                      (4 (return sum))))))))

main loop は拡張された loop 構文でもっと簡潔にすることが出来ます。上記の関数は非常に手続き型として自然です。 C言語 のようなプログラミング言語でも同様のアルゴリズムを使うことができるでしょう。それでも、 Lisp は組み込みの複素数のサポートがあり、 C言語 の switch とは違い、返り値を返す case 文もあるために、ちょっとした利点があるのです。

float 関数は割り算の結果を小数型に変換するということに注意してください。こうしなければ、 Lisp は分数で計算を行い、結果が見栄えのいいものになりません。(1000桁の分母を持つ分数の見栄えがいいと思うのでなければですが。)この関数を Lisp で読み込むとすぐに試してみることが出来ます。

 CL-USER> (integrate-square-path (lambda (x) (/ 1 x)) #C(-1 -1) 2 0.001)
 #C(-0.0019999794 6.2832413)

これには、 2πi が結果になると予想する理屈が関連しています。

ここで、留数定理の推論によれば、 0 でない面積の線積分を行う限りにおいて、その面積には極がある、ということになります。ここで私たちが求める物を提供してくれる先ほど定義した関数の一つ上の階層のシンプルな関数を書くとこうなります。

(defun pole-p (f start length precision)
  "与えられた正方形が f の極を持つかどうか、ということと
もし極があるのなら、その中心を返す。"
  (when (> (abs (integrate-square-path f start length precision)) (sqrt precision))
    (+ start (/ (complex length length) 2))))

この返り値は、次の関数の再帰を終了する場合に役に立ちます。次の関数では正方形を四つの正方形に分け、この課題を終えるために相互再帰を使用します。

            
(defun find-pole (f start length precision)
  "正方形の面積にある函数の極を探し出す。"
  (when (> (* precision 10) length) 
    (return-from find-pole (pole-p f start length precision)))
  (let ((h (float (/ length 2))))
    (flet ((check-pole (start)
             (when (pole-p f start h precision)
               (find-pole f start h precision))))
      (or (check-pole start)
          (check-pole (+ start (complex 0 h)))
          (check-pole (+ start (complex h 0)))
          (check-pole (+ start (complex h h)))))))

これは、関数型プログラミングがいかにコードの行数を節約できるかという実例になります。ここでは小さな補助関数を定義します。これは、そのレキシカル環境のために f, start, h そして precision の値がわかっているので、この関数に渡すのは正方形の右上の角の座標だけでいいのです。また、 or マクロの機能が余分な分岐を省いてくれます。しかし、 check-pole と find-pole が違った状況で何を返すのか、返り値が制御構造にどのような影響を及ぼすのかを苦労して理解するのはいい練習になるでしょう。

最後に、函数 f の根を見つけるために、 1/f の極を見つけなければなりません。これは非常に簡単に出来ます。

(defun find-root (f start length precision)
  "正方形の面積内の函数の根を見つける"
  (find-pole (lambda (x) (/ 1 (funcall f x))) start length precision))

では試してみましょう。函数 f(x)=x²+2 は ±sqrt(2)*i という二つの複素数の根を持っています。ここまで定義したプログラムが上部の根を見つけられるか見てみましょう。

 CL-USER> (find-root (lambda (x) (+ (* x x) 2)) #C(-1 0) 5 0.0001) 
 #C(-5.493164E-4 1.4138794)

正しい答えを返しているみたいですね。