6

The Tribonacci sequence satisfies

$$T(n) = T(n-1) + T(n-2) + T(n-3)$$

with $T(0)=0$, $T(1)=1$, $T(2)=1$. I need to calculate $T(y) \mod 10000$ for $y > 2^{40}$.

How can I make this faster? I know that this is periodic in $(\mathbb{Z}/10000\mathbb{Z})^3$, but I can't find the period.

Any suggestions? My program needs a lot of time to calculate such $T(y)$.

Chris Taylor
  • 28,955
florian
  • 63
  • 1
    Write the recursion in the form $v_{n+1} = A v_n$ where $v_n$ is the $3 \times 3$ vector with entries $T(n+2), T(n+1), T(n)$, then use binary exponentiation (http://en.wikipedia.org/wiki/Exponentiation_by_squaring) to evaluate $A$ to a large power. – Qiaochu Yuan Jun 16 '11 at 09:21
  • Also, the methods in http://math.stackexchange.com/questions/42880/why-are-the-periods-of-these-permutations-often-1560/42885#42885 can be used to find $d$ such that the period divides $d$, although more work is necessary because the modulus isn't squarefree... – Qiaochu Yuan Jun 16 '11 at 09:23
  • 2
    Converting the recurrence to matrix exponentiation is indeed the standard approach for these problems as Qiaochu says, but in this particular case the period is rather small relative to $10000^3$ — it's only $124000$ and can be found by computer in less than a second (by keeping track of the last three values, etc, and stopping when you reach (0,1,1) again). So you could also just use the period: $T(y) \equiv T(y\bmod 124000)$ and calculate the latter the normal way. Depends on the time constraints, of course. – ShreevatsaR Jun 16 '11 at 09:57
  • 2
    The characteristic polynomial is $p(x)=x^3-x^2-x-1$. This is irreducible mod 5, so we know that a root mod 5 is in the field $GF(125)$. Therefore the period mod 5 must be a factor of $124=4\cdot31$, and must be larger than 4. So only 31 needs to be checked, and turns out actually to be a period. Then we can proceed and compute the period mod 25, 125 ... In each step either the period stays the same or is multiplied by 5. Then we can do powers of 2 in the same way starting with the observation $p(x)(x+1)=x^4+1\pmod 2$, so mod 2 the period is 4. Chinese remainder theorem in the end. – Jyrki Lahtonen Jun 16 '11 at 11:03
  • My comment was just building on the stuff behind the link given in Qiaochu's comment. The familiar trick of computing the order modulo a prime power, when the order modulo a prime is known. – Jyrki Lahtonen Jun 16 '11 at 11:09
  • On repeated squaring of the shift matrix see also here and here (for Fib's, but the analogy for Trib's or any constant coef linear recursion is clear). – Bill Dubuque Jun 16 '11 at 15:46

2 Answers2

3

Call $U(n) = (T(n),T(n+1),T(n+2))$. The recurrence relation means that for all n, $U(n+1) = f(U(n))$ where $f$ is the linear transformation that sends $(a,b,c)$ to $(b,c,a+b+c)$.

Thus, in order to compute $T(n)$, instead of computing every $T(i)$ for every i, you can simply compute the linear transformation $f^n$, apply it to $U(0) = (0,1,1)$ to get $U(n) = (T(n),T(n+1),T(n+2))$.

To compute $f^{2^{40}}$ modulo 10000, write $f$ as a matrix with coefficients in $\mathbb{Z}/10000\mathbb{Z}$, and square it 40 times.

mercio
  • 50,180
  • 1
    He actually wants $T(y)$ for some $y > 2^{40}$ (not just $y = 2^{40}$), but otherwise this answer is fine. To compute $f^y$ we use exponentiation by squaring, as Qiaochu said. – ShreevatsaR Jun 16 '11 at 09:45
1

You're right that the sequence is periodic, and its period is less than $(10^4)^3$.

The following pseudocode will calculate the periodicity. The argument $m$ is the number that you are taking the modulus with respect to.

def TribonacciPeriod(m):
    a = 1; b = 1; c = 2 // manually do one iteration
    n = 1
    while (a != 0 or b != 1 or c != 1):
        tmp = (a + b + c) mod m
        a = b; b = c; c = tmp
        n += 1
    return n

This is guaranteed to terminate and return a value less than $m^3$.

Chris Taylor
  • 28,955
  • You also need to prove that the first time after 0 the sequence hits 0 is actually the period. Is that even true? – ShreevatsaR Jun 16 '11 at 09:34
  • No, this is wrong. Your code returns 991 and indeed $T(991) \equiv 0 \mod 10000$, but the next three numbers are $T(992) \equiv 3360$, $T(993) \equiv 1841$ and $T(994) \equiv 5201$. The actual period is much larger: it's $124000$. – ShreevatsaR Jun 16 '11 at 09:52
  • You're right, thanks. I edited the algorithm to make it correct (which has the unfortunate side effect that the comments no longer make sense - sorry about that.) – Chris Taylor Jun 16 '11 at 10:01
  • Now there's a theoretical problem with this approach (although it works for the specific initial values 0, 1, 1 and m=10000). The period could be as large as $m^3$ in the worst case, which for $m=10000$ is a bit too large for computer. :-) Unless one can prove that the period is always smaller… or, well, it happens to work for this particular problem so we don't have to think about what happens in general. – ShreevatsaR Jun 16 '11 at 10:06
  • 1
    Interestingly (in the sense that I wouldn't have guessed this) it appears that $\varphi(m)$, the worst-case run time to calculate the period of the Tribonacci sequence mod $k$ for $k<m$, scales like $\varphi(m)\sim m^2$, rather than $m^3$. Bonus points for an explanation of why this is :) – Chris Taylor Jun 16 '11 at 10:58
  • Hi, I did this and found a period all 124000 numbers. Could it be that F(2**40) mod 10000 is the same as F(31776) mod 10000? – florian Jun 16 '11 at 12:45
  • Yes, that is correct. – Chris Taylor Jun 16 '11 at 13:08