確かに、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;
}
