Common Lisp/Advanced topics/Numbers/Example 1

Problem: given a function f defined on complex numbers and a square area that contains only one root of the function, find this root.

We will use residue theorem for the function 1/f(x). First we need to be able to calculate integrals on a square path.

(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))))))))

The main loop could be made more concise with the use of extended loop syntax. The above function is very procedural in nature. You would use the same algorithm in C-like programming languages. Still, Lisp has a little advantage due to its native complex number support and the fact that case returns value, unlike switch in C.

Note the use of float function that converts the result of division to float. Without it, Lisp would operate with rational numbers, and the result won't be pretty (unless you find rationals with 1000-digit denominators pretty). As soon as the function is loaded into Lisp, you can test it:

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

This correlates with the theory that predicts 2πi being the result.

Now, the corollary of the residue theorem states that there is a pole in the area iff the path integral that goes around that area is not zero. We can write a simple function on top of the previous one that provides what we need:

(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))))

The return value would be useful in the recursion-terminating case of the next function, which divides a square into 4 smaller squares and uses indirect recursion to complete its task.

            
(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)))))))

Now, this is an example how functional programming can save lines of code. We define a small helper function that uses its lexical environment to know the values of f, start, h and precision, so the only thing we need to pass to it is the upper-right corner of the square. The power of or macro saved some superfluous branching as well. This function, while elegant, is quite hard to understand. It's a good exercise to work out what do check-pole and find-pole return in different situations and how their return values affect the control flow.

Finally, to find a root of function f we need to find a pole of 1/f. This is quite easy to do.

(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))

Let's test it. f(x)=x²+2 has two complex roots: ±sqrt(2)*i. Let's see if our program can find the upper one:

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

Looks like the right answer!