数学板子(后面会陆续更新的)
线性筛逆元
void init()
{
inv[0] = inv[1] = fac[0] = fac[1] = INV[0] = INV[1] = 1;
for (int i = 2; i < N; i++)
{
INV[i] = 1ll * (mod - mod / i) * INV[mod % i] % mod;
inv[i] = 1ll * inv[i - 1] * INV[i] % mod;
fac[i] = 1ll * fac[i - 1] * i % mod;
}
}
分解质因数
#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 10;
int primes[N],cnt;
bool st[N];
void init() {
st[1] = true;
st[0] = true;
for(int i=2;i<N;i++)
{
if(!st[i]) primes[cnt ++ ] = i;
for(int j=0;primes[j] * i < N ; j ++ )
{
st[primes[j] * i] = true;
if(i % primes[j] == 0)
break;
}
}
}
int main(){
init();
int T;
cin >> T;
while(T -- )
{
int n;
cin >> n;
for(int i=0;primes[i]*primes[i]<=n;i++)
{
int x = primes[i];
if(n % x == 0)
{
int c = 0;
while(n % x == 0) c ++ ,n /= x;
cout << x << " " << c << endl;
}
}
if(n > 1) cout << n << " " << 1 << endl;
cout << endl;
}
return 0;
}
NTT
- 多项式乘法
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int g = 3, gi = 332748118, mod = 998244353; // g为原根,gi为g的逆元
const int M = 1 << 22;
int tot,bit;
int a[M], b[M],rev[M];
int n,m;
int qmi(int a, int b) //快速幂,a为底数,k为指数
{
int res = 1;
while (b)
{
if (b & 1) res = 1ll * res * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return res;
}
void NTT(int *a, int n, int type) // NTT,type=1时系数表示法转点值表示法,否则点值表示法转系数表示法
{
for (int i = 0; i < n; ++i) //先将a变成最后的样子
if (i < rev[i])
swap(a[i], a[rev[i]]);
for (int i = 1; i < n; i <<= 1)
{ //在这之前NTT与FFT无异
int gn = qmi(type ? g : gi, (mod - 1) / (i << 1)); //单位原根g_n
for (int j = 0; j < n; j += (i << 1)) //枚举具体区间,j也就是区间右端点
{
int g0 = 1;
for (int k = 0; k < i; ++k, g0 = 1ll * g0 * gn % mod) //合并,记得要取模
{
int x = a[j + k], y = 1ll * g0 * a[i + j + k] % mod;
a[j + k] = (x + y) % mod;
a[i + j + k] = (x - y + mod) % mod;
}
}
}
}
int main()
{
cin >> n >> m;
while((1 << bit) < n + m + 1) bit ++ ;
tot = 1 << bit;
for(int i=0;i<tot;i++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
memset(a,0,sizeof a);
memset(b,0,sizeof b);
for(int i=0;i<=n;i++) cin >> a[i];
for(int i=0;i<=m;i++) cin >> b[i];
NTT(a, tot, 1);
NTT(b, tot, 1);
for (int i = 0; i < tot; ++i)
a[i] = 1ll * a[i] * b[i] % mod; // O(n)乘法
NTT(a, tot, 0); //点值表示法转系数表示法
int inv = qmi(tot, mod - 2); // inv为len的逆元(费马小定理求逆元)
for (int i = 0; i <= n + m; ++i) //输出
printf("%d ",1ll * a[i] * inv % mod);
puts("");
return 0;
}
- 多项式求逆和多项式快速幂
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 400003, mod = 998244353, G = 3, Gi = 332748118;
int n, k, A[N]; // n表示有多少项 0 ~ n-1,k表示快速幂次数
int ans[N];
int rev[N];
int Ln[N];
int Exp[N];
inline int read(){
int ch = getchar(), res = 0;
while(ch < '0' || ch > '9') ch = getchar();
while(ch >= '0' && ch <= '9'){
res = (res * 10ll + ch - '0') % mod;
ch = getchar();
}
return res;
}
int qmi(int a, int b){
int res = 1;
while(b){
if(b & 1) res = (LL) res * a % mod;
a = (LL) a * a % mod;
b >>= 1;
}
return res;
}
void NTT(int *A, int limit, int type){
for(int i = 0;i < limit;i ++)
if(i < rev[i]) swap(A[i], A[rev[i]]);
for(int mid = 1;mid < limit;mid <<= 1){
int Wn = qmi(type == 1 ? G : Gi, (mod - 1) / (mid << 1));
for(int j = 0;j < limit;j += mid << 1){
int w = 1;
for(int k = 0;k < mid;k ++, w = (LL) w * Wn % mod){
int x = A[j + k], y = (LL) w * A[j + k + mid] % mod;
A[j + k] = (x + y) % mod;
A[j + k + mid] = (x - y + mod) % mod;
}
}
}
if(type == -1){
int inv = qmi(limit, mod - 2);
for(int i = 0;i < limit;i ++) A[i] = (LL) A[i] * inv % mod;
}
}
void poly_inv(int *A, int deg){
static int tmp[N];
if(deg == 1){
ans[0] = qmi(A[0], mod - 2);
return;
}
poly_inv(A, deg + 1 >> 1);
int limit = 1, L = -1;
while(limit <= (deg << 1)){limit <<= 1; L ++;}
for(int i = 0;i < limit;i ++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
for(int i = 0;i < deg;i ++) tmp[i] = A[i];
for(int i = deg;i < limit;i ++) tmp[i] = 0;
NTT(tmp, limit, 1); NTT(ans, limit, 1);
for(int i = 0;i < limit;i ++) ans[i] = (2 - (LL) ans[i] * tmp[i] % mod + mod) % mod * ans[i] % mod;
NTT(ans, limit, -1);
for(int i = deg;i < limit;i ++) ans[i] = 0;
}
void poly_Ln(int *A, int deg){
static int tmp[N];
poly_inv(A, deg);
for(int i = 1;i < deg;i ++) tmp[i - 1] = (LL) i * A[i] % mod;
tmp[deg - 1] = 0;
int limit = 1, L = -1;
while(limit <= (deg << 1)){limit <<= 1; L ++;}
for(int i = 0;i < limit;i ++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
NTT(ans, limit, 1); NTT(tmp, limit, 1);
for(int i = 0;i < limit;i ++) Ln[i] = (LL) ans[i] * tmp[i] % mod;
NTT(Ln, limit, -1);
for(int i = deg + 1;i < limit;i ++) Ln[i] = 0;
for(int i = deg;i;i --) Ln[i] = (LL) Ln[i - 1] * qmi(i, mod - 2) % mod;
Ln[0] = 0;
for(int i = 0;i < limit;i ++) ans[i] = tmp[i] = 0;
}
void poly_Exp(int *A, int deg){
if(deg == 1){
Exp[0] = 1;
return;
}
poly_Exp(A, deg + 1 >> 1);
poly_Ln(Exp, deg);
for(int i = 0;i < deg;i ++) Ln[i] = (A[i] + (i == 0) - Ln[i] + mod) % mod;
int limit = 1, L = -1;
while(limit <= (deg << 1)){limit <<= 1; L ++;}
for(int i = 0;i < limit;i ++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
NTT(Exp, limit, 1); NTT(Ln, limit, 1);
for(int i = 0;i < limit;i ++) Exp[i] = (LL) Exp[i] * Ln[i] % mod;
NTT(Exp, limit, -1);
for(int i = deg;i < limit;i ++) Exp[i] = 0;
for(int i = 0;i < limit;i ++) Ln[i] = ans[i] = 0;
}
void poly_ni(int n){ // 多项式求逆
memset(A,0,sizeof A);
memset(ans,0,sizeof ans);
memset(rev,0,sizeof rev);
memset(Exp,0,sizeof Exp);
memset(Ln,0,sizeof Ln);
for(int i = 1;i < n;i ++) {
int x = read();
A[i] = -x + mod;
}
A[0] = 1 % mod;
poly_inv(A,n);
for(int i=0;i<n;i++) printf("%d ",ans[i]);
}
void poly_mi(int n,int k) { // 多项式快速幂
memset(A,0,sizeof A);
memset(ans,0,sizeof ans);
memset(rev,0,sizeof rev);
memset(Exp,0,sizeof Exp);
memset(Ln,0,sizeof Ln);
for(int i = 0;i < n;i ++) A[i] = read();
poly_Ln(A, n);
for(int i = 0;i < n;i ++) A[i] = (LL) Ln[i] * k % mod, Ln[i] = 0;
poly_Exp(A, n);
for(int i = 0;i < n;i ++) printf("%d ", Exp[i]);
puts("");
}
int main(){
n = read();
poly_ni(n);
}
- 分治NTT
分治NTT处理这个a,b,a_0 = 1,b_0 = 0
// 分治ntt处理ai和bi
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N = 262144, K = 17, mod = 998244353, G = 3; // K设为log(N)
int n, i, x;
int pos[N + 5], A[N + 5], B[N + 5], W[N + 5], g[K + 5], ng[K + 5], inv[N + 5], inv2;
int p[N + 5], invp[N + 5], cost[N + 5], w[N + 5], pre[N + 5];
int ea[N + 5], eb[N + 5], fa[N + 5], fb[N + 5];
int qmi(int a, int b)
{
int res = 1;
while (b)
{
if (b & 1)
res = 1ll * res * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return res;
}
inline void NTT(int *a, int n, int t)
{
for (int i = 1; i < n; i++)
if (i < pos[i])
swap(a[i], a[pos[i]]);
for (int d = 0; (1 << d) < n; d++)
{
int m = 1 << d, m2 = m << 1, _w = t == 1 ? g[d] : ng[d];
for (int i = 0; i < n; i += m2)
for (int w = 1, j = 0; j < m; j++)
{
int &A = a[i + j + m], &B = a[i + j], t = 1LL * w * A % mod;
A = B - t;
if (A < 0)
A += mod;
B = B + t;
if (B >= mod)
B -= mod;
w = 1LL * w * _w % mod;
}
}
if (t == -1)
for (int i = 0, j = inv[n]; i < n; i++)
a[i] = 1LL * a[i] * j % mod;
}
void solve(int l, int r)
{
if (l == r)
{
if (l)
{
ea[l+1]=((ea[l]-(1LL-p[l])*fa[l]%mod*pre[l])%mod+mod)*invp[l]%mod; // fa[l]就是ntt处理出来的那个求和
eb[l+1]=((eb[l]-cost[l]-(1LL-p[l])*fb[l]%mod*pre[l])%mod+mod)*invp[l]%mod;
}
return;
}
int mid = (l + r) >> 1;
solve(l, mid);
int k = 1;
while (k < r - l + 1) k <<= 1;
for (int i = 0; i < k; i++)
A[i] = B[i] = W[i] = 0;
for (int i = l; i <= mid; i++)
{
A[i - l] = ea[i];
B[i - l] = eb[i];
}
for (int i = 1; i <= r - l; i++)
W[i] = w[i];
int j = __builtin_ctz(k) - 1;
for (int i = 0; i < k; i++)
pos[i] = pos[i >> 1] >> 1 | ((i & 1) << j);
NTT(A,k,1);
NTT(B,k,1);
NTT(W,k,1);
for (int i = 0; i < k; i++)
{
A[i] = 1ll * A[i] * W[i] % mod;
B[i] = 1ll * B[i] * W[i] % mod;
}
NTT(A,k,-1);
NTT(B,k,-1);
for (int i = mid + 1; i <= r; i++)
{
fa[i] = (fa[i] + A[i - l]) % mod;
fb[i] = (fb[i] + B[i - l]) % mod;
}
solve(mid + 1, r);
}
int main()
{
// freopen("data.in","r",stdin);
// freopen("data.out","w",stdout);
for (g[K] = qmi(G, (mod - 1) / N), ng[K] = qmi(g[K], mod - 2), i = K - 1; ~i; i--)
g[i] = 1LL * g[i + 1] * g[i + 1] % mod, ng[i] = 1LL * ng[i + 1] * ng[i + 1] % mod;
for (int i = 1; i <= N; i++)
inv[i] = qmi(i, mod - 2);
inv2 = inv[2];
int T;
scanf("%d", &T);
while (T--)
{
scanf("%d", &n);
for (int i = 0; i < n; i++)
{
int x;
scanf("%d%d", &x, &cost[i]);
p[i] = 1ll * x * inv[100] % mod;
invp[i] = 100ll * inv[x] % mod;
}
for(int i=1;i<n;i++) {
scanf("%d",&w[i]);
pre[i] = (pre[i - 1] + w[i]) % mod;
}
for(i=1;i<n;i++)pre[i]=qmi(pre[i],mod-2);
for (i = 0; i <= n; i++) fa[i] = fb[i] = 0;
ea[0] = 1, eb[0] = 0;
ea[1] = 1, eb[1] = (mod-cost[0])%mod;
solve(0,n-1);
int res = 1ll * (mod - eb[n]) * qmi(ea[n],mod - 2) % mod;
printf("%d\n",res);
}
return 0;
}
- 分治ntt小板子
#include <bits/stdc++.h>
#define ll long long
#define Poly vector<int>
using namespace std;
const int g = 3, gi = 332748118, mod = 998244353; // g为原根,gi为g的逆元
const int M = 1 << 22;
int tot,bit;
int rev[M];
int qmi(int a, int b) //快速幂,a为底数,k为指数
{
int res = 1;
while (b)
{
if (b & 1) res = 1ll * res * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return res;
}
void NTT(Poly &a, int type) // NTT,type=1时系数表示法转点值表示法,否则点值表示法转系数表示法
{
int n = a.size();
for (int i = 0; i < n; ++i) //先将a变成最后的样子
if (i < rev[i])
swap(a[i], a[rev[i]]);
for (int i = 1; i < n; i <<= 1)
{ //在这之前NTT与FFT无异
int gn = qmi(type ? g : gi, (mod - 1) / (i << 1)); //单位原根g_n
for (int j = 0; j < n; j += (i << 1)) //枚举具体区间,j也就是区间右端点
{
int g0 = 1;
for (int k = 0; k < i; ++k, g0 = 1ll * g0 * gn % mod) //合并,记得要取模
{
int x = a[j + k], y = 1ll * g0 * a[i + j + k] % mod;
a[j + k] = (x + y) % mod;
a[i + j + k] = (x - y + mod) % mod;
}
}
}
}
Poly operator*(Poly a,Poly b) {
int n = a.size();
int m = b.size();
n --;
m --;
bit = 0;
while((1 << bit) < n + m + 1) bit ++ ;
tot = 1 << bit;
for(int i=0;i<tot;i++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
a.resize(tot);
b.resize(tot);
NTT(a,1);
NTT(b,1);
for (int i = 0; i < tot; ++i)
a[i] = 1ll * a[i] * b[i] % mod; // O(n)乘法
NTT(a, 0); //点值表示法转系数表示法
int inv = qmi(tot, mod - 2); // inv为len的逆元(费马小定理求逆元)
a.resize(n + m + 1);
for (int i = 0; i <= n + m; ++i) //输出
a[i] = 1ll * a[i] * inv % mod;
return a;
}
Poly operator+(const Poly &a,const Poly &b) {
int tot = max(a.size(),b.size());
Poly res(tot,0);
for(int i=0;i<tot;i++) {
if(i < a.size()) res[i] = (res[i] + a[i]) % mod;
if(i < b.size()) res[i] = (res[i] + b[i]) % mod;
}
return res;
}
int a[M],n,m;
Poly odd;
Poly even;
pair<Poly,Poly> merge(int l,int r) {
if(l == r) {
odd.resize(a[l] + 1);
even.resize(1);
even[0] = 1;
for(int i=0;i<a[l];i++) {
odd[i] = 0;
}
odd[a[l]] = 1;
return {Poly(odd),Poly(even)};
}
int mid = l + r >> 1;
auto L = merge(l,mid);
auto R = merge(mid + 1,r);
return {L.first * R.second + L.second * R.first,L.first * R.first + L.second * R.second};
}
int main() {
scanf("%d%d",&n,&m);
for(int i=1;i<=n;i++) scanf("%d",&a[i]);
auto res = merge(1,n).first;
res.resize(m + 1);
printf("%d\n",res[m]);
return 0;
}
- 分治ntt大板子(天使姐姐)
#include <bits/stdc++.h>
#define fp(i, a, b) for (int i = (a), i##_ = (b) + 1; i < i##_; ++i)
#define fd(i, a, b) for (int i = (a), i##_ = (b) - 1; i > i##_; --i)
#define MUL(a, b) ((ll)(a) * (b) % P) // 大模数时这里ll要改成__int128
#define ADD(a, b) (((a) += (b)) >= P ? (a) -= P : 0) // (a += b) %= P
#define SUB(a, b) (((a) -= (b)) < 0 ? (a) += P: 0) // ((a -= b) += P) %= P
#define int long long
using namespace std;
const int N = 3e5 + 5, P = 998244353; //模数
using ll = int64_t;
using Poly = vector<int>;
/*---------------------------------------------------------------------------*/
class Cipolla {
int P, I2 {};
using pll = pair<ll, ll>;
#define X first
#define Y second
ll mul(ll a, ll b) const {
return a * b % P;
}
pll mul(pll a, pll b) const {
return {(a.X * b.X + I2 * a.Y % P * b.Y) % P,
(a.X * b.Y + a.Y * b.X) % P
};
}
template<class T> T POW(T a, int b, T x) {
for (; b; b >>= 1, a = mul(a, a))
if (b & 1) x = mul(x, a);
return x;
}
public:
Cipolla(int p = 0) : P(p) {}
pair<int, int> sqrt(int n) {
int a = rand(), x;
if (!(n %= P)) return {0, 0};
if (POW(n, (P - 1) >> 1, 1ll) == P - 1) return {-1, -1};
while (POW(I2 = ((ll) a * a - n + P) % P, (P - 1) >> 1, 1ll) == 1) a =
rand();
x = (int) POW(pll {a, 1}, (P + 1) >> 1, {1, 0}).X;
if (2 * x > P) x = P - x;
return {x, P - x};
}
#undef X
#undef Y
};
/*---------------------------------------------------------------------------*/
Poly getInv(int L) {
Poly inv(L);
inv[1] = 1;
fp(i, 2, L - 1)
inv[i] = MUL((P - P / i), inv[P % i]);
return inv;
}
ll POW(ll a, ll b = P - 2, ll x = 1) {
for (; b; b >>= 1, a = MUL(a, a))
if (b & 1) x = MUL(x, a);
return x;
}
auto inv = getInv(N); // NOLINT
/*---------------------------------------------------------------------------*/
namespace NTT {
const int g = 3; //原根
Poly Omega(int L) {
int wn = POW(g, P / L);
Poly w(L);
w[L >> 1] = 1;
fp(i, L / 2 + 1, L - 1) w[i] = MUL(w[i - 1], wn);
fd(i, L / 2 - 1, 1) w[i] = w[i << 1];
return w;
}
auto W = Omega(1 << 21); // NOLINT 初始化 2的次幂大于等于NTT长度
void DIF(int *a, int n) {
for (int k = n >> 1; k; k >>= 1)
for (int i = 0, y; i < n; i += k << 1)
for (int j = 0; j < k; ++j)
y = a[i + j + k], a[i + j + k] = MUL(a[i + j] - y + P, W[k +
j]), ADD(a[i + j], y);
}
void IDIT(int *a, int n) {
for (int k = 1; k < n; k <<= 1)
for (int i = 0, x, y; i < n; i += k << 1)
for (int j = 0; j < k; ++j)
x = a[i + j], y = MUL(a[i + j + k], W[k + j]),
a[i + j + k] = x - y < 0 ? x - y + P : x - y, ADD(a[i + j],
y);
int Inv = P - (P - 1) / n;
fp(i, 0, n - 1) a[i] = MUL(a[i], Inv);
reverse(a + 1, a + n);
}
}
/*---------------------------------------------------------------------------*/
namespace Polynomial {
// basic operator
int norm(int n) {
return 1 << (32 - __builtin_clz(n - 1));
}
void norm(Poly &a) {
if (!a.empty()) a.resize(norm(a.size()), 0);
else a =
{0};
}
void DFT(Poly &a) {
NTT::DIF(a.data(), a.size());
}
void IDFT(Poly &a) {
NTT::IDIT(a.data(), a.size());
}
Poly &dot(Poly &a, Poly &b) {
fp(i, 0, a.size() - 1) a[i] = MUL(a[i], b[i]);
return a;
}
// mul / div int
Poly &operator*=(Poly &a, int b) {
for (auto &x : a) x = MUL(x, b);
return
a;
}
Poly operator*(Poly a, int b) {
return a *= b;
}
Poly operator*(int a, Poly b) {
return b * a;
}
Poly &operator/=(Poly &a, int b) {
return a *= POW(b);
}
Poly operator/(Poly a, int b) {
return a /= b;
}
// Poly add / sub
Poly &operator+=(Poly &a, Poly b) {
a.resize(max(a.size(), b.size()));
fp(i, 0, b.size() - 1) ADD(a[i], b[i]);
return a;
}
Poly operator+(Poly a, Poly b) {
return a += b;
}
Poly &operator-=(Poly &a, Poly b) {
a.resize(max(a.size(), b.size()));
fp(i, 0, b.size() - 1) SUB(a[i], b[i]);
return a;
}
Poly operator-(Poly a, Poly b) {
return a -= b;
}
// Poly mul
Poly operator*(Poly a, Poly b) {
int n = a.size() + b.size() - 1, L = norm(n);
if (a.size() <= 8 || b.size() <= 8) {
Poly c(n);
fp(i, 0, a.size() - 1) fp(j, 0, b.size() - 1)
c[i + j] = (c[i + j] + (ll) a[i] * b[j]) % P;
return c;
}
a.resize(L), b.resize(L);
DFT(a), DFT(b), dot(a, b), IDFT(a);
return a.resize(n), a;
}
// Poly inv
Poly Inv2k(Poly a) { // |a| = 2 ^ k
int n = a.size(), m = n >> 1;
if (n == 1) return {POW(a[0])};
Poly b = Inv2k(Poly(a.begin(), a.begin() + m)), c = b;
b.resize(n), DFT(a), DFT(b), dot(a, b), IDFT(a);
fp(i, 0, n - 1) a[i] = i < m ? 0 : P - a[i];
DFT(a), dot(a, b), IDFT(a);
return move(c.begin(), c.end(), a.begin()), a;
}
Poly Inv(Poly a) {
int n = a.size();
norm(a), a = Inv2k(a);
return a.resize(n), a;
}
// Poly div / mod
Poly operator/(Poly a,Poly b) {
int k = a.size() - b.size() + 1;
if (k < 0) return {0};
reverse(a.begin(), a.end());
reverse(b.begin(), b.end());
b.resize(k), a = a * Inv(b);
a.resize(k), reverse(a.begin(), a.end());
return a;
}
pair<Poly, Poly> operator%(Poly a, const Poly& b) {
Poly c = a / b;
a -= b * c, a.resize(b.size() - 1);
return {c, a};
}
// Poly calculus
Poly deriv(Poly a) {
fp(i, 1, a.size() - 1) a[i - 1] = MUL(i, a[i]);
return a.pop_back(), a;
}
Poly integ(Poly a) {
a.push_back(0);
fd(i, a.size() - 1, 1) a[i] = MUL(inv[i], a[i - 1]);
return a[0] = 0, a;
}
// Poly ln
Poly Ln(Poly a) {
int n = a.size();
a = deriv(a) * Inv(a);
return a.resize(n - 1), integ(a);
}
// Poly exp
Poly Exp(Poly a) {
int n = a.size(), k = norm(n);
Poly b = {1}, c, d;
a.resize(k);
for (int L = 2; L <= k; L <<= 1) {
d = b, b.resize(L), c = Ln(b), c.resize(L);
fp(i, 0, L - 1) c[i] = a[i] - c[i] + (a[i] < c[i] ? P : 0);
ADD(c[0], 1), DFT(b), DFT(c), dot(b, c), IDFT(b);
move(d.begin(), d.end(), b.begin());
}
return b.resize(n), b;
}
// Poly sqrt
Poly Sqrt(Poly a) {
int n = a.size(), k = norm(n);
a.resize(k);
Poly b = {(new Cipolla(P))->sqrt(a[0]).first, 0}, c;
for (int L = 2; L <= k; L <<= 1) {
b.resize(L), c = Poly(a.begin(), a.begin() + L) * Inv2k(b);
fp(i, L / 2, L - 1) b[i] = MUL(c[i], (P + 1) / 2);
}
return b.resize(n), b;
}
// Poly pow
Poly Pow(Poly &a, int b) {
return Exp(Ln(a) * b); // a[0] = 1
}
Poly Pow(Poly a, int b1, int b2) { // b1 = b % P, b2 = b % phi(P) and b >= n iff a[0] > 0
int n = a.size(), d = 0, k;
while (d < n && !a[d]) ++d;
if ((ll) d * b1 >= n) return Poly(n);
a.erase(a.begin(), a.begin() + d);
k = POW(a[0]), norm(a *= k);
a = Pow(a, b1) * POW(k, P - 1 - b2);
a.resize(n), d *= b1;
fd(i, n - 1, 0) a[i] = i >= d ? a[i - d] : 0;
return a;
}
// Get [x ^ k](f / g)
int divAt(Poly f, Poly g, ll k) {
int n = max(f.size(), g.size()), m = norm(n);
for (; k; k >>= 1) {
f.resize(m * 2, 0), DFT(f);
g.resize(m * 2, 0), DFT(g);
fp(i, 0, 2 * m - 1) f[i] = MUL(f[i], g[i ^ 1]);
fp(i, 0, m - 1) g[i] = MUL(g[2 * i], g[2 * i + 1]);
g.resize(m), IDFT(f), IDFT(g);
for (int i = 0, j = k & 1; i < n; i++, j += 2) f[i] = f[j];
f.resize(n), g.resize(n);
}
return f[0];
}
// Get a[k] by a[n] = sum c[i] * a[n - i]
int LinearRecur(Poly a, Poly c, ll k) {
c[0] = P - 1, a = a * c, a.resize(c.size() - 1);
return divAt(a, c, k);
}
//Poly cyc循环卷积
Poly cyc(Poly a, Poly b) {
int n = a.size();
reverse(a.begin(), a.end());
b.insert(b.end(), b.begin(), b.end());
a = a * b;
fp(i, 0, n - 1) a[i] = a[i + n - 1];
return a.resize(n), a;
}
//Poly subpoly差卷积
Poly subpoly(Poly a, Poly b) {
int n = a.size();
int m = b.size();
reverse(b.begin(), b.end());
Poly c = a * b;
fp(i, 0, n - 1) c[i] = c[i + m - 1];
return c.resize(n), c;
}
}
/*---------------------------------------------------------------------------*/
Poly BerlekampMassey(Poly &a) {
Poly c, las, cur;
int k, delta, d, w;
fp(i, 0, a.size() - 1) {
d = P - a[i];
fp(j, 0, c.size() - 1) d = (d + (ll) a[i - j - 1] * c[j]) % P;
if (!d) continue;
if (c.empty()) {
k = i, delta = d, c.resize(i + 1);
continue;
}
cur = c, w = POW(delta, P - 2, P - d);
if (c.size() < las.size() + i - k) c.resize(las.size() + i - k);
SUB(c[i - k - 1], w);
fp(j, 0, las.size() - 1) c[i - k + j] = (c[i - k + j] + (ll) w * las[j])
% P;
if (cur.size() <= las.size() + i - k) las = cur, k = i, delta = d;
}
return c;
}
/*---------------------------------------------------------------------------*/
using namespace Polynomial;
signed main() {
Poly a(3);
a[0] = 1;
a[1] = 1;
a[2] = 1;
cout << (a * a)[3] << endl;
return 0;
}
- 分治ntt大板子(jls)
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
constexpr int mod = 998244353;
inline int read(){
int ch = getchar(), res = 0;
while(ch < '0' || ch > '9') ch = getchar();
while(ch >= '0' && ch <= '9'){
res = (res * 10ll + ch - '0') % mod;
ch = getchar();
}
return res;
}
int norm(int x) {
if (x < 0) {
x += mod;
}
if (x >= mod) {
x -= mod;
}
return x;
}
template<class T>
T power(T a, int b) {
T res = 1;
for (; b; b /= 2, a *= a) {
if (b % 2) {
res *= a;
}
}
return res;
}
struct Z { // 带模的整数
int x;
Z(int x = 0) : x(norm(x)) {}
int val() const {
return x;
}
Z operator-() const {
return Z(norm(mod - x));
}
Z inv() const {
assert(x != 0);
return power(*this, mod - 2);
}
Z &operator*=(const Z &rhs) {
x = i64(x) * rhs.x % mod;
return *this;
}
Z &operator+=(const Z &rhs) {
x = norm(x + rhs.x);
return *this;
}
Z &operator-=(const Z &rhs) {
x = norm(x - rhs.x);
return *this;
}
Z &operator/=(const Z &rhs) {
return *this *= rhs.inv();
}
friend Z operator*(const Z &lhs, const Z &rhs) {
Z res = lhs;
res *= rhs;
return res;
}
friend Z operator+(const Z &lhs, const Z &rhs) {
Z res = lhs;
res += rhs;
return res;
}
friend Z operator-(const Z &lhs, const Z &rhs) {
Z res = lhs;
res -= rhs;
return res;
}
friend Z operator/(const Z &lhs, const Z &rhs) {
Z res = lhs;
res /= rhs;
return res;
}
friend istream &operator>>(istream &is, Z &a) {
i64 v;
is >> v;
a = Z(v);
return is;
}
friend ostream &operator<<(ostream &os, const Z &a) {
return os << a.val();
}
};
vector<int> rev;
vector<Z> roots{0, 1};
void dft(vector<Z> &a) {
int n = a.size();
if (int(rev.size()) != n) {
int k = __builtin_ctz(n) - 1;
rev.resize(n);
for (int i = 0; i < n; i ++) {
rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
}
}
for (int i = 0; i < n; i ++) {
if (rev[i] < i) {
swap(a[i], a[rev[i]]);
}
}
if (int(roots.size()) < n) {
int k = __builtin_ctz(roots.size());
roots.resize(n);
while ((1 << k) < n) {
Z e = power(Z(3), (mod - 1) >> (k + 1));
for (int i = 1 << (k - 1); i < (1 << k); i ++) {
roots[i << 1] = roots[i];
roots[i << 1 | 1] = roots[i] * e;
}
k ++;
}
}
for (int k = 1; k < n; k *= 2) {
for (int i = 0; i < n; i += 2 * k) {
for (int j = 0; j < k; j ++) {
Z u = a[i + j], v = a[i + j + k] * roots[k + j];
a[i + j] = u + v, a[i + j + k] = u - v;
}
}
}
}
void idft(vector<Z> &a) {
int n = a.size();
reverse(a.begin() + 1, a.end());
dft(a);
Z inv = (1 - mod) / n;
for (int i = 0; i < n; i ++) {
a[i] *= inv;
}
}
struct Poly {
vector<Z> a;
Poly() {}
Poly(const vector<Z> &a) : a(a) {}
Poly(const initializer_list<Z> &a) : a(a) {}
int size() const {
return a.size();
}
void resize(int n) {
a.resize(n);
}
Z operator[](int idx) const {
if (idx < size()) {
return a[idx];
} else {
return 0;
}
}
Z &operator[](int idx) {
return a[idx];
}
Poly mulxk(int k) const {
auto b = a;
b.insert(b.begin(), k, 0);
return Poly(b);
}
Poly modxk(int k) const {
k = min(k, size());
return Poly(vector<Z>(a.begin(), a.begin() + k));
}
Poly divxk(int k) const {
if (size() <= k) {
return Poly();
}
return Poly(vector<Z>(a.begin() + k, a.end()));
}
friend Poly operator+(const Poly &a, const Poly &b) {
vector<Z> res(max(a.size(), b.size()));
for (int i = 0; i < int(res.size()); i ++) {
res[i] = a[i] + b[i];
}
return Poly(res);
}
friend Poly operator-(const Poly &a, const Poly &b) {
vector<Z> res(max(a.size(), b.size()));
for (int i = 0; i < int(res.size()); i ++) {
res[i] = a[i] - b[i];
}
return Poly(res);
}
friend Poly operator*(Poly a, Poly b) {
if (a.size() == 0 || b.size() == 0) {
return Poly();
}
int sz = 1, tot = min(5000000, a.size() + b.size() - 1);
while (sz < tot) {
sz *= 2;
}
a.a.resize(sz);
b.a.resize(sz);
dft(a.a);
dft(b.a);
for (int i = 0; i < sz; i ++) {
a.a[i] = a[i] * b[i];
}
idft(a.a);
a.resize(tot);
return a;
}
friend Poly operator*(Z a, Poly b) {
for (int i = 0; i < int(b.size()); i ++) {
b[i] *= a;
}
return b;
}
friend Poly operator*(Poly a, Z b) {
for (int i = 0; i < int(a.size()); i ++) {
a[i] *= b;
}
return a;
}
Poly &operator+=(Poly b) {
return (*this) = (*this) + b;
}
Poly &operator-=(Poly b) {
return (*this) = (*this) - b;
}
Poly &operator*=(Poly b) {
return (*this) = (*this) * b;
}
Poly deriv() const {
if (a.empty()) {
return Poly();
}
vector<Z> res(size() - 1);
for (int i = 0; i < size() - 1; i ++) {
res[i] = (i + 1) * a[i + 1];
}
return Poly(res);
}
Poly integr() const {
vector<Z> res(size() + 1);
for (int i = 0; i < size(); i ++) {
res[i + 1] = a[i] / (i + 1);
}
return Poly(res);
}
Poly inv(int m) const {
Poly x{a[0].inv()};
int k = 1;
while (k < m) {
k *= 2;
x = (x * (Poly{2} - modxk(k) * x)).modxk(k);
}
return x.modxk(m);
}
Poly log(int m) const {
return (deriv() * inv(m)).integr().modxk(m);
}
Poly exp(int m) const {
Poly x{1};
int k = 1;
while (k < m) {
k *= 2;
x = (x * (Poly{1} - x.log(k) + modxk(k))).modxk(k);
}
return x.modxk(m);
}
Poly pow(int k, int m) const {
int i = 0;
while (i < size() && a[i].val() == 0) {
i ++;
}
if (i == size() || 1LL * i * k >= m) {
return Poly(vector<Z>(m));
}
Z v = a[i];
auto f = divxk(i) * v.inv();
return (f.log(m - i * k) * k).exp(m - i * k).mulxk(i * k) * power(v, k);
}
Poly sqrt(int m) const {
Poly x{1};
int k = 1;
while (k < m) {
k *= 2;
x = (x + (modxk(k) * x.inv(k)).modxk(k)) * ((mod + 1) / 2);
}
return x.modxk(m);
}
Poly mulT(Poly b) const {
if (b.size() == 0) {
return Poly();
}
int n = b.size();
reverse(b.a.begin(), b.a.end());
return ((*this) * b).divxk(n - 1);
}
};
const int N = 1e6 + 10;
int n,k,P[N];
bool st[N];
int fac[N],inv[N];
vector<Z> f[N];
int qmi(int a,int b) {
int res = 1;
while(b) {
if(b & 1) res = 1ll * res * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return res;
}
void init() {
fac[0] = inv[0] = 1;
for(int i=1;i<N;i++) {
fac[i] = 1ll * fac[i - 1] * i % mod;
inv[i] = 1ll * inv[i - 1] * qmi(i,mod - 2) % mod;
}
}
int C(int a,int b) {
if(a < b) return 0;
return 1ll * fac[a] * inv[a - b] % mod * inv[b] % mod;
}
int get(int n,int m) {
if(n >= 2*m) {
return (C(n - m - 1,m - 1) + C(n - m,m) )% mod;
}
return 0;
}
Poly merge(int l,int r) {
if(l == r) return Poly(f[l]);
int mid = l + r >> 1;
return merge(l,mid) * merge(mid + 1,r);
}
void solve() {
scanf("%d%d",&n,&k);
for(int i=1;i<=n;i++) {
scanf("%d",&P[i]);
st[i] = false;
}
int tot = 0;
for(int i=1;i<=n;i++) {
if(st[i]) continue;
int cnt = 0;
tot ++ ;
int x = i;
while(!st[x]) {
cnt ++ ;
st[x] = true;
x = P[x];
}
f[tot].resize(cnt / 2 + 1);
f[tot][0] = 1;
for(int j=1;j<=cnt/2;j++) {
{
f[tot][j] = get(cnt,j);
}
}
}
Poly ans = merge(1,tot);
ans.resize(k + 1);
printf("%d\n",ans[k]);
}
int main() {
// freopen("data.in","r",stdin);
init();
int T;
scanf("%d",&T);
while (T --) {
solve();
}
return 0;
}
线性基(求最大异或和)
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 100010;
int n;
LL a[N];
int main()
{
scanf("%d", &n);
for (int i = 0; i < n; i ++ ) scanf("%lld", &a[i]);
int k = 0;
for (int i = 62; i >= 0; i -- )
{
for (int j = k; j < n; j ++ )
if (a[j] >> i & 1)
{
swap(a[j], a[k]);
break;
}
if (!(a[k] >> i & 1)) continue;
for (int j = 0; j < n; j ++ )
if (j != k && (a[j] >> i & 1))
a[j] ^= a[k];
k ++ ;
if (k == n) break;
}
LL res = 0;
for (int i = 0; i < k; i ++ ) res ^= a[i];
printf("%lld\n", res);
return 0;
}
MIN25求1-n的质数和
#include <bits/stdc++.h>
using namespace std;
const int N = 1000010;
typedef long long LL;
namespace Min25 {
int prime[N], id1[N], id2[N], flag[N], ncnt, m;
LL g[N], sum[N], a[N], T, n;
inline int ID(LL x) {
return x <= T ? id1[x] : id2[n / x];
}
inline LL calc(LL x) {
return x * (x + 1) / 2 - 1;
}
inline LL f(LL x) {
return x;
}
inline void init() {
m = ncnt = 0;
T = sqrt(n + 0.5);
for (int i = 2; i <= T; i++) {
if (!flag[i]) prime[++ncnt] = i, sum[ncnt] = sum[ncnt - 1] + i;
for (int j = 1; j <= ncnt && i * prime[j] <= T; j++) {
flag[i * prime[j]] = 1;
if (i % prime[j] == 0) break;
}
}
for (LL l = 1; l <= n; l = n / (n / l) + 1) {
a[++m] = n / l;
if (a[m] <= T) id1[a[m]] = m; else id2[n / a[m]] = m;
g[m] = calc(a[m]);
}
for (int i = 1; i <= ncnt; i++)
for (int j = 1; j <= m && (LL)prime[i] * prime[i] <= a[j]; j++)
g[j] = g[j] - (LL)prime[i] * (g[ID(a[j] / prime[i])] - sum[i - 1]);
}
inline LL solve(LL x) {
if (x <= 1) return x;
return n = x, init(), g[ID(n)];
}
}
int main() {
LL n; scanf("%lld", &n);
printf("%lld\n", Min25::solve(n));
}
数学小结论
- (a + b) = (a ^ b) + 2 * (a & b)
- $\lfloor x / y \rfloor$ = z,已知x和z求y的范围是[x / (z + 1) + 1,x / z]
- 网格中从(0,0)出发向下和向右走,到达(x,y)的方案数未C(x,x + y)
- 组合数求和:$ \sum_{i=0}^k\binom{n+i}n=\binom{n+k+1}{n+1} $
- 错排公式: D(n) = (n-1) [D(n-2) + D(n-1)],特殊地,D(1) = 0, D(2) = 1.
- 博弈论公平组合游戏可以视为一颗有向数,一个节点的sg(u) = 子节点mex{sg(v)},可以将多个状态拆分成,多个局面下的一个状态,并且sg(root_1) ^ sg(root_2) …。
- n 的链选出 m 个互不相邻的点的方案数$ \binom{n-m+1}{m} $ (n >= 2 * m - 1)
- n 的环选出 m 个互不相邻的点的方案数$ \binom{n-m-1}{m - 1} + \binom{n-m}{m}$ * (n >= 2 * m)