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 divi

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

StackExchange Computer Science Q#3109, answer score: 4