モンモール数を包除原理で計算する

　モンモール数を計算するときは漸化式を使えば計算できるが、実は包除原理を使っても計算できるのではないかと思ったので試してみた。

　予想どおり包除原理でも計算できるようだ。
計算量は漸化式を使うとO(N)、包除原理を使うとO(N²)なので漸化式の方が高速。

#include <iostream>
#include <cassert>

using namespace std;

const long long MOD = 1e9 + 7;
const int N = 1000;
long long f[N];
long long comb[N][N];
long long x[N];
long long y[N];

// 漸化式で解く
void solveWithRec() {
  x[0] = 1, x[1] = 0;
  for (int i = 2; i < N; i++)
    x[i] = (i-1) * (x[i-1] + x[i-2]) % MOD;
}

// 包除原理で解く
void solveWithIncExc() {
  f[0] = 1;
  for (int i = 1; i < N; i++)
    f[i] = i * f[i-1] % MOD;

  comb[0][0] = 1;
  for (int i = 1; i < N; i++) {
    comb[i][0] = comb[i][i] = 1;
    for (int j = 1; j < i; j++)
      comb[i][j] = (comb[i-1][j] + comb[i-1][j-1]) % MOD;
  }
  for (int i = 0; i < N; i++) {
    y[i] = f[i];
    for (int j = 1; j <= i; j++) {
      if (j % 2)
        y[i] -= comb[i][j] * f[i-j] % MOD;
      else
        y[i] += comb[i][j] * f[i-j] % MOD;
      y[i] = (y[i] % MOD + MOD) % MOD;
    }
  }
}

int main(int argc, char *argv[])
{
  solveWithRec();
  solveWithIncExc();

  for (int i = 0; i < N; i++)
    assert(x[i] == y[i]);

  cout << "same!" << endl;
  
  return 0;
}