What is the fastest way to generate the n-th Motzkin number?

Question

I want to generate all the Motzkin Number and store in an array. The formula is given as follows:

My current implementation is just way too slow:

void generate_slow() {
    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;
    }
}

void generate_slightly_faster() {
    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 l = 0, r = i - 2; l <= r; ++l, --r) {
            if (l != r) {
                result = (result + (2 * (mm[l] * mm[r]) % MODULO)) % MODULO;
            }
            else {
                result = (result + ((mm[l] * mm[r]) % MODULO)) % MODULO;
            }
        }
        mm[i] = result;
    }
}

Besides, I'm stuck with finding a closed form for the recurrence matrix so that I can apply exponentiation squaring. Can anyone suggest a better algorithm? Thanks.
Edit I can't apply the second formula because division doesn't apply when modulo a number. The maximum of n is 10,000 which is beyond the range of 64-bit integer, so the answer is modulo with a larger number m where m = 10^14 + 7. Larger integer library is not allowed.

Answer 1

Indeed you can use the second formula. The division can be done with the modular multiplicative inverse. Even if the modular number is not prime, which makes it difficult, it is also possible (I found some helpful hints in the discussion to the MAXGAME challenge):

Prime factorise MOD as = 43 * 1103 * 2083 * 1012201. Compute all quantities modulo each of these primes and then use chinese remainder theorem to find out the value modulo MOD. Beware, as here divison is also involved, for each quantity one would need to maintain the highest powers of each of these primes which divides them as well.

Following C++ program prints the first 10000 Motzkin numbers modulo 100000000000007:

#include <iostream>
#include <stdexcept>

// Exctended Euclidean algorithm: Takes a, b as input, and return a
// triple (g, x, y), such that ax + by = g = gcd(a, b)
// (http://en.wikibooks.org/wiki/Algorithm_Implementation/Mathematics/
// Extended_Euclidean_algorithm)
void egcd(int64_t a, int64_t b, int64_t& g, int64_t& x, int64_t& y) {
    if (!a) {
        g = b; x = 0; y = 1;
        return;
    }
    int64_t gtmp, xtmp, ytmp;
    egcd(b % a, a, gtmp, ytmp, xtmp);
    g = gtmp; x = xtmp - (b / a) * ytmp; y = ytmp;
}

// Modular Multiplicative Inverse
bool modinv(int64_t a, int64_t mod, int64_t& ainv) {
    int64_t g, x, y;
    egcd(a, mod, g, x, y);
    if (g != 1)
        return false;
    ainv = x % mod;
    if (ainv < 0)
        ainv += mod;
    return true;
}

// returns (a * b) % mod
// uses Russian Peasant multiplication
// (http://stackoverflow.com/a/12171020/237483)
int64_t mulmod(int64_t a, int64_t b, int64_t mod) {
    if (a < 0) a += mod;
    if (b < 0) b += mod;
    int64_t res = 0;
    while (a != 0) {
        if (a & 1) res = (res + b) % mod;
        a >>= 1;
        b = (b << 1) % mod;
    }
    return res;
}

// Takes M_n-2 (m0) and M_n-1 (m1) and returns n-th Motzkin number
// all numbers are modulo mod
int64_t motzkin(int64_t m0, int64_t m1, int n, int64_t mod) {
    int64_t tmp1 = ((2 * n + 3) * m1 + (3 * n * m0));
    int64_t tmp2 = n + 3;

    // return 0 if mod divides tmp1 because:
    // 1. mod is prime
    // 2. if gcd(tmp2, mod) != 1 --> no multiplicative inverse!
    // --> 3. tmp2 is a multiple from mod
    // 4. tmp2 divides tmp1 (Motzkin numbers aren't floating point numbers)
    // --> 5. mod divides tmp1
    // --> tmp1 % mod = 0
    // --> (tmp1 * tmp2^(-1)) % mod = 0
    if (!(tmp1 % mod))
        return 0;

    int64_t tmp3;
    if (!modinv(tmp2, mod, tmp3))
        throw std::runtime_error("No multiplicative inverse");
    return (tmp1 * tmp3) % mod;
}

int main() {
    const int64_t M    = 100000000000007;
    const int64_t MD[] = { 43, 1103, 2083, 1012201 }; // Primefactors
    const int64_t MX[] = { M/MD[0], M/MD[1], M/MD[2], M/MD[3] };
    int64_t e1[4];

    // Precalculate e1 for the Chinese remainder algo
    for (int i = 0; i < 4; i++) {
        int64_t g, x, y;
        egcd(MD[i], MX[i], g, x, y);
        e1[i] = MX[i] * y;
        if (e1[i] < 0)
            e1[i] += M;
    }

    int64_t m0[] = { 1, 1, 1, 1 };
    int64_t m1[] = { 1, 1, 1, 1 };
    for (int n = 1; n < 10000; n++) {

        // Motzkin number for each factor
        for (int i = 0; i < 4; i++) {
            int64_t tmp = motzkin(m0[i], m1[i], n, MD[i]);
            m0[i] = m1[i];
            m1[i] = tmp;
        }

        // Chinese remainder theorem
        int64_t res = 0;
        for (int i = 0; i < 4; i++) {
            res += mulmod(m1[i], e1[i], M);
            res %= M;
        }
        std::cout << res << std::endl;
    }

    return 0;
}