而且, 这题用NTT + FFT貌似不可行.
NTT模数不对, FFT爆精度.
首先想一个暴力算法:
设f[i][j] 表示选i个, 在模p 意义下是j的方案数. f[i][j]=f[i−1][k]×l(k+l)mod 复杂度O(nm)
然后发现一个优化, 可以只记录一下每个数模p的结果, 记录一下cnt
然后可以优化到差不多O(np)
然后又发现一个性质 f[k] * f[j] = f[j + k] 然后类比 x^a * x^b = x^{a + b} 如何快速计算x^{a + b}呢
快速幂就好了.
同理, 计算f 也可以快速幂.
然后数组可以用滚动数组优化.
最终复杂度是O(\log n * p) 再算上一些奇奇怪怪的预处理.
#include <bits/stdc++.h>
using namespace std;
// #define int long long
#define ll long long
const int N = 2e2 + 10;
const int L = 2e7 + 10;
const int mod = 20170408;
ll F[N], G[N], f[N], g[N], n, m, p;
int v[L], q[L / 10], cnt[N];
inline void multiply(ll *a, ll *b)
{
static ll res[N];
for (int i = 0; i < p; i++)
for (int j = 0; j < p; j++)
{
res[i + j] = (a[i] * b[j] % mod + res[i + j]) % mod;
}
for (int i = 0; i < p; i++)
{
a[i] = (res[i] + res[i + p]) % mod;
res[i] = res[i + p] = 0;
}
}
signed main()
{
cin >> n >> m >> p;
F[1] = G[1] = f[1] = g[1] = 1;
for (int i = 2; i <= m; i++)
{
f[i % p]++;
g[i % p]++;
F[i % p]++;
G[i % p]++;
if (v[i] == 0)
{
v[i] = i;
// cout << i << ' ';
G[i % p]--;
g[i % p]--;
q[++q[0]] = i;
}
for (int j = 1; j <= q[0]; j++)
{
if (q[j] > v[i] || q[j] * i > m) break;
v[i * q[j]] = q[j];
}
}
n--;
while (n)
{
if (n & 1)
{
multiply(F, f);
multiply(G, g);
}
n >>= 1;
multiply(f, f);
multiply(g, g);
}
cout << (F[0] - G[0] + mod) % mod << endl;
return 0;
}