# constructor for Z_phi numbers --- creates a "Z_phi" S3 object
Z_phi <- function(a, b=0) {
z <- list(a = a, b = b)
class(z) <- "Z_phi"
z
}
# render it nicely on the screen
format.Z_phi <- function(x, ...) {
sign <- if (x$b >= 0) '+' else '-'
paste0(x$a, sign, abs(x$b), 'φ')
}
print.Z_phi <- function(x, ...) {
cat(format(x), "\n")
}Computing big Fibonacci numbers using the Golden Ratio, abstract algebra, and S3 classes in R
Programming
R
Math
This post is a bit interdisciplinary. It stems from a back and forth that I had on X about ways of computing the Fibonacci numbers. I learned something new in that interaction, and I thought that it would be a wonderful opportunity to demonstrate some cool math and some cool R programming. If you already know the math, you might learn something about R. If you already know about R, you might learn something about math.
One of my favourite facts is that there exists a closed form expression for the \(n\)-th Fibonacci number.
\[F_n = \frac{\varphi^n - (1-\varphi)^n}{\sqrt{5}}\]
where \(\varphi = \frac{1+\sqrt{5}}{2}\) is the Golden Ratio. This formula is known as Binet’s formula. It was very surprising to me to learn that there should be a closed-form expression for the n-th Fibonacci number, let alone that this mess of very irrational-looking numbers can somehow sum to an integer.
If you want to use this in practice, you do have a bit of a problem: working with irrational numbers like \(\varphi\) and \(\sqrt{5}\) will inevitably lead to precision problems. But what I’ve recently learned, and what this post is about, is that you can use the algebraic properties of \(\varphi\) to avoid dealing with irrational numbers at all.
One way to see how the Fibonacci numbers connect to the Golden Ratio is by thinking about the equation that traditionally defines \(\varphi\): \(\varphi\) is the number which is one less than its square, that is, \(\varphi^2 = \varphi + 1\). Just by the basic properties of exponents, we can rewrite this same equation as \(\varphi^2 = \varphi^1 + \varphi^0\). This looks a bit Fibonacci-ish: the second power of \(\varphi\) is the sum of its first and zero-th powers.
What about \(\varphi^3\)? \[\begin{align} \varphi^3 = \varphi \times \varphi^2 = \varphi \times (\varphi^1 + \varphi^0) = \varphi^2 + \varphi^1 \end{align}\]
The third power of \(\varphi\) is the sum of its first and second powers. Very Fibonacci-ish again. In fact the same reasoning will show that for any integer \(n\), we will have \(\varphi^n = \varphi^{n-1} + \varphi^{n-2}\). So the exponents \(n\) on \(\varphi\) are subject to the Fibonacci recurrence in this way.
I’ve known all of this for a long time, but what I recently learned is that you can use all of this to compute Fibonacci numbers without ever computing the value of \(\varphi^n\) itself, or even dealing with any irrational numbers at all.
The way this is going to work is to consider numbers of the form \(a + b\varphi\), where \(a\) and \(b\) are integers. We’re actually going to consider a whole sequence of these numbers, defined by \(f_0 = 0 + \varphi\), and \(f_n=(f_0)^n = f_0 \times f_{n-1}\).
Here’s where I’m going to use a little bit of R. I’m going to create an S3 class to represent this type of number. I’ll call these Z_phi numbers. If you’ve never seen S3 classes before, a very good overview is given in Hadley Wickham’s Advanced R book. If you know even a little bit about object oriented programming from another language, it’ll look a bit weird, but you should be able to follow along. If you don’t know any programming at all, just go with it; I think this should all be actually quite readable.
Let’s try it out.
print(Z_phi(3, 4))3+4φ print(Z_phi(1, -1))1-1φ Neat. Next, we have to give them a way to add and multiply. Adding should be easy: for two of these Z_phi numbers \(z_1 = a_1 + b_1 \varphi\) and \(z_2 = a_2 + b_2 \varphi\), the coefficients should just sum straightforwardly: \(z_1 + z_2 = (a_1 + a_2) + (b_1 + b_2)\varphi\).
Multiplication requires just a bit of simple algebra to figure out. (The third line follows from the fact that \(\varphi^2 = \varphi + 1\)).
\[\begin{align} z_1 \times z_2 &= (a_1 + b_1 \varphi) \times (a_2 + b_2 \varphi)\\ &= a_1a_2 + a_1 b_2 \varphi + b_1 a_2 \varphi + b_1 b_2 \varphi^2 \\ &= a_1a_2 + (a_1 b_2 + b_1 a_2) \varphi + b_1 b_2 (\varphi + 1) \\ &= a_1 a_2 + b_1 b_2 + (a_1 b_2 + b_1 a_2 + b_1 b_2) \varphi \end{align} \]
Note that this is another Z_phi number since \(a_1 a_2 + b_1 b_2\) and \(a_1 b_2 + b_1 a_2 + b_1 b_2\) are both integers. This fact is actually quite important, and nice. It means that we can add and multiply any two Z_phi numbers and we get another Z_phi number. In math—specifically, abstract algebra—we say that this makes the set of Z_phi objects a ring, which is the name for a set of objects that you can add and multiply together. The name of this particular ring is \(\mathbb{Z}[\varphi]\), which is why I’ve given these objects this name in my code.
We can implement these operations in R as follows — if you’ve never seen this before, it’s going to look a bit crazy because R is weird sometimes, but this is really useful to know how to do.
"+.Z_phi" <- function(e1, e2) {
Z_phi(e1$a + e2$a, e1$b + e2$b)
}
"*.Z_phi" <- function(e1, e2) {
a1 <- e1$a
b1 <- e1$b
a2 <- e2$a
b2 <- e2$b
Z_phi(a1 * a2 + b1 * b2, a1 * b2 + b1 * a2 + b1 * b2)
}This lets us add and multiply Z_phi numbers just like any other numbers.
z1 <- Z_phi(1, 2)
z2 <- Z_phi(3, 4)
print(z1 + z2)4+6φ print(z1 * z2)11+18φ With this, we can execute the plan I laid out above. Recall, we were going to look at powers of \(f_1 = \varphi = 0 + 1\varphi\). Let’s look at the first 10.
f_0 <- f_i <- Z_phi(0, 1)
for (i in 1:10) {
cat(str_glue("φ^{i} = {format(f_i)}"), "\n")
f_i <- f_i * f_0
}φ^1 = 0+1φ
φ^2 = 1+1φ
φ^3 = 1+2φ
φ^4 = 2+3φ
φ^5 = 3+5φ
φ^6 = 5+8φ
φ^7 = 8+13φ
φ^8 = 13+21φ
φ^9 = 21+34φ
φ^10 = 34+55φ Hopefully, something is jumping out at you. Both the sequence of \(a\)’s and the sequence of \(b\)’s form the Fibonacci sequence! To be specific, we can say that the \(n\)-th Fibonacci number is the \(b\) coefficient of \(\varphi^n\), when it’s written as a Z_phi number. I won’t provide a complete proof for why this is true—there should be enough here to put together the missing pieces if you want to convince yourself.
To me this feels like a cheat code. Using Binet’s formula, we can compute \(F_n\) in one fell swoop, but it requires doing arithmetic with the irrational number \(\varphi\), which can potentially be a bit troublesome. But here, we have a way of computing \(F_n\) that is somehow in terms of \(\varphi\), but does not actually involve computing or doing arithmetic using the value of \(\varphi\). The only adding and multiplying that takes place is on integers. We compute \(F_n\) by taking \(\varphi\) to the power of \(n\), but we’re using the algebraic properties of these Z_phi numbers to somehow sidestep the irrationality entirely.
Where it really starts to look like a cheat code is when you realize that now you can use the properties of exponentiation to cut down the number of steps it takes to calculate \(F_n\). The classic way of calculating \(F_n\) is just by recursion. Here’s an example implementation.
fib_classic <- function(n, zero = 0L, one = 1L) {
if (n == 0) {
return(zero)
} else if (n == 1) {
return(one)
}
v <- c(zero, one)
while (n > 1) {
v <- c(v[2], sum(v))
n <- n - 1
}
v[2]
}
fib_classic(10)[1] 55This is pretty straightforward. (Ignore these weird zero and one arguments for now; I’ll talk about that later.)
Anyway, this will take \(n\) steps to calculate \(F_n\). Now, you might think that this exponentiation method will take \(n\) steps as well, since it involves multiplying \(\varphi\) by itself \(n\) times, but we can shorten this with exponentiation by squaring. The idea here is that, for example, \(\varphi^{10} = ((\varphi^2)^2)^2 \times \varphi^2\), which can be accomplished in only around \(\log_2(n)\) steps. Here’s an implementation.
"^.Z_phi" <- function(z, n) {
stopifnot(n >= 1)
if (n == 1) {
z
} else {
switch(n %% 2 + 1, (z * z)^(n %/% 2L), z * (z * z)^((n-1L) %/% 2L))
}
}
Z_phi(0, 1)^1034+55φ With all of this scaffolding in place, we can write a really concise and neat Fibonacci function. (Again, ignore the slightly unexpected phi argument for now; I’ll come back to it).
fib_phi <- function(n, phi = Z_phi(0, 1)) {
z <- phi^n
z$b
}
fib_phi(10)[1] 55I just think that’s so cool. We’re using the algebraic properties of the Golden Ratio to compute Fibonacci numbers in an apparently closed form, \(\varphi ^ n\), but we’re doing it without having to actually perform any arithmetic on the Golden Ratio itself.
Now, what does this buy us? Is it faster? Based on how I’ve described what’s going on, you’d probably expect that it is. I did. Let’s test it out by calculating the 50th Fibonacci number.
bench::mark(classic = fib_classic(50), phi = fib_phi(50))# A tibble: 2 × 6
expression min median `itr/sec` mem_alloc `gc/sec`
<bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
1 classic 11µs 12.3µs 77984. 0B 93.7
2 phi 28.6µs 31.9µs 30501. 9KB 73.4Huh. Even though the Z_phi approach uses fewer steps, it still takes more than twice as long on average to find \(F_{50}\). Why?
I’m not an expert on this stuff, but I believe the answer is overhead. The classic method takes more steps than the \(\varphi\) method, but it’s still only 50 integer additions. That’s going to be fast. On the other hand, the phi method involves all kinds of extra baggage: dealing with the S3 class, putting things in lists and pulling them back out, multiplication. 50 is just too small for the extra investment to pay off.
But we should be able to see some evidence that the phi method is better if we look at how the time changes for increasing values of \(n\).
df <- map(
seq(10, 200, by = 10),
\(x) {
bench::mark(classic = fib_classic(x), phi = fib_phi(x)) |>
mutate(n = x)
}
) |>
list_rbind()
df |>
transmute(method = as.character(expression), n, median) |>
ggplot(aes(x = n, y = median, color = method)) +
geom_point() +
labs(y = expression("Median time to compute " * F[n]), x = "n",
title = expression("Comparison of methods up to "*F[200])) +
ggthemes::theme_few() +
bench::scale_y_bench_time(NULL)
Sure enough, you can see that the phi method overtakes the classic method somewhere around \(F_{200}\). More importantly though, there’s a clear difference in the shape of these two curves: the classic method is a very obviously straight line, whereas the phi method looks potentially logarithmic, as the above reasoning suggests it should be.
But there is a big problem with what I’ve just done. The Fibonacci numbers grow really fast, and in fact both of these methods are producing numbers larger than R knows how to work with precisely. They already disagree on the value of \(F_{100}\), and in fact neither one finds the correct value of \(F_{100} = 354224848179261915075\).
print(fib_classic(100), digits = 20)[1] 354224848179261997056print(fib_phi(100), digits = 20)[1] 354224848179261865984This is where the extra arguments that I added to the functions come in. To work with numbers this big, we’ll need a package that can deal with arbitrary-precision integer arithmetic. I’ll use the {bignum} package for this. By specifying the zero and one arguments to fib_classic as arbitrary-precision integers, the result will be an arbitrary-precision integer as well, and we get the right answer.
library(bignum)
format(fib_classic(100, biginteger(0), biginteger(1)), notation = "dec")[1] "354224848179261915075"Now, here’s something quite neat. I don’t actually have to write a new Z_phi implementation. Everything I’ve done will work just fine if I supply bigintegers to the arguments of Z_phi.
Z_phi(biginteger(0), biginteger(1))^100218922995834555169026+354224848179261915075φ Thus I can just supply a biginteger-valued phi to the fib_phi function. Here’s a comparison of the speed of the two functions at computing \(F_{100}\).
big_zero <- biginteger(0)
big_one <- biginteger(1)
bench::mark(classic = fib_classic(100, big_zero, big_one),
phi = fib_phi(100, phi = Z_phi(big_zero, big_one)))# A tibble: 2 × 6
expression min median `itr/sec` mem_alloc `gc/sec`
<bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
1 classic 12.86ms 13.96ms 65.9 51.3KB 34.5
2 phi 4.01ms 4.75ms 194. 16.5KB 24.3Now we’re talking! The Z_phi method takes around a third of the time to calculate \(F_{100}\) compared to the classic method.
In fact, when we switch to using arbitrary precision integers, fib_phi catches up much faster than it did using base numeric variables—it overtakes somewhere around \(F_{20}\).
df <- map(
seq(10, 200, 10),
\(x) {
bench::mark(
classic = fib_classic(x, big_zero, big_one),
phi = fib_phi(x, phi = Z_phi(big_zero, big_one))
) |>
mutate(n = x)
}) |>
list_rbind()
df |>
transmute(method = as.character(expression), n, median) |>
ggplot(aes(x = n, y = median, color = method)) +
geom_point() +
labs(y = expression("Median time to compute " * F[n]), x = "n",
title = expression("Comparison of methods up to "*F[200])) +
ggthemes::theme_few() +
bench::scale_y_bench_time(base = NULL)
Once again I have to say, this feels like a cheat code. We’re somehow calculating the exact value of the n-th Fibonacci number by doing some Golden Ratio magic, without actually having to know the value Golden Ratio, in logarithmic time.
Is 200 that big? Here’s the number.
format(fib_phi(200, Z_phi(big_zero, big_one)), notation = 'dec')[1] "280571172992510140037611932413038677189525"The 200th Fibonacci number is pretty big, but the logarithmic curve suggests we can go way bigger without having to wait too long. Let’s try a really big one.
bench::mark(phi = fib_phi(10000, phi = Z_phi(big_zero, big_one)))# A tibble: 1 × 6
expression min median `itr/sec` mem_alloc `gc/sec`
<bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
1 phi 13.5ms 14.2ms 70.4 112KB 20.9And here it is (checked against this random Gist with the first 10,000 Fibonacci numbers, and matches).
cat(format(fib_phi(10000, phi = Z_phi(big_zero, big_one)), notation = 'dec'))33644764876431783266621612005107543310302148460680063906564769974680081442166662368155595513633734025582065332680836159373734790483865268263040892463056431887354544369559827491606602099884183933864652731300088830269235673613135117579297437854413752130520504347701602264758318906527890855154366159582987279682987510631200575428783453215515103870818298969791613127856265033195487140214287532698187962046936097879900350962302291026368131493195275630227837628441540360584402572114334961180023091208287046088923962328835461505776583271252546093591128203925285393434620904245248929403901706233888991085841065183173360437470737908552631764325733993712871937587746897479926305837065742830161637408969178426378624212835258112820516370298089332099905707920064367426202389783111470054074998459250360633560933883831923386783056136435351892133279732908133732642652633989763922723407882928177953580570993691049175470808931841056146322338217465637321248226383092103297701648054726243842374862411453093812206564914032751086643394517512161526545361333111314042436854805106765843493523836959653428071768775328348234345557366719731392746273629108210679280784718035329131176778924659089938635459327894523777674406192240337638674004021330343297496902028328145933418826817683893072003634795623117103101291953169794607632737589253530772552375943788434504067715555779056450443016640119462580972216729758615026968443146952034614932291105970676243268515992834709891284706740862008587135016260312071903172086094081298321581077282076353186624611278245537208532365305775956430072517744315051539600905168603220349163222640885248852433158051534849622434848299380905070483482449327453732624567755879089187190803662058009594743150052402532709746995318770724376825907419939632265984147498193609285223945039707165443156421328157688908058783183404917434556270520223564846495196112460268313970975069382648706613264507665074611512677522748621598642530711298441182622661057163515069260029861704945425047491378115154139941550671256271197133252763631939606902895650288268608362241082050562430701794976171121233066073310059947366875Cool stuff!