**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!