Common Lisp/高階主題/數字/示例 1

問題：給定定義在複數上的函式f，以及只包含函式的一個根的正方形區域，找到這個根。

我們將使用 留數定理 用於函式 1/f(x)。首先我們需要能夠計算正方形路徑上的積分。

(defun integrate-square-path (f start length precision)
  "f is the function to integrate.
   Start is a complex number: the lower left corner of the square.
   Length is the length of the side of square.
   Precision is the distance between two points that we consider acceptable."
  (let* ((sum 0) ;;The result would be summed there
         (n (ceiling (/ length precision))) ;;How many points on each side
         (step (float (/ length n))) ;;Distance between points
         (j 0) ;;index
         (side 0) ;;The number of side: from 0 to 3
         (d (complex step 0)) ;;Complex difference between two points
         (cur start)) ;;Current position 
    (loop (incf sum (* (funcall f cur) d)) ;;Increment the sum
          (incf cur d) ;;Change the position
          (incf j) ;;Increment the index
          (when (= j n) ;;Time to change the side
            (setf j 0)  
            (incf side)
            (setf d (case side  ;;Change the direction
                      (1 (complex 0 step))
                      (2 (complex (- step) 0))
                      (3 (complex 0 (- step)))
                      (4 (return sum))))))))

主迴圈可以使用擴充套件的loop語法使其更簡潔。上面的函式本質上是過程式的。您將在類似 C 的程式語言中使用相同的演算法。但是，Lisp 由於其原生複數支援以及case返回值的事實（不像 C 中的switch），因此具有一定的優勢。

請注意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 相吻合。

現在，留數定理的推論指出，當圍繞該區域的路徑積分不為零時，該區域存在極點。我們可以基於前面的函式編寫一個簡單的函式，它提供了我們需要的功能

(defun pole-p (f start length precision)
  "Tests, whether the given square contains a pole of f and
returns the center of the square if there is a pole after all"
  (when (> (abs (integrate-square-path f start length precision)) (sqrt precision))
    (+ start (/ (complex length length) 2))))

返回值將在下一個函式的遞迴終止情況下很有用，該函式將一個正方形分成 4 個更小的正方形並使用間接遞迴來完成其任務。

            
(defun find-pole (f start length precision)
  "Finds a pole of the function in the square area"
  (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)
  "Finds a root of the function in the square area"
  (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)

看起來是正確的答案！