[快速数论变换 NTT]

时间:2023-03-09 02:46:04
[快速数论变换 NTT]

先粘一个模板。这是求高精度乘法的

#include <bits/stdc++.h>
#define maxn 1010
using namespace std; char s[maxn]; typedef long long ll;
ll A[maxn], B[maxn]; const int md = 998244353, G = 3; int n; ll power_mod(ll a, ll b){
ll ret = 1;
while(b > 0){
if(b & 1)ret = ret * a % md;
b >>= 1;
a = a * a % md;
}return ret;
} void NTT(ll A[], int n, int type){
for(int i = 0, j = 0; i < n; i ++){
if(i < j)swap(A[i], A[j]);
for(int t = n >> 1; (j ^= t) < t; t >>= 1);
} for(int k = 2; k <= n; k <<= 1){
ll wn = type > 0 ? power_mod(G, (md - 1) / k) : power_mod(G, (md - 1) - (md - 1) / k);
for(int i = 0; i < n; i += k){
ll w = 1;
for(int j = 0; j < k >> 1; j ++){
ll T = w * A[i+j+(k>>1)] % md;
A[i+j+(k>>1)] = (A[i+j] - T) % md;
A[i+j] = (A[i+j] + T) % md;
w = w * wn % md;
}
}
} if(type < 0){
ll inv = power_mod(n, md - 2);
for(int i = 0; i < n; i ++)
A[i] = A[i] * inv % md;
}
} int main(){
freopen("mul.in", "r", stdin);
freopen("mul.out", "w", stdout);
scanf("%s", s);
int n1 = strlen(s);
for(int i = n1 - 1; ~i; i --)
A[n1-i-1] = s[i] ^ 48;
scanf("%s", s);
int n2 = strlen(s);
for(int i = n2 - 1; ~i; i --)
B[n2-i-1] = s[i] ^ 48;
for(n = 1; n <= n1 + n2; n <<= 1);
NTT(A, n, 1), NTT(B, n, 1);
for(int i = 0; i < n; i ++)
A[i] = A[i] * B[i] % md;
NTT(A, n, -1);
for(int i = 0; i < n; i ++){
A[i] = (A[i] + md) % md;
A[i+1] += A[i] / 10, A[i] %= 10;
}
while(n && A[n] == 0) n --;
for(int i = n; i >= 0; i --)
printf("%lld", A[i]);
return 0;
}

[SDOI 2015]序列统计

取指标变成加法问题,构造生成函数F(x),F(x)的n次方就是答案。

注意指标要mod M-1

还有注意0的问题,如果x!=0就不用管数据中的0,否则需要两种情况求和(然而数据中x并不等于0我就没写)

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define maxn 50010
using namespace std;
typedef long long ll;
const int md = 1004535809, G = 3; int N, M, X, S, a[maxn], Log[maxn], n;
ll F[maxn], H[maxn], ret[maxn]; int p[maxn], primes, vis[maxn]; ll power_mod(ll a, ll b, ll md){
ll ret = 1;
while(b > 0){
if(b & 1) ret = ret * a % md;
b >>= 1;
a = a * a % md;
}return ret;
} void pre_prime(){
int n = 20000;
for(int i = 2; i <= n; i ++){
if(!vis[i])p[++ primes] = i;
for(int j = 1; j <= primes; j ++){
if(i * p[j] > n)break;
vis[i * p[j]] = true;
if(i % p[j] == 0)break;
}
}
} int getG(int M){
pre_prime();
int ret = 2;
for(; ; ret ++){
bool flag = true;
for(int i = 1; i <= primes; i ++){
if(M-1 < p[i])break;
if((M-1) % p[i] == 0){
if(power_mod(ret, (M - 1) / p[i], M) == 1){
flag = false;
break;
}
}
}
if(flag)return ret;
}
} void NTT(ll A[], int n, int type){
for(int i = 0, j = 0; i < n; i ++){
if(i > j)swap(A[i], A[j]);
for(int t = n >> 1; (j ^= t) < t; t >>= 1);
} for(int k = 2; k <= n; k <<= 1){
ll wn = type > 0 ? power_mod(G, (md-1)/k, md) : power_mod(G, md-1-(md-1)/k, md);
for(int i = 0; i < n; i += k){
ll w = 1;
for(int j = 0; j < k >> 1; j ++){
ll T = w * A[i+j+(k>>1)];
A[i+j+(k>>1)] = (A[i+j] - T + md) % md;
A[i+j] = (A[i+j] + T) % md;
w = w * wn % md;
}
}
} if(type < 0){
ll inv = power_mod(n, md - 2, md);
for(int i = 0; i < n; i ++)
A[i] = A[i] * inv % md;
}
} void CF(ll A[]){
for(int i = 0; i < n; i ++)H[i] = F[i];
NTT(A, n, 1), NTT(H, n, 1);
for(int i = 0; i < n; i ++)
A[i] = A[i] * H[i] % md;
NTT(A, n, -1);
for(int i = (M - 1); i < n; i ++)
(A[i % (M - 1)] += A[i]) %= md, A[i] = 0;
} void power(int b){
ret[0] = 1;
while(b){
if(b & 1)CF(ret);
b >>= 1;
CF(F);
}
} int main(){
scanf("%d%d%d%d", &N, &M, &X, &S);
for(int i = 1; i <= S; i ++)
scanf("%d", &a[i]); for(n = 1; n <= M; n <<= 1); n <<= 1; ll nw = 1, Gn = getG(M);
for(int i = 0; i < M-1; i ++)
Log[nw] = i, nw = nw * Gn % M; for(int i = 1; i <= S; i ++)
if(a[i])F[Log[a[i]]] = 1; power(N);
printf("%lld\n", (ret[Log[X]] + md) % md);
return 0;
}

[COGS 2287]疯狂的机器人

这里有Catalan数的链接

#include <bits/stdc++.h>
#define maxn 300010
using namespace std;
typedef long long ll;
int n;
const int md = 998244353, G = 3; ll power_mod(ll a, ll b = md - 2){
ll ret = 1;
while(b > 0){
if(b & 1)ret = ret * a % md;
b >>= 1;
a = a * a % md;
}return ret;
} ll fac[maxn], inv[maxn], f[maxn]; void NTT(ll* A, int n, int type){
for(int i = 0, j = 0; i < n; i ++){
if(i > j)swap(A[i], A[j]);
for(int t = n >> 1; (j ^= t) < t; t >>= 1);
} for(int k = 2; k <= n; k <<= 1){
ll wn = power_mod(G, type > 0 ? (md-1) / k : md-1 - (md-1) / k);
for(int i = 0; i < n; i += k){
ll w = 1;
for(int j = 0; j < k >> 1; j ++){
ll T = w * A[i+j+(k>>1)] % md;
A[i+j+(k>>1)] = (A[i+j] - T) % md;
A[i+j] = (A[i+j] + T) % md;
w = w * wn % md;
}
}
} if(type == -1){
ll inv = power_mod(n);
for(int i = 0; i < n; i ++)
A[i] = A[i] * inv % md;
}
} int main(){
freopen("crazy_robot.in", "r", stdin);
freopen("crazy_robot.out", "w", stdout);
scanf("%d", &n);
fac[0] = inv[0] = 1;
for(int i = 1; i <= 2 * n; i ++)
fac[i] = fac[i-1] * i % md;
inv[2 * n] = power_mod(fac[2 * n]);
for(int i = 2 * n - 1; i >= 1; i --)
inv[i] = inv[i+1] * (i+1) % md;
for(int i = 0; i <= n; i ++){
if(i & 1){f[i] = 0;continue;}
f[i] = fac[i] * inv[i / 2] % md * inv[i / 2 + 1] % md * inv[i] % md;
} int _N = 1;
for(; _N <= n + n; _N <<= 1);
NTT(f, _N, 1);
for(int i = 0; i < _N; i ++)
f[i] = f[i] * f[i] % md;
NTT(f, _N, -1);
ll ans = 0;
for(int i = 0; i <= n; i ++)
(ans += inv[n - i] * f[i]) %= md;
cout << (ans * fac[n] % md + md) % md << endl;
return 0;
}

[NOI十连测3007]卡常大法好。

#include <cstring>
#include <cstdio>
#define maxn 150010 const int md = 998244353, G = 3; int R;
#define fastcall __attribute__((optimize("-Os")))
#define IL __inline__ __attribute__((always_inline))
fastcall IL int mul_mod(int a, int b){
__asm__ __volatile__ ("\tmull %%ebx\n\tdivl %%ecx\n" :"=d"(R):"a"(a),"b"(b),"c"(md));
return R;
} int n, N; int f[maxn], g[maxn], h[maxn], fac[maxn], inv[maxn], pw[maxn], tmp[maxn]; fastcall IL int power_mod(int a, long long b = md - 2){
int ret = 1;
while(b > 0){
if(b & 1)ret = mul_mod(ret, a);
b >>= 1;
a = mul_mod(a, a);
}return ret;
} fastcall IL void swap(int& a, int& b){
a ^= b ^= a ^= b;
} fastcall IL void NTT(int* A, int n, int type){
for(int i = 0, j = 0; i < n; i ++){
if(i > j)swap(A[i], A[j]);
for(int t = n >> 1; (j ^= t) < t; t >>= 1);
} for(int k = 2; k <= n; k <<= 1){
int wn = power_mod(G, type > 0 ? (md-1)/k : (md-1)-(md-1)/k);
for(int i = 0, m = k >> 1; i < n; i += k){
int w = 1;
for(int j = 0; j < k >> 1; j ++){
int T = mul_mod(w, A[i+j+m]);
A[i+j+m] = (A[i+j] - T + md) % md;
A[i+j] = (A[i+j] + T) % md;
w = mul_mod(w, wn);
}
}
} if(type < 0){
int inv = power_mod(n);
for(int i = 0; i < n; i ++)
A[i] = mul_mod(A[i], inv);
}
} fastcall void Getinv(int* a, int* b, int n){
if(n == 1){b[0] = power_mod(a[0]); return;}
Getinv(a, b, n >> 1);
memcpy(tmp, a, sizeof(a[0]) * n);
NTT(tmp, n << 1, 1);
NTT(b, n << 1, 1);
for(int i = 0; i < n << 1; i ++)
b[i] = mul_mod(b[i], (2ll - mul_mod(b[i], tmp[i]) + md) % md);
NTT(b, n << 1, -1);
memset(b + n, 0, sizeof(a[0]) * n);
} int F[30010][16];
int K; fastcall void cdq(int l, int r){
if(l == r)return;
int mid = l + r >> 1, len = r - l + 1;
cdq(l, mid);
for(n = 1; n <= len; n <<= 1);
for(int i = 0; i < n; i ++)f[i] = h[i] = 0;
for(int i = l; i <= mid; i ++)
f[i-l] = mul_mod(F[i][K-1] + F[i][K], inv[i]);
h[0] = 0;
for(int i = 1; i < n; i ++)
h[i] = mul_mod(g[i], inv[i-1]);
if(n<=32) {
for(int i = 0; i < n; i ++)tmp[i] = 0;
for(int i = 0 ; i < n ; ++ i) {
for(int j = 0 ; j <= i ; ++ j) {
tmp[i] += mul_mod(h[j], f[i-j]);
if(tmp[i]>=md) tmp[i] -= md;
}
}
for(int i = 0 ; i < n ; ++ i) h[i] = tmp[i];
} else {
NTT(f, n, 1), NTT(h, n, 1);
for(int i = 0; i < n; i ++)
h[i] = mul_mod(h[i], f[i]);
NTT(h, n, -1);
}
for(int i = mid + 1; i <= r; i ++){
F[i][K] += mul_mod(h[i-l], fac[i-1]);
if(F[i][K] >= md)F[i][K] -= md;
}
cdq(mid + 1, r);
} int main(){
N = 30000;
for(n = 1; n <= N; n <<= 1);
fac[0] = inv[0] = 1;
for(int i = 1; i <= n; i ++)fac[i] = mul_mod(fac[i - 1], i);
inv[n] = power_mod(fac[n]);
for(int i = n - 1; i >= 1; i --)inv[i] = mul_mod(inv[i + 1], i + 1);
for(long long i = 0; i <= n; i ++)pw[i] = power_mod(2, i * (i - 1) / 2);
for(int i = 1; i < n; i ++)h[i] = mul_mod(pw[i], inv[i - 1]);
for(int i = 0; i < n; i ++)f[i] = mul_mod(pw[i], inv[i]);
Getinv(f, g, n);
n <<= 1;
NTT(g, n, 1);
NTT(h, n, 1);
for(int i = 0; i < n; i ++)
g[i] = mul_mod(g[i], h[i]);
NTT(g, n, -1);
for(int i = 1; i <= N; i ++)
g[i] = mul_mod(g[i], fac[i - 1]); for(int i = 0; i <= N; i ++)F[i][0] = pw[i];
for(int i = 1; i <= 15; i ++)K = i, cdq(0, N); static int S[20][20];
S[0][0] = 1;
for(int i = 1; i <= 15; i ++)
for(int j = 1; j <= i; j ++){
S[i][j] = mul_mod(S[i-1][j], j) + S[i-1][j-1];
if(S[i][j] >= md)S[i][j] -= md;
} int test, n, m;
scanf("%d", &test);
while(test --){
scanf("%d%d", &n, &m);
int ans = 0;
for(int i = 0; i <= m; i ++){
ans += mul_mod(mul_mod(S[m][i], F[n][i]), fac[i]);
if(ans >= md)ans -= md;
}
printf("%d\n", ans);
}
return 0;
}