Indeed you can use the second formula. The division can be done with the modular multiplicative inverse. Even if the modular number is not prime, which makes it difficult, it is also possible (I found some helpful hints in the discussion to the 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.
Following C++ program prints the first 10000 Motzkin numbers modulo 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
// (https://mcmap.net/q/267374/-ways-to-do-modulo-multiplication-with-primitive-types)
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;
}
m. Furthermore, it's not possible to use large integer library. – Deccan