# How to compute linear recurrence using matrix with fraction coefficients?

What I'm trying to do is generate Motzkin numbers mod a large number $10^{14} + 7$ (not prime), and it needs to compute the $n$th Motzkin number as fast as possible. From Wikipedia, the formula for the $n$th Motzkin number is defined as following:

\qquad \displaystyle \begin{align} M_{n+1} &= M_n + \sum_{i=0}^{n-1} M_iM_{n-1-i} \\ &= \frac{2n+3}{n+3}M_n + \frac{3n}{n+3}M_{n-1} \end{align}

My initial approach is to use the second formula which is obviously faster, but the problem I ran into is the division since modular arithmetic rule doesn't apply.

void generate_motzkin_numbers() {
motzkin = 1;
motzkin = 1;
ull m0 = 1;
ull m1 = 1;
ull numerator;
ull denominator;
for (int i = 2; i <= MAX_NUMBERS; ++i) {
numerator = (((2*i + 1)*m1 + 3*(i - 1)*m0)) % MODULO;
denominator = (i + 2);
motzkin[i] = numerator/denominator;
m0 = m1;
m1 = motzkin[i];
}
}


Then I tried the second formula, but the running time is horribly slow because the summation:

void generate_motzkin_numbers_nested_recurrence() {
mm = 1;
mm = 1;
mm = 2;
mm = 4;
mm = 9;
ull result;
for (int i = 5; i <= MAX_NUMBERS; ++i) {
result = mm[i - 1];
for (int k = 0; k <= (i - 2); ++k) {
result = (result + ((mm[k] * mm[i - 2 - k]) % MODULO)) % MODULO;
}
mm[i] = result;
}
}


Next, I'm thinking of using matrix form which eventually can be speed up using exponentiation squaring technique, in other words $M_{n+1}$ can be computed as follows: $$M_{n+1} = \begin{bmatrix} \dfrac{2n+3}{n+3} & \dfrac{3n}{n+3} \\ 1 & 0\end{bmatrix}^n \cdot \begin{bmatrix} 1 \\ 1\end{bmatrix}$$ With exponentiation by squaring, this method running time is $O(\log(n))$ which I guess the fastest way possible, where MAX_NUMBERS = 10,000. Unfortunately, again the division with modular is killing me. After apply the modulo to the numerator, the division is no longer accurate. So my question is, is there another technique to compute this recurrence modulo a number? I'm think of a dynamic programming approach for the summation, but I still think it's not as fast as this method. Any ideas or suggestions would be greatly appreciated.

• Sorry if I misunderstood, but can't you just use a lookup table instead of calculating them? The OEIS sequence A001006 lists them. Will you need to represent values of $M_n$ for large $n$? Meaning the numbers grow very quickly, and you can only represent the few first ones with 32/64-bit integers.
– Juho
Aug 9, 2012 at 19:15
• @Juho: Lookup table won't work because the limit of the file size, I need to store 10,000 numbers which I don't think it's possible. It surely pass 64 bit integer, that's why the answer is modulo with $m$. Aug 9, 2012 at 19:18
• @Chan: That file will be smaller than 1MB, so what's the problem? (Using machine-word-size ints for numbers of this kind is a crazy idea, by the way.)
– Raphael
Aug 9, 2012 at 21:19
• @Raphael: The file size is only 50000 bytes, while this is just one part of the problem, unless I make all my variables 1-2 characters, but I guess that's not the purpose of problem. By the way, the matrix formula given above is incorrect. I'm trying to find another way. Aug 9, 2012 at 21:43
• @Chan: By the way, there are several non-recursive formulae on the OEIS page. Alas, not a closed one as far as I can see, but you may be able to reduce your problem to one you can solve (quickly). Are asymptotics enough or do you need the precise numbers?
– Raphael
Aug 10, 2012 at 7:01

Your matrix formulation doesn't work- you need to use different $n$s in the matrices.
One idea to rescue your original approach is as follows. First, factor $10^{14} + 7 = 98794607 \times 1012201$ (only the second factor is prime). It is enough to compute $M_n$ modulo the two factors. The first factor is smaller than $2^{32}$, so you can store the first $10000$ Motzkin numbers modulo it using only 40KB. The second factor is prime and larger than $10003$, so you can use modular division, in which you multiply with the inverse (see here for how to calculate the inverse). This only works since $n+3$ will actually have an inverse, and that's why we factored $10^{14} + 7$.