数论板子（待续）
一、质数：
1、试除法判质数：
bool is_prime(int n){ if(n < 2) return false; for(int i = 2; i <= n / i; i ++){ if(n % i == 0) return false; } return true; }
2、试除法分解质因数
void divide(int n){ for(int i = 2; i <= n / i; i ++){ if(n % i == 0){ int s = 0; while(n % i == 0){ n /= i; s ++; } v.push_back({i, s}); } } if(n > 1){ v.push_back({n, 1}); } }
3、朴素筛法
把所有数的倍数全部筛掉
int primes[N], cnt; bool st[N]; void get_primes(int n){ for (int i = 2; i <= n; i++){ if (!st[i]) primes[cnt++] = i; for (int j = i + i; j <= n; j += i){ st[j] = true; } } }
4、埃氏筛法
把质数的倍数全部筛掉
void get_primes(int n) { for (int i = 2; i <= n; i++) { if (!st[i]) { primes[cnt++] = i; for (int j = i + i; j <= n; j += i) { st[j] = true; } } } }
5、线性筛法：
每个合数只会被最小质因子筛掉
void get_primes(int n) { for (int i = 2; i <= n; i++) { if (!st[i]) { primes[cnt++] = i; } for (int j = 0; primes[j] <= n / i; j++) { st[i * primes[j]] = true; if (i % primes[j] == 0) { break; } } } }
二、约数
1、试除法求所有约数
void get_d(int n) { for(int i = 1; i <= n / i; i ++){ if(n % i == 0){ v.push_back(i); if(i != n / i){ v.push_back(n / i); } } } }
2、约数个数
设n = p1k1p2k2 p3k3...... piki
则ans = (k1 + 1) * (k2 + 1) * ... * (ki + 1)
三、欧拉函数
１、求欧拉函数
int phi(int n){ int res = n; for(int i = 2; i <= n / i; i ++){ if(n % i == 0){ res = res / i * (i - 1); while(n % i == 0){ n /= i; } } } if(res > 1){ res = res / n * (n - 1); } return res; }
２、欧拉筛
void get_eulers(int n) { for (int i = 2; i <= n; i++) { if (!st[i]) { primes[cnt++] = i; phi[i] = i - 1; } for (int j = 0; primes[j] <= n / i; j++) { st[i * primes[j]] = true; if (i % primes[j] == 0) { phi[i * primes[j]] = primes[j] * phi[i]; break; } phi[i * primes[j]] = (primes[j] - 1) * phi[i]; } } }
四、快速幂
ll qmi(ll a, ll b, ll p){ ll res = 1; while(b){ if(b & 1){ res = res * a % p; } a = a * a % p; b >>= 1; } return res; }
五、欧几里得
１．欧几里得
ll gcd(ll a, ll b){ return b ? gcd(b, a % b) : a; }
２．扩展欧几里得
ax + by = c
d = gcd(a, b)
x = x0 + k * b / d
y = y0 - k * a / d
ll exgcd(ll a, ll b, ll &x, ll &y){ if(!b){ x = 1, y = 0; return a; } ll d = exgcd(b, a % b, y, x); y -= a / b * x; return d; }
六、中国剩余定理
１．mi 两两互质
ll exgcd(ll a, ll b, ll &x, ll &y){ if(!b){ x = 1, y = 0; return a; } ll d = exgcd(b, a % b, y, x); y -= a / b * x; return d; } ll china(int n, ll M){ ll ans = 0; for(int i = 0; i < n; i ++){ ll Mi = M / m[i]; ll ti, x; exgcd(Mi, m[i], ti, x); ans += a[i] * Mi * ti; } return (ans % M + M) % M; }
2. 不互质
ll exgcd(ll a, ll b, ll &x, ll &y) { if (!b) { x = 1, y = 0; return a; } ll d = exgcd(b, a % b, y, x); y -= a / b * x; return d; } ll n, m[N], a[N], m1, m2, a1, a2, x, y, g, c; ll china() { m1 = m[0], a1 = a[0]; for (int i = 1; i < n; i++) { m2 = m[i], a2 = a[i]; c = a2 - a1; g = exgcd(m1, m2, x, y); if (c % g) { a1 = -1; break; } x *= c / g; ll t = m2 / g; x = (x % t + t) % t; a1 += m1 * x; m1 *= m2 / g; } return a1; }
七、高斯消元
1. 求浮点数解（保证有唯一解）
void gauss() { int c, r; for (c = 1, r = 1; c <= n; c++) { int t = r; for (int i = r; i <= n; i++) { if (fabs(b[i][c]) > fabs(b[t][c])) { t = i; } } for (int i = c; i <= n + 1; i++) { swap(b[t][i], b[r][i]); } for (int i = n + 1; i >= c; i--) { b[r][i] /= b[r][c]; } for (int i = r + 1; i <= n; i++) { if (fabs(b[i][c]) > eps) { for (int j = n + 1; j >= c; j--) { b[i][j] -= b[r][j] * b[i][c]; } } } r++; } for (int i = n; i >= 1; i--) { for (int j = i + 1; j <= n; j++) { b[i][n + 1] -= b[i][j] * b[j][n + 1]; } } }
2.异或（开关问题）
int gauss(){ int r, c; for(r = 1, c = 1; c <= n; c ++){ int t = r; for(int i = r + 1; i <= n; i ++){ if(a[i][c]){ t = i; break; } } if(!a[t][c]) continue; for(int i = c; i <= n + 1; i ++){ swap(a[t][i], a[r][i]); } for(int i = r + 1; i <= n; i ++){ if(a[i][c]){ for(int j = n + 1; j >= c; j --){ a[i][j] ^= a[r][j]; } } } r ++; } int res = 1; if(r < n + 1){ for(int i = r; i <= n; i ++){ if(a[i][n + 1]) return -1; res *= 2; } } return res; }
八、整数分块
for(int l = 1, r; l <= m; l = r + 1){ r = m / (m / l); res = res + (ll)(sum[r] - sum[l - 1]) * (m / l) * (m / l); }
九、快速傅里叶变换 FFT
const int N = 3e6 + 10;
const double PI = acos(-1.0);
int n, m;
ll ans[N];
int rev[N];
struct Complex{
double x, y;
Complex(double _x = 0.0, double _y = 0.0){
x = _x, y = _y;
}
Complex operator + (const Complex & b){
return Complex(x + b.x, y + b.y);
}
Complex operator - (const Complex &b){
return Complex(x - b.x, y - b.y);
}
Complex operator * (const Complex &b){
return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
}
};
void change(Complex y[], int len){
for(int i = 0; i < len; i ++){
rev[i] = rev[i >> 1] >> 1;
if(i & 1) rev[i] |= (len >> 1);
}
for(int i = 0; i < len; i ++){
if(i < rev[i]) swap(y[i], y[rev[i]]);
}
}
void fft(Complex y[], int len, int on){
change(y, len);
for(int h = 2; h <= len; h <<= 1){
Complex wn(cos(2 * PI / h), sin(2 * PI * on / h));
for(int j = 0; j < len; j += h){
Complex w(1, 0);
for(int k = j; k < j + h / 2; k ++){
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w * wn;
}
}
}
if(on == -1){
for(int i = 0; i < len; i ++){
y[i].x /= len;
}
}
}
FFT求快速幂
void qmi(int b){ while(b){ if(b & 1){ fft(x1, len, 1); fft(a, len, 1); for(int i = 0; i < len; i ++) x1[i] = x1[i] * a[i]; fft(x1, len, -1); fft(a, len, -1); for(int i = 0; i < len; i ++) mid[i] = (int)(x1[i].x + 0.5); for(int i = 0; i < len; i ++){ mid[i + 1] += mid[i] / 10; mid[i] %= 10; } for(int i = 0; i < len; i ++){ x1[i] = Complex(mid[i], 0); } for(int i = 0; i < len; i ++) mid[i] = (int)(a[i].x + 0.5); for(int i = 0; i < len; i ++){ mid[i + 1] += mid[i] / 10; mid[i] %= 10; } for(int i = 0; i < len; i ++){ a[i] = Complex(mid[i], 0); } } fft(a, len, 1); for(int i = 0; i < len; i ++){ a[i] = a[i] * a[i]; } fft(a, len, -1); for(int i = 0; i < len; i ++) mid[i] = (int)(a[i].x + 0.5); for(int i = 0; i < len; i ++){ mid[i + 1] += mid[i] / 10; mid[i] %= 10; } for(int i = 0; i < len; i ++){ a[i] = Complex(mid[i], 0); } b >>= 1; } }
十、快速数论变换 NTT
void ntt(ll y[], int len, int on){ change(y, len); for(int h = 2; h <= len; h <<= 1){ ll gn = qmi(3, (mod - 1) / h, mod); //3是原根 if(on == -1) gn = qmi(gn, mod - 2, mod); for(int j = 0; j < len; j += h){ ll g = 1; for(int k = j; k < j + h / 2; k ++){ ll u = y[k]; ll t = g * y[k + h / 2] % mod; y[k] = (u + t) % mod; y[k + h / 2] = (u - t + mod) % mod; g = g * gn % mod; } } } if(on == -1){ ll inv = qmi(len, mod - 2, mod); for(int i = 0; i < len; i ++){ y[i] = y[i] * inv % mod; } } }
NTT求快速幂
void ntt_qmi(ll b, int len){ int mm = m - 1; while(b){ ntt(a, len, 1); if(b & 1){ ntt(res, len, 1); for(int i = 0; i < len; i ++){ res[i] = res[i] * a[i] % mod; } ntt(res, len, -1); for(int i = mm; i < len; i ++){ res[i % mm] = (res[i % mm] + res[i]) % mod; res[i] = 0; } } for(int i = 0; i < len; i ++){ a[i] = a[i] * a[i] % mod; } ntt(a, len, -1); for(int i = mm; i < len; i ++){ a[i % mm] = (a[i % mm] + a[i]) % mod; a[i] = 0; } b >>= 1; } }
三模数NTT+快速幂
#include <map>
#include <set>
#include <array>
#include <queue>
#include <stack>
#include <cmath>
#include <vector>
#include <cstdio>
#include <cstring>
#include <sstream>
#include <iostream>
#include <stdlib.h>
#include <algorithm>
#include <unordered_map>
using namespace std;
typedef long long ll;
typedef pair<int, int> PII;
#define sd(a) scanf("%d", &a)
#define sdd(a, b) scanf("%d%d", &a, &b)
#define slld(a) scanf("%lld", &a)
#define slldd(a, b) scanf("%lld%lld", &a, &b)
#define m1 998244353
#define m2 469762049
#define m3 1004535809
const int N = 3e2 + 10;
const int M = 2e7 + 20;
const int mod = 20170408;
const int INF = 0x3f3f3f3f;
const double PI = acos(-1.0);
const int Mod[] = {998244353, 469762049, 1004535809};
int n, m, p;
int rev[N];
ll vis[N], h[N];
int primes[M], cnt = 0;
bool st[M];
void get(int n)
{
st[1] = true;
for (int i = 2; i <= n; i++)
{
if (!st[i])
primes[cnt++] = i;
for (int j = 0; primes[j] <= n / i; j++)
{
st[i * primes[j]] = true;
if (i % primes[j] == 0)
{
break;
}
}
}
}
ll qmi(ll a, ll b, ll p)
{
ll res = 1;
while (b)
{
if (b & 1)
res = res * a % p;
a = a * a % p;
b >>= 1;
}
return res;
}
void change(ll y[], int len)
{
for (int i = 0; i < len; i++)
{
rev[i] = rev[i >> 1] >> 1;
if (i & 1)
rev[i] |= (len >> 1);
}
for (int i = 0; i < len; i++)
{
if (i < rev[i])
swap(y[i], y[rev[i]]);
}
}
void ntt(ll y[], int len, int on, ll MOD)
{
change(y, len);
for (int h = 2; h <= len; h <<= 1)
{
ll wn = qmi(3, (MOD - 1) / h, MOD);
if (on == -1)
wn = qmi(wn, MOD - 2, MOD);
for (int j = 0; j < len; j += h)
{
ll w = 1;
for (int k = j; k < j + h / 2; k++)
{
ll u = y[k];
ll t = w * y[k + h / 2] % MOD;
y[k] = (u + t) % MOD;
y[k + h / 2] = (u - t + MOD) % MOD;
w = w * wn % MOD;
}
}
}
if (on == -1)
{
ll inv = qmi(len, MOD - 2, MOD);
for (int i = 0; i < len; i++)
{
y[i] = y[i] * inv % MOD;
}
}
}
ll mult(ll a, ll b, ll p)
{
ll res = 0;
while (b)
{
if (b & 1)
res = (res + a) % p;
a = (a + a) % p;
b >>= 1;
}
return res;
}
ll A[N], B[N], C[N], D[N];
void mul(ll a[], ll b[], ll res[], ll len)
{
memcpy(A, a, sizeof(A));
memcpy(B, a, sizeof(B));
memcpy(C, a, sizeof(C));
memcpy(D, b, sizeof(D));
ntt(A, len, 1, Mod[0]);
ntt(D, len, 1, Mod[0]);
for (int i = 0; i < len; i++)
{
A[i] = A[i] * D[i] % Mod[0];
}
ntt(A, len, -1, Mod[0]);
memcpy(D, b, sizeof(D));
ntt(B, len, 1, Mod[1]);
ntt(D, len, 1, Mod[1]);
for (int i = 0; i < len; i++)
{
B[i] = B[i] * D[i] % Mod[1];
}
ntt(B, len, -1, Mod[1]);
memcpy(D, b, sizeof(D));
ntt(C, len, 1, Mod[2]);
ntt(D, len, 1, Mod[2]);
for (int i = 0; i < len; i++)
{
C[i] = C[i] * D[i] % Mod[2];
}
ntt(C, len, -1, Mod[2]);
ll M12 = 1ll * m1 * m2;
ll inv2 = qmi(m2, m1 - 2, m1);
ll inv1 = qmi(m1, m2 - 2, m2);
ll mul2 = 1ll * m2 * inv2 % M12;
ll mul1 = 1ll * m1 * inv1 % M12;
ll inv = qmi(M12 % m3, m3 - 2, m3);
ll m12 = M12 % mod;
ll c1, c2, c3, c4, q;
for (int i = 0; i <= (p << 1); i++)
{
c1 = A[i], c2 = B[i], c3 = C[i];
c4 = (mult(c1, mul2, M12) + mult(c2, mul1, M12)) % M12;
q = ((c3 - c4) % m3 + m3) % m3 * inv % m3;
res[i] = (q * m12 % mod + c4) % mod;
}
for (int i = p; i < len; i++)
{
res[i % p] = (res[i % p] + res[i]) % mod;
res[i] = 0;
}
}
ll res[N];
void qmi_ntt(ll y[], int len, int n)
{
memset(res, 0, sizeof(res));
res[0] = 1;
while (n)
{
if (n & 1)
{
mul(res, y, res, len);
}
mul(y, y, y, len);
n >>= 1;
}
}
ll mid[N], ans[3], ans1, ans2;
void solve()
{
cin >> n >> m >> p;
get(m);
for (int i = 1; i <= m; i++)
{
vis[i % p]++;
if (st[i])
h[i % p]++;
}
int len = 1;
while (len <= p + p - 1)
len <<= 1;
qmi_ntt(vis, len, n);
ans1 = res[0];
qmi_ntt(h, len, n);
ans1 = (ans1 - res[0] + mod) % mod;
cout << ans1 << "\n";
}
int main()
{
#ifdef ONLINE_JUDGE
#else
freopen("/home/jungu/code/in.txt", "r", stdin);
// freopen("/home/jungu/桌面/11.21/2/in9.txt", "r", stdin);
#endif
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
int T = 1;
// sd(T);
// cin >> T;
while (T--)
{
solve();
}
return 0;
}
