August 19, 2019

leonardo's recurrance

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!

Content by © Jared Davis 2019-2020

Powered by Hugo & Kiss.