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)

看起来是正确的答案！