生成第 n 个 Motzkin 数的最快方法是什么？

Question

我想生成所有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。不允许使用更大的整数库。

Answer 1

确实可以使用第二个公式。可以使用modular multiplicative inverse 进行除法。即使模数不是素数，这很困难，也有可能（我在discussion 到MAXGAME challenge 中找到了一些有用的提示）：

将 MOD 因式分解为 = 43 * 1103 * 2083 * 1012201。以每个素数为模计算所有量，然后使用中国剩余定理找出模 MOD 的值。请注意，这里也涉及除法，因为每个量都需要保持这些素数中的每一个的最高幂，这也将它们除法。

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

Answer 2

警告：以下代码错误，因为它使用整数除法（例如 5/2 = 2 而不是 2.5）。随意修复它！

这是使用动态规划的好机会。计算斐波那契数非常相似。

sample code:

cache = {}
cache[0] = 1
cache[1] = 1

def motzkin(n):
    if n in cache:
        return cache[n]
    else:
        result = 3*n*motzkin(n - 2)/(n + 3) + (2*n + 3)*motzkin(n - 1)/(n + 3)
        cache[n] = result
        return result

for i in range(10):
    print i, motzkin(i)

print motzkin(1000)

"""
0 1
1 1
2 2
3 4
4 9
5 21
6 53
7 134
8 346
9 906
75794754010998216916857635442484411813743978100571902845098110153309261636322340168650370511949389501344124924484495394937913240955817164730133355584393471371445661970273727286877336588424618403572614523888534965515707096904677209192772199599003176027572021460794460755760991100028703368873821893050902166740481987827822643139384161298315488092901472934255559058881743019252022468893544043541453423967661847226330177828070589283132360685783010085347614855435535263090005810
"""

问题是因为这些数字变得如此之大，如果你想把它们全部存储在缓存中会耗尽内存。那么最好使用 for 循环记住前两个术语。如果您想找到许多数字的 motzkin 数，我建议您先对数字进行排序，然后在 for 循环中处理每个数字时，输出结果。

编辑：我创建了一个循环版本，但得到的结果与我之前的递归函数不同。至少其中一个肯定是错的！！希望您仍然可以看到它是如何工作的并且可以修复它！

def motzkin2(numbers):
    numbers.sort() #assumes no duplicates
    up_to = 0
    if numbers[0] == 0:
        yield 1
        up_to += 1
    if 1 in numbers[:2]:
        yield 1
        up_to += 1

    max_ = numbers[-1]
    m0 = 1
    m1 = 1
    for n in range(3, max_ + 1):
        m2 = 3*n*m0/(n + 3) + (2*n + 3)*m1/(n + 3)
        if n == numbers[up_to]:
            yield n, m2
            up_to += 1
        m0, m1 = m1, m2



for pair in motzkin2([9,1,3,7, 1000]):
    print pair

"""
1
(3, 2)
(7, 57)
(9, 387)
(1000, 32369017020536373226194869003219167142048874154652342993932240158930603189131202414912032918968097703139535871364048699365879645336396657663119183721377260183677704306107525149452521761041198342393710275721776790421499235867633215952014201548763282500175566539955302783908853370899176492629575848442244003609595110883079129592139070998456707801580368040581283599846781393163004323074215163246295343379138928050636671035367010921338262011084674447731713736715411737862658025L)
"""