Pretty much the most interesting blog on the Internet.— Prof. Steven Landsburg

Once you get past the title, and the subtitle, and the equations, and the foreign quotes, and the computer code, and the various hapax legomena, a solid 50% English content!—The Proprietor

Wednesday, August 26, 2015

Tribonacci Coda in Mathematica

One question one might have after going through the Haskell Tribonacci problem is why it should be necessary to go through all that complication with recursive matrix multiplication in order to calculate such a simple function. Surely, one might say, there must be a mathematical formula into which you just plug in \(n\) and out pops the corresponding Tribonacci number.

The answer is: Yes, there is, but, well ... Have a look at it:

$$\begin{multline} \left(-\frac{1}{264} \left(1+i \sqrt{3}\right) \sqrt[3]{26136-4488 \sqrt{33}}-\frac{\left(1-i \sqrt{3}\right) \sqrt[3]{99+17 \sqrt{33}}}{4\ 33^{2/3}}\right) \\ \left(\frac{1}{3}-\frac{1}{6} \left(1-i \sqrt{3}\right) \sqrt[3]{19-3 \sqrt{33}}-\frac{1}{6} \left(1+i \sqrt{3}\right) \sqrt[3]{19+3 \sqrt{33}}\right)^{n+1} + \\ \left(-\frac{1}{264} \left(1-i \sqrt{3}\right) \sqrt[3]{26136-4488 \sqrt{33}}-\frac{\left(1+i \sqrt{3}\right) \sqrt[3]{99+17 \sqrt{33}}}{4\ 33^{2/3}}\right) \\ \left(\frac{1}{3}-\frac{1}{6} \left(1+i \sqrt{3}\right) \sqrt[3]{19-3 \sqrt{33}}-\frac{1}{6} \left(1-i \sqrt{3}\right) \sqrt[3]{19+3 \sqrt{33}}\right)^{n+1} + \\ \left(\frac{1}{132} \sqrt[3]{26136-4488 \sqrt{33}}+\frac{1}{66} \sqrt[3]{33 \left(99+17 \sqrt{33}\right)}\right) \\ \left(\frac{1}{3}+\frac{1}{3} \sqrt[3]{19-3 \sqrt{33}}+\frac{1}{3} \sqrt[3]{19+3 \sqrt{33}}\right)^{n+1} \\ \end{multline}$$

That is quite a monster, but note:

  1. Even though the formula contains complex numbers, the result is evidently real because the first two terms are complex conjugates of each other and the third term purely real.
  2. The functional dependence on \(n\) is actually quite simple; it is just a sum of three constants, each taken to the \(n+1\)th power and multiplied by a different set of constants.
  3. The exponentiated terms are the eigenvalues of \(\hat{a}\) (i.e., the roots of its characteristic polynomial ). The approximate, absolute value of the first two roots is \(0.737353\). Hence, for sufficiently large \(n\), they become vanishingly small. For the third root, it is \(1.83929\) and it grows exponentially. So to use this formula on large \(n\), you only need calculate the third term and round to the nearest integer.
  4. The previous observation that, for large \(n\), the result appears to have \(0.264649 n\) decimal digits is explained by \(log_{10}1.83929=0.264649\).

If you only want to calculate, let's say, the magnitude and leading digits of a high Tribonacci number, this formula is likely your best bet. Just remember that your intrinsic floating-point functions will start to fail around \(n \approx 1000\) if you aren't careful. For the trailing digits, you'd of course use the matrix method, but mod down each matrix multiplication by \(\approx 10^{16}\) to prevent the elements from becoming unnecessarily large.

If anyone is curious how I came up with this formula: I used Mathematica (one of the few pieces of software I hold in genuine awe). While Mathematica can solve a shocking number of symbolic problems on its own, I mostly use it in what one might call Mathematician's Helper or Smart Blackboard mode. In that mode, it is quite superior to the pencil-and-paper method. One still specifies all the transformations manually, but Mathematica keeps track of all tedious details which prevents the—in my case, sadly all-too-common—forgotten or miscopied terms or conditions, and furnishes helpful numeric figures and visualization aids at the press of a button.

Here is the Mathematica code which helped me discover this equation:

a = {{1, 1, 1}, {1, 0, 0}, {0, 1, 0}}; es = Eigensystem[Transpose[a]]; d = DiagonalMatrix[es[[1]]]; v = es[[2]]; u = FullSimplify[Inverse[es[[2]]]]; (* A quick check that a is indeed the matrix product of u, d, and v *) FullSimplify[u.d.v] == a (* Another quick check that v and u are indeed inverses *) FullSimplify[v.u] == IdentityMatrix[3] (* The Mathematica equivalent of the Haskell function *) tribonacci3[n_] := MatrixPower[a, n][[1, 1]] (* The new formula based on the diagonalized matrix *) tribonacci4[n_] := FullSimplify[(u.DiagonalMatrix[Map[#^n &, es[[1]]]].v)[[1, 1]]] (* One last check that they are the same *) Table[tribonacci3[n], {n, 1, 20}] == Table[tribonacci4[n], {n, 1, 20}] (* Produce the formula in nice TeX format *) tribonacci4[n] // ToRadicals // TeXForm

Update: Edited to clean up the Mathematica code and note one more helpful property of the formula.