我想生成所有 Motzkin Number并存储在一个数组中。公式如下:

我目前的实现速度太慢了:

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;
    }
}

此外，我一直在寻找递归矩阵的封闭形式，以便我可以应用指数平方。谁能建议一个更好的算法？谢谢。
编辑 我无法应用第二个公式，因为除法在对数字取模时不适用。 n 的最大值为 10,000，超出了 64 位整数的范围，因此答案是对更大的数 m 取模，其中 m = 10^14 + 7。不允许使用更大的整数库。

最佳答案

确实可以使用第二个公式。可以使用 modular multiplicative inverse 完成除法.即使模数不是质数，这使得它变得困难，它也是可能的(我在 discussion 到 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.

以下 C++ 程序打印前 10000 个 Motzkin 数模 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;
}

关于c++ - 生成第 n 个 Motzkin 数的最快方法是什么？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/11893036/