I have finished reading *An Imaginary Tale* by Paul Nahin. About three quarters of the way through the book, I started fiddling with an example of applied complex analysis: representing a particular family of recurrence relations as an equation in complex numbers. I learned quite a bit along the way, and I even wrote some toy code to try out the efficiency of the approach.

## Solving Recurrence Relations

The genearlized Fibonacci sequence is

$$u_n = pu_{n + 1} + qu_{n + 2}$$

The common Fibonacci sequence has parameters $p = q = 1$; this can be solved for any values of $q$ or $p$.

### Solve the Quadradic

Nahin gives an example with $p = 4$ and $q = -1$. Assume $u{n} = kz^{n}$, both $k$ and $z$ constants. The recurrence becomes

$$kz^{n+2} = 4kz^{n+1} - 8kz^{n}$$

Divide through by $kz^{n}$ to isolate $z^{2}$ on the left hand

$$z^{2} = 4z - 8$$

then solve the quadratic

$$
\begin{align}
z & = \frac{4 \pm \sqrt{16 - 32}}{2} \\\\\

z & = 2 \pm i2
\end{align}
$$

### Write Out Initial Conditions

Now for the recurrence. The initial conditions are $u_{0} = u_{1} = 1$:

$$
\begin{align}
n = 0 & \colon k_1 + k_2 = 1\\\\\

n = 1 & \colon k_1 (2 + 2i) + k_2 (2 - 2i) = 1
\end{align}
$$

### Solve Simultaneous Equations

Here, I diverge from Nahin’s demonstration. He uses the exponetial form - partly for editorial reasons, I assume, given he had already demonstrated that $e^{i \pi}$ is a cleaner way to represent complex numbers. In most cases, it is; yet this calculation works much more cleanly with the algebraic form.

Our goal here is to solve two simultaneous equations. Recall that $i * i = -1$, and you can see

From here, observing $k_1 = 1 - k_2$ and $k_2 = 1 - k_1$, we get

and

### Convert to Exponential Form

Now we need to come back to Nahin’s approach, as it will be necessary to move forward. To convert a complex number $a + bi$ to exponential form $re^{i \theta}$, take $r = \sqrt{a^2 + b^2}$ and $\theta = \tan^{-1}{b/a}$. So

and

Giving us

as a conjugate pair. We’ll place these back into the general formula to obtain

### Simplify using Euler’s Formula

Euler’s Formula states that $e^{i \theta} = \cos \theta + i \sin \theta$. If we set the messy parts in the above equation to $c = 2^{3n/2 - 2} \sqrt{5}$ and $ \theta = n \pi / 4 + \tan^-1(1/2)$, then we have

Then using the summation formulat for $\cos \theta + \phi$, we get

and trigonometric function tables show us that

and

Substitute these into our working equation

### Express Simplified form in Code

The above is fairly simple two write and moderately efficient to run in code.

```
from numpy import cos, sin, pi, power, rint
def leonardo(n):
return rint(
power(2, (n * 3 / 2) - 1) * (2 * cos(n * pi / 4) - sin(n * pi / 4))
)
sequence(11) # returns -98304.0
```

### Compare with Naive Recursive and Iterative Recursive Forms

Is all that thinking on the part of the programmer worth it? Let’s do a quick benchmark against a naive recursive implementation and an iterative approach. The following benchmark code is written in Go. I have checked the accuracy of each function up to $n = 32$, at which point each begins to break down for various reasons beyond the scope of this exercise.

Here’s the code:

```
package leonardo_benchmark
import (
"math"
"testing"
)
// NaiveRecursive runs the Leonardo algorithm as a direct transliteration of its statement u_n = 4*u_n1 - 8*u_n2 into code
func NaiveRecursive(n int) int {
if n < 2 {
return 1
}
solution := 4*NaiveRecursive(n-1) - 8*NaiveRecursive(n-2)
return solution
}
// Iterative runs the Leonardo algorithm as an iteratively expressed transliteration of its original statement
func Iterative(n int) int {
solution := 1
u2 := 1
u1 := 1
if n < 2 {
return solution
}
for i := 2; i <= n; i++ {
u2 = u1
u1 = solution
solution = 4*u1 - 8*u2
}
return solution
}
// QuarterPi equals one-fourth PI
const QuarterPi = math.Pi * 0.25
// Analytic runs the Leonardo algorithm as a direct transliteration of our simplified analytic expression into code
func Analytic(n int) int {
power := (float64(3*n) / 2) - 1
konst := math.Pow(2, power)
theta := float64(n) * QuarterPi
solution := konst * (2*math.Cos(theta) - math.Sin(theta))
return int(math.Round(solution))
}
func BenchmarkNaiveRecursive(b *testing.B) {
for n := 0; n < b.N; n++ {
NaiveRecursive(32)
}
}
func BenchmarkIterative(b *testing.B) {
for n := 0; n < b.N; n++ {
Iterative(32)
}
}
func BenchmarkAnalytic(b *testing.B) {
for n := 0; n < b.N; n++ {
Analytic(32)
}
}
```

and some results on my MacBook Air:

```
jared$ GOPATH=$GOPATH:$(pwd) go test -bench=.
goos: darwin
goarch: amd64
BenchmarkNaiveRecursive-4 100 14022483 ns/op
BenchmarkIterative-4 50000000 25.9 ns/op
BenchmarkAnalytic-4 20000000 74.8 ns/op
PASS
ok _/Users/jared/dev/numerical-algorithms 4.347s
```

That’s a bit disappointing. Let’s try removing some of the higher-level code from the math package and see if we can improve:

```
func Analytic(n int) int {
power := (float64(3*n)/2 - 1)
konst := 2 << uint(power)
theta := float64(n) * QuarterPi
solution := float64(konst) * (2*math.Cos(theta) - math.Sin(theta))
return int(math.Round(solution) / 2)
}
```

```
jared$ GOPATH=$GOPATH:$(pwd) go test -bench=.
goos: darwin
goarch: amd64
BenchmarkNaiveRecursive-4 100 14841424 ns/op
BenchmarkIterative-4 50000000 26.3 ns/op
BenchmarkAnalytic-4 30000000 35.1 ns/op
PASS
ok _/Users/jared/dev/numerical-algorithms 4.695s
```

Performing the $2^{3/2n - 1}$ operation via bit-shifting more than doubled the performance. It’s still not as successful as the iterative solution, however, and it loses quite a bit of readability.

## What We’ve Learned

If you started out reading this with high hopes, you should by now realize in retrospect that an iterative solution would always perform better than an analytic solution. An analytic solution, while powerful for human thinking and understanding, requires much more overhead in programming language support; analytic equations also suffer from the physical characteristics of binary computers.

On the other hand, this exercise teaches quite a bit:

- How to view the linkage between iterative and analytic procedures
- How to manipulate infinite sequences with bounded functions (eg. trigonometric functions)
- How to translate pure mathematics to popular programming languages
- How to use a small portion of the
`numpy`

library in Python and`math`

package in Go. - How to benchmark Go code
- (If you click “view source” on the mathematics equations) how to write LaTex

and of course, when I say, “what we’ve learned” I mean what I have learned; I wrote this post over the course of several days to make sure I had the most straightforward calculations and the most readable code I could write (without spending *too much* time). I hope you’ve enjoyed reading!