# 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[0] = 1;
motzkin[1] = 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[0] = 1;
mm[1] = 1;
mm[2] = 2;
mm[3] = 4;
mm[4] = 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 '12 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$. – Chan Aug 9 '12 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 '12 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. – Chan Aug 9 '12 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 '12 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$.