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