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

Posted 2023-02-22

【中文标题】生成第 n 个 Motzkin 数的最快方法是啥？

【英文标题】：What is the fastest way to generate the n-th Motzkin number?生成第 n 个 Motzkin 数的最快方法是什么？

【发布时间】：2012-08-09 22:40:06

【问题描述】：

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

【问题讨论】：

你可能想让你的标题更有趣一点；）

@therefromhere：谢谢，会的。

我不明白。您是否实现了仅依赖于 n、M_n 和 M_n-1 的 M_n+1 表达式？那应该很快。

@user827992：已经问过了，但我想这个问题与编程技术更相关。从第二个公式可以看出，数学公式是微不足道的。

你说“问题”——这是作业还是面试问题，还是一些编程挑战？

【参考方案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://***.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;

【讨论】：

我不确定它是否发生在所考虑的素数上，但是“如果 mod 除以 tmp1，则返回 0 的原因是：... (tmp1 * tmp2^(-1)) % mod = 0”是错误的。比如M_9 = 835不能被11整除，计算为M_9 = (19*323 + 24*127)/11 = (5*11*167)/11 = 5*167。

这里也有问题。 M_84 不能被 43 整除，(2*83*M_83 + 3*83*M_82) 的因式分解为[(2,1),(5,1),(19,1),(43,1),(15887,1),(42020639,1),(349259661165016424944007,1)]。

@DanielFischer：恐怕你是对的，我的推理是错误的。我将尝试更正该功能。感谢您指出我的错误。

【参考方案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)
"""

【讨论】：

这段代码对任何人有用还是我应该删除它？看到我忘记了模块化划分：X。您可以使用 Euclids 算法计算出该 mod 中的逆来修复它：X