How to compute linear recurrence using matrix with fraction coefficients?

Question

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.

Answer 1

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$.