6

すべてのモツキン数を生成して配列に格納したいと思います。式は次のとおりです。 ここに画像の説明を入力してください

私の現在の実装は遅すぎます:

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

その上、べき乗二乗を適用できるように、漸化式行列の閉じた形を見つけることに固執しています。誰かがより良いアルゴリズムを提案できますか?ありがとう。
編集数値を法とする場合は除算が適用されないため、2番目の式を適用できません。最大値nは10,000で、64ビット整数の範囲を超えているため、答えはモジュロであり、数値mm= 10 ^ 14+7です。これより大きな整数ライブラリは許可されていません。

4

2 に答える 2

1

確かに、2 番目の式を使用できます。除算はモジュラー乗法逆数で行うことができます。モジュラー数が素数でなくても可能ですが、それも可能です ( MAXGAME チャレンジの議論でいくつかの役立つヒントを見つけました)。

MOD を 43 * 1103 * 2083 * 1012201 として素因数分解します。これらの各素数を法としてすべての量を計算し、中国剰余定理を使用して MOD を法とする値を見つけます。ここでは除算も含まれているので注意してください。それぞれの量について、それらを分割するこれらの素数のそれぞれの最高の累乗を維持する必要があります。

次の C++ プログラムは、100000000000007 を法とする最初の 10000 個のモツキン数を出力します。

#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;
}
于 2012-09-01T22:21:10.293 に答える
0

警告: 次のコードは、整数除算を使用しているため間違っています (例: 2.5 ではなく 5/2 = 2)。お気軽に修理してください!

これは、動的計画法を使用する良い機会です。フィボナッチ数を計算するのと非常によく似ています。

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
"""

問題は、これらの数値が非常に大きくなるため、すべてをキャッシュに格納すると、非常に高くしたい場合にメモリが不足することです。次に、前の 2 つの用語を覚えている for ループを使用することをお勧めします。多くの数値のモツキン数を見つけたい場合は、最初に数値を並べ替えてから、for ループで各数値に近づくにつれて結果を出力することをお勧めします。

編集: ループ バージョンを作成しましたが、以前の再帰関数とは異なる結果が得られました。そのうちの少なくとも 1 つが間違っている必要があります。うまくいけば、あなたはまだそれがどのように機能するかを見て、それを修正することができます!

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)
"""
于 2012-08-10T00:10:51.673 に答える