Number 数论
作者：互联网
Number 数论
素数筛
欧拉筛
对if(i % prime[j] == 0) break;的解释
当i % prime[j] == 0时有 i = k * prime[j]; 若j++有 i * prime[j + 1] = k * prime[j] * prime[j + 1] 也是prime[j]的因子，导致重复筛
const int maxn = 1e7 + 8;
vector<int > prime(maxn, 0);
bitset<maxn> vis(0);
void getPrime() {
for(int i = 2;i < maxn;i++) {
if(vis[i] == 0) prime[++prime[0]] = i;
for(int j = 1;j <= prime[0];j++) {
if(i * prime[j] >= maxn) break;
vis[i * prime[j]] = 1;
if(i % prime[j] == 0) break;
}
}
}
埃式筛
void Prime(){
vector<int > vis(maxn, 0); //初始化都是素数
vis[0] = vis[1] = 1; //0 和 1不是素数
for (int i = 2; i <= maxn; i++) {
if (!vis[i]) { //如果i是素数，让i的所有倍数都不是素数
for (int j = i * i; j <= maxn; j += i) {
vis[j] = 1;
}
}
}
}
欧拉函数
用处：1：求1-n内与n互质的数量
2：求合数的逆元
3：欧拉降幂
int phi(int x) {
int ans = x;
for(int i = 2;i * i <= x;i++) {
if(x % i == 0) {
ans = ans / i * (i - 1);
while(x % i == 0) x /= i;
}
}
if(x > 1) ans = ans / x * (x - 1);
return ans;
}
线性筛欧拉函数
const int maxn = 100005;
vector<int > prime(maxn, 0), vis(maxn, 0), phi(maxn, 0);
void getPrime() {
phi[1] = 1;
for(int i = 2;i < maxn;i++) {
if(vis[i] == 0) prime[++prime[0]] = i, phi[i] = i - 1;
for(int j = 1;j <= prime[0];j++) {
if(i * prime[j] >= maxn) break;
vis[i * prime[j]] = 1;
if(i % prime[j] != 0) phi[i * prime[j]] = phi[i] * (prime[j] - 1);
//对于任意的a，b若gcd(a, b) == 1 -> phi(a * b) = phi(a) * phi(b)
else {
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
}
}
}
任意区间素数打表
const int M=1e6+5,N=(1<<16);
int prime[N+5],is[N+5],tot;
void getPrime(){
is[1]=1;//mmpaaaaaaaaaaaaaaaa
for(int i=2;i<=N;++i){
if(!is[i]) prime[++tot]=i;
for(int j=1;j<=tot&&(ll)i*prime[j]<=N;++j){
is[i*prime[j]]=1;
if(i%prime[j]==0) break;
}
}
}
int issp[M];
int main(){
getPrime();
int T,cas=0;scanf("%d",&T);
while(T--){
m(issp,0);
int a,b;scanf("%d%d",&a,&b);
if(a<=N&&b<=N) {
int ans=0;
for(int i=a;i<=b;++i) if(!is[i]) ++ans;
printf("Case %d: %d\n",++cas,ans);
continue;
}
int ans=0;
if(a<=N) {
for(int i=a;i<=N;++i) if(!is[i]) ++ans;
a=N+1;
}
for(int i=1;i<=tot;++i) {
ll l=a/prime[i]*prime[i];//左端点l.
if(l<a) l+=prime[i];
if(l==prime[i]) l+=prime[i];
for(ll j=l;j<=b;j+=prime[i]) issp[j-a]=1;
}
for(int j=a;j<=b;++j) if(!issp[j-a]) ++ans;
printf("Case %d: %d\n",++cas,ans);
}
}
大数素数测试
Miller_Rabin
#include <bits/stdc++.h>
using namespace std;
using ll = __int128;
const int S = 8; //随机算法判定次数一般 8～10 就够了
// 计算 ret = (a*b)%c a,b,c < 2^63
__int128 mult_mod(__int128 a, __int128 b, __int128 c) {
a %= c;
b %= c;
__int128 ret = 0;
__int128 tmp = a;
while (b) {
if (b & 1) {
ret += tmp;
if (ret > c)
ret -= c;//直接取模慢很多
}
tmp <<= 1;
if (tmp > c)
tmp -= c;
b >>= 1;
}
return ret;
}
// 计算 ret = (a^n)%mod
__int128 pow_mod(__int128 a, __int128 n, __int128 mod) {
__int128 ret = 1;
__int128 temp = a % mod;
while (n) {
if (n & 1)
ret = mult_mod(ret, temp, mod);
temp = mult_mod(temp, temp, mod);
n >>= 1;
}
return ret;
}
bool check(__int128 a, __int128 n, __int128 x, __int128 t) {
__int128 ret = pow_mod(a, x, n);
__int128 last = ret;
for (int i = 1; i <= t; i++) {
ret = mult_mod(ret, ret, n);
if (ret == 1 && last != 1 && last != n - 1)
return true;//合数
last = ret;
}
if (ret != 1)
return true;
else
return false;
}
//**************************************************
// Miller_Rabin 算法
// 是素数返回 true,(可能是伪素数)
// 不是素数返回 false
//**************************************************
bool MR(__int128 n) {
if (n < 2)
return false;
if (n == 2)
return true;
if ((n & 1) == 0)
return false;//偶数
__int128 x = n - 1;
__int128 t = 0;
while ((x & 1) == 0) {
x >>= 1;
t++;
}
srand(time(NULL));
for (int i = 0; i < S; i++) {
__int128 a = rand() % (n - 1) + 1;
if (check(a, n, x, t))
return false;
}
return true;
}
inline __int128 read() {
__int128 x = 0, f = 1;
char ch = getchar();
while (ch < '0' || ch > '9') {
if (ch == '-')
f = -1;
ch = getchar();
}
while (ch >= '0' && ch <= '9') {
x = x * 10 + ch - '0';
ch = getchar();
}
return x * f;
}
int main() {
__int128 n = read();
if (MR(n))
puts("Yes");
else
puts("No");
return 0;
}
Miller_Rabin(Easy ver)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll mul(ll a, ll b, ll p) {
return (__int128)a * b % p;
}
ll qp(ll x, ll y, ll p) {
ll z = 1;
for (; y; y >>= 1, x = mul(x, x, p)) {
if (y & 1) z = mul(z, x, p);
}
return z;
}
bool MR(ll x, ll b) {
for (ll k = x - 1; k; k >>= 1) {
ll cur = qp(b, k, x);
if (cur != 1 && cur != x - 1) return 0;
if (k & 1 || cur == x - 1) return 1;
}
return true;
}
bool test(ll x) {
if (x == 1) return 0;
static ll p[] = {2, 3, 5, 7, 17, 19, 61};
for (ll y : p) {
if (x == y) return 1;
if (!MR(x, y)) return 0;
}
return 1;
}
int main() {
ll n;
while (scanf("%lld", &n) != EOF) {
puts(test(n) ? "Y" : "N");
}
return 0;
}
大数分解
逆元
Exgcd求a在b意义下的逆元
i f ( g c d ( a , b ) = = 1 ) if(gcd(a, b)==1) if(gcd(a,b)==1)
int exgcd(int a, int b, int &x, int &y) {
if(b == 0) {
x = 1;
y = 0;
return a;
}
else {
int gcd = exgcd(b, a % b, x, y);
int t = x;
x = y;
y = t - a / b * y;
return gcd;
}
}
int getInv(int a, int b){
int x, y;
exgcd(a, b, x, y);
return (x + b) % b;
}
费马小定理求a在mod p意义下的逆元(p是素数)
a ∗ a p − 2 = 1 ( m o d p ) a * a^{p - 2} = 1{(modp)} a∗ap−2=1(modp)
int qpow(int a, int b, int mod) {
int ans = 1;
while(b) {
if(b & 1) ans = ans * a % mod;
b >>= 1, a = a * a % mod;
}
return ans;
}
int getInv(int x, int mod) {
return qpow(x, mod - 2, mod);
}
费马小定理求a在mod p意义下的逆元(p不是素数)
a ∗ a p h i ( p ) − 1 = 1 ( m o d p ) a*a^{phi(p) - 1} = 1(modp) a∗aphi(p)−1=1(modp)
int qpow(int a, int b, int mod) {
int ans = 1;
while(b) {
if(b & 1) ans = ans * a % mod;
b >>= 1, a = a * a % mod;
}
return ans;
}
int phi(int x) {
int ans = x;
for(int i = 2;i * i <= x;i++) {
if(x % i == 0) {
ans = ans * (i - 1) / i;
while(x % i == 0) x /= i;
}
}
if(x > 1) ans = ans * (x - 1) / x;
return ans;
}
int getInv(int x, int mod) {
return qpow(x, phi(mod) - 1, mod);
}
O(n)求逆元
vector<ll> inv(3000006, 0);
void liner_inv(int n, int p) {
inv[1] = 1;
for(int i = 2;i <= n;i++) inv[i] = 1ll * (p - p / i) * inv[p % i] % p;
}
O(n + log(mod))求任意n个数的逆元
const ll mod = 1e9 + 7;
const ll p = 998244353;
const ll maxn = 1e5 + 7;
vector<ll> a(maxn, 0), pre(maxn, 1), pre_inv(maxn, 0), inv(maxn, 0);
// 数组元素 前缀积 前缀积的逆元 每个数的逆元
ll qpow(ll a, ll b) {
ll ret = 1;
while(b) {
if(b & 1) ret = ret * a % mod;
b >>= 1, a = a * a % mod;
}
return ret;
}
void solve(int n) {//数组长度参数
for(int i = 1;i <= n;i++) cin >> a[i];
for(int i = 1;i <= n;i++) pre[i] = pre[i - 1] * a[i] % mod;
pre_inv[n] = qpow(pre[n], mod - 2);
for(int i = n;i >= 1;i--) pre_inv[i - 1] = pre_inv[i] * a[i] % mod;
for(int i = 1;i <= n;i++) inv[i] = pre_inv[i] * pre[i - 1] % mod;
}
小数循环节
int c(int tmp) {
while(tmp % 2 == 0) tmp /= 2;
while(tmp % 5 == 0) tmp /= 5;
int len = 0, e = 1;
if(tmp == 1) return -1;
while(1) {
e = e * 10 % tmp;
if(e == 1) break;
len++;
}
return len;
}
组合数
C n m = n ! m ! ∗ ( n − m ) ! C_n^m = \frac{n!}{m!*(n - m)!} Cnm=m!∗(n−m)!n!
正常版本
const int maxn = 2e5 + 7;
vector<ll> fac(maxn, 1);
void init(ll p) {//if module is different
for(ll i = 1;i <= 200000;i++) fac[i] = fac[i - 1] * i % p;
}
ll qpow(ll a, ll b, ll p) {
ll ret = 1;
while(b) {
if(b & 1) ret = ret * a % p;
b >>= 1;
a = a * a % p;
}
return ret;
}
ll getInv(ll x, ll p) {
if(x == 1) return 1;
return qpow(x, p - 2, p);
}
ll c(ll n, ll m, ll p) {
if(!m) return 1;
if(m > n) return 0;
return fac[n] * getInv(fac[m], p) % p * getInv(fac[n - m], p) % p;
}
Mint版本
constexpr int mod = 1000000007;
constexpr int maxn = 1e5 + 6;
// assume -mod <= x < 2mod
int norm(int x) {if(x < 0) x += mod; if (x >= mod) x -= mod; return x; }
template<class T>T power(T a, ll b){T res = 1;for (; b; b /= 2, a *= a)if (b & 1)res *= a;return res;}
struct Mint {
int x;Mint(int x = 0) : x(norm(x)){}
int val() const {return x;}
Mint operator-() const {return Mint(norm(mod - x));}
Mint inv() const {assert(x != 0);return power(*this, mod - 2);}
Mint &operator*=(const Mint &rhs) { x = ll(x) * rhs.x % mod; return *this;}
Mint &operator+=(const Mint &rhs) { x = norm(x + rhs.x); return *this;}
Mint &operator-=(const Mint &rhs) { x = norm(x - rhs.x); return *this;}
Mint &operator/=(const Mint &rhs) {return *this *= rhs.inv();}
friend Mint operator*(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res *= rhs; return res;}
friend Mint operator+(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res += rhs; return res;}
friend Mint operator-(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res -= rhs; return res;}
friend Mint operator/(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res /= rhs; return res;}
};
vector<Mint> fac(maxn, 1);
void init() {
for(ll i = 1;i < maxn;i++) fac[i] = fac[i - 1] * i;
}
Mint c(ll n, ll m) {
return fac[n] / fac[m] / fac[n - m];
}
容斥原理
CRT
using ll = long long;
ll qpow(ll a, ll b, ll mod) {
ll ret = 1;
while(b) {
if(b & 1) ret = ret * a % mod;
b >>= 1;
a = a * a % mod;
}
return ret;
}
ll getInv(ll x, ll mod) {
return qpow(x, mod - 2, mod);
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
ll n, mul = 1, sum = 0;
cin >> n;
vector<ll> a(n + 1, 0), mu(n + 1, 0);
for(int i = 1;i <= n;i++) cin >> a[i] >> mu[i];
for(int i = 1;i <= n;i++) mul *= a[i];
for(int i = 1;i <= n;i++) {
ll lcm = mul / a[i];
ll val = getInv(lcm, a[i]) * mu[i] * lcm;
sum += val;
}
cout << sum % mul << '\n';
return 0;
}
Lucas/ExLucas
Lucas
const int maxn = 2e5 + 7;
vector<ll> fac(maxn, 1);
void init(ll p) {//if module is different
for(ll i = 1;i <= 200000;i++) fac[i] = fac[i - 1] * i % p;
}
ll qpow(ll a, ll b, ll p) {
ll ret = 1;
while(b) {
if(b & 1) ret = ret * a % p;
b >>= 1;
a = a * a % p;
}
return ret;
}
ll getInv(ll x, ll p) {
if(x == 1) return 1;
return qpow(x, p - 2, p);
}
ll c(ll n, ll m, ll p) {
if(!m) return 1;
if(m > n) return 0;
return fac[n] * getInv(fac[m], p) % p * getInv(fac[n - m], p) % p;
}
ll lucas(ll n, ll m, ll p) {
if(n == 0) return 1;
return lucas(n / p, m / p, p) * c(n % p, m % p, p) % p;
}
Mint version
constexpr int mod = 1000000007;
const int maxn = 1e5 + 6;
// assume -mod <= x < 2mod
int norm(int x) {if(x < 0) x += mod; if (x >= mod) x -= mod; return x; }
template<class T>T power(T a, ll b){T res = 1;for (; b; b /= 2, a *= a)if (b & 1)res *= a;return res;}
struct Mint {
int x;Mint(int x = 0) : x(norm(x)){}
int val() const {return x;}
Mint operator-() const {return Mint(norm(mod - x));}
Mint inv() const {assert(x != 0);return power(*this, mod - 2);}
Mint &operator*=(const Mint &rhs) { x = ll(x) * rhs.x % mod; return *this;}
Mint &operator+=(const Mint &rhs) { x = norm(x + rhs.x); return *this;}
Mint &operator-=(const Mint &rhs) { x = norm(x - rhs.x); return *this;}
Mint &operator/=(const Mint &rhs) {return *this *= rhs.inv();}
friend Mint operator*(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res *= rhs; return res;}
friend Mint operator+(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res += rhs; return res;}
friend Mint operator-(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res -= rhs; return res;}
friend Mint operator/(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res /= rhs; return res;}
};
vector<Mint> fac(maxn, 1);
void init() {
for(ll i = 1;i < maxn;i++) fac[i] = fac[i - 1] * i;
}
Mint c(ll n, ll m) {
return fac[n] / fac[m] / fac[n - m];
}
Mint lucas(ll n, ll m) {
if(m == 0) return 1;
if(n == 0) return 1;
if(n < mod && m < mod) return c(n, m);
return lucas(n / mod, m / mod) * c(n % mod, m % mod);
}
Gcd
Exgcd
int exgcd(int a, int b, int &x, int &y) {
if(b == 0) {
x = 1;
y = 0;
return a;
}
else {
int gcd = exgcd(b, a % b, x, y);
int t = x;
x = y;
y = t - a / b * y;
return gcd;
}
}
矩阵快速幂
实例：求斐波那契数列的第n项
$$
\begin{bmatrix}
{}&{}&{}&{}\
{}&{}&{}&{}\
{}&{}&{}&{}\
{}&{}&{}&{}\
\end{bmatrix}
\begin{bmatrix}
{}&{}&{}\
{}&{}&{}\
{}&{}&{}\
\end{bmatrix}
\begin{bmatrix}
{}&{}\
{}&{}\
\end{bmatrix}
$$
$$
\begin{bmatrix}
{f_{}}&{f_{}}&{f_{}}&{f_{}}\
{f_{}}&{f_{}}&{f_{}}&{f_{}}\
{f_{}}&{f_{}}&{f_{}}&{f_{}}\
{f_{}}&{f_{}}&{f_{}}&{f_{}}\
\end{bmatrix}
\begin{bmatrix}
{f_{}}&{f_{}}&{f_{}}\
{f_{}}&{f_{}}&{f_{}}\
{f_{}}&{f_{}}&{f_{}}\
\end{bmatrix}
\begin{bmatrix}
{f_{}}&{f_{}}\
{f_{}}&{f_{}}\
\end{bmatrix}
$$
$$
\begin{bmatrix}
{f_{}}\
{f_{}}\
{f_{}}\
{f_{}}\
\end{bmatrix}
\begin{bmatrix}
{f_{}}\
{f_{}}\
{f_{}}\
\end{bmatrix}
\begin{bmatrix}
{f_{}}\
{f_{}}\
\end{bmatrix}
$$
$$
f_{i} = a*f_{i - 1} + i:
\begin{bmatrix}
{a}&{1}&{0}\
{0}&{1}&{1}\
{0}&{0}&{1}\
\end{bmatrix}
\cdot
\begin{bmatrix}
{f_{i-1}}\
{{i}}\
{{1}}\
\end{bmatrix}
\begin{bmatrix}
{f_{i}}\
{{i+1}}\
{{1}}\
\end{bmatrix}
$$
$$
f_{i} = f_{i - 1} + f_{i - 2}:
\begin{bmatrix}
{1}&{1}\
{1}&{0}\
\end{bmatrix}
\cdot
\begin{bmatrix}
{f_{n-1}}\
{f_{n-2}}\
\end{bmatrix}
\begin{bmatrix}
{f_{n}}\
{f_{n-1}}\
\end{bmatrix}
$$
#include<bits/stdc++.h>
using namespace std;
using ll = long long;
int n_; //Size of Matrix
const ll p = 1e9 + 7;
const int N = 10;
struct matrix {
ll a[N][N]{};
matrix init() { //to E
matrix E;
for(register int i = 0;i < n_;i++) E.a[i][i] = 1;
return E;
}
matrix f1_() { //f(n) = f(n - 1) + f(n - 2), f[n] = qpow(A, n - 2);
matrix ret;
ret.a[0][0] = 1, ret.a[0][1] = 1;
ret.a[1][0] = 1, ret.a[1][1] = 0;
return ret;
}
matrix f2_() { //f(n) = f(n - 1) + f(n - 2) + f(n - 3);
matrix ret;
ret.a[0][0] = 1, ret.a[0][1] = 1, ret.a[0][2] = 1;
ret.a[1][0] = 1, ret.a[1][1] = 0, ret.a[1][2] = 0;
ret.a[2][0] = 0, ret.a[2][1] = 1, ret.a[2][2] = 0;
return ret;
}
matrix f3_() {
matrix ret;
//to construct ...
return ret;
}
matrix operator * (const matrix& b) {
matrix ret;
for(register int i = 0;i < n_;i++)
for(register int j = 0;j < n_;j++)
for(register int k = 0;k < n_;k++)
ret.a[i][j] = (ret.a[i][j] + a[i][k] * b.a[k][j]) % p;
return ret;
}
matrix qpow(matrix a, ll b) {
matrix ans = init();
while(b) {
if(b & 1) ans = ans * a;
b >>= 1, a = a * a;
}
return ans;
}
void show() {
for(register int i = 0;i < n_;i++)
for(register int j = 0;j < n_;j++)
cout << a[i][j] << " \n"[j == n_ - 1];
return ;
}
};
ll qpow_(ll a, ll b, ll mod) {
ll ret = 1;
while(b){
if(b & 1) ret = ret * a % mod;
b >>= 1, a = a * a % mod;
}
return ret;
}
void Fib_n() {
int n; cin >> n;
if(n == 1 || n == 2) {
cout << 1 << '\n';
return ;
}
matrix f = f.f1_();
n_ = 2;
f = f.qpow(f, n - 2);
cout << f.a[0][0] + f.a[0][1] << '\n';
}
int main() {
Fib_n();
return 0;
}
约瑟夫环
FFT
const int MAXN = 1e6 + 10;
const double Pi = acos(-1.0);
struct complex{
double x, y;
complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
}a[MAXN], b[MAXN];
complex operator + (complex a, complex b) { return complex(a.x + b.x , a.y + b.y);}
complex operator - (complex a, complex b) { return complex(a.x - b.x , a.y - b.y);}
complex operator * (complex a, complex b) { return complex(a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x);}
int l, r[MAXN], ans[MAXN];
int limit = 1;
void fft(complex *A, int type) {
for(int i = 0; i < limit; i++) if (i < r[i]) swap(A[i], A[r[i]]);
for (int mid = 1; mid < limit; mid <<= 1) {
complex Wn( cos(Pi / mid) , type * sin(Pi / mid) );
for (int R = mid << 1, j = 0; j < limit; j += R) {
complex w(1, 0);
for (int k = 0; k < mid; k++, w = w * Wn) {
complex x = A[j + k], y = w * A[j + mid + k];
A[j + k] = x + y;
A[j + mid + k] = x - y;
}
}
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
string aa, bb;
cin >> aa >> bb;
int n = aa.size(), m = bb.size();
for (int i = 0; i < n; i++) a[i].x = aa[i] ^ 0x30;
for (int i = 0; i < m; i++) b[i].x = bb[i] ^ 0x30;
while (limit <= n + m) limit <<= 1, l++;
for (int i = 0; i < limit; i++) r[i] = ( r[i >> 1] >> 1 ) | ( (i & 1) << (l - 1));
fft(a, 1);
fft(b, 1);
for (int i = 0; i <= limit; i++) a[i] = a[i] * b[i];
fft(a, -1);
for (int i = 0; i <= n + m; i++) ans[i] = (int)(a[i].x / limit + 0.5);
for(int i = n + m;i >= 1;i--) ans[i - 1] += ans[i] / 10, ans[i] %= 10;
for(int i = 0;i <= n + m - 2;i++) cout << ans[i]; cout << '\n';
}
NTT
#include <bits/stdc++.h>
using namespace std;
constexpr int P = 998244353;
vector<int> rev, roots{0, 1};
int power(int a, int b) {
int res = 1;
for(; b; b >>= 1, a = 1ll * a * a % P) if(b & 1) res = 1ll * res * a % P;
return res;
}
void dft(vector<int> &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) {
int e = power(3, (P - 1) >> (k + 1));
for(int i = 1 << (k - 1); i < (1 << k); ++i) {
roots[2 * i] = roots[i];
roots[2 * i + 1] = 1ll * roots[i] * e % P;
}
++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) {
int u = a[i + j];
int v = 1ll * a[i + j + k] * roots[k + j] % P;
int x = u + v;
if(x >= P) x -= P;
a[i + j] = x;
x = u - v;
if(x < 0) x += P;
a[i + j + k] = x;
}
}
}
}
void idft(vector<int> &a) {
int n = a.size();
reverse(a.begin() + 1, a.end());
dft(a);
int inv = power(n, P - 2);
for(int i = 0; i < n; ++i) a[i] = 1ll * a[i] * inv % P;
}
struct Poly {
vector<int> a;
Poly(){}
Poly(int a0) { if (a0) a = {a0}; }
Poly(const vector<int> &a1) : a(a1) {
while(!a.empty() && !a.back()) a.pop_back();
}
int size() const { return a.size(); }
int operator[](int idx) const { if (idx < 0 || idx >= size()) return 0; 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<int>(a.begin(), a.begin() + k));
}
Poly divxk(int k) const {
if (size() <= k) return Poly();
return Poly(vector<int>(a.begin() + k, a.end()));
}
friend Poly operator+(const Poly a, const Poly &b) {
vector<int> res(max(a.size(), b.size()));
for (int i = 0; i < int(res.size()); ++i) {
res[i] = a[i] + b[i];
if (res[i] >= P) res[i] -= P;
}
return Poly(res);
}
friend Poly operator-(const Poly a, const Poly &b) {
vector<int> res(max(a.size(), b.size()));
for (int i = 0; i < int(res.size()); ++i) {
res[i] = a[i] - b[i];
if (res[i] < 0) res[i] += P;
}
return Poly(res);
}
friend Poly operator*(Poly a, Poly b) {
int sz = 1, tot = 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] = 1ll * a[i] * b[i] % P;
idft(a.a);
return Poly(a.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<int> res(size() - 1);
for (int i = 0; i < size() - 1; ++i) res[i] = 1ll * (i + 1) * a[i + 1] % P;
return Poly(res);
}
Poly integr() const {
if(a.empty()) return Poly();
vector<int> res(size() + 1);
for (int i = 0; i < size(); ++i) res[i + 1] = 1ll * a[i] * power(i + 1, P - 2) % P;
return Poly(res);
}
Poly inv(int m) const {
Poly x(power(a[0], P - 2));
int k = 1;
while(k < m) {
k *= 2;
x = (x * (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 * (1 - x.log(k) + modxk(k))).modxk(k);
}
return x.modxk(m);
}
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)) * ((P + 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);
}
vector<int> eval(vector<int> x) const {
if (size() == 0) return vector<int>(x.size(), 0);
const int n = max(int(x.size()), size());
vector<Poly> q(4 * n);
vector<int> ans(x.size());
x.resize(n);
function<void(int, int, int)> build = [&](int p, int l, int r) {
if (r - l == 1) {
q[p] = vector<int>{1, (P - x[l]) % P};
}
else {
int m = (l + r) / 2;
build(2 * p, l, m);
build(2 * p + 1, m, r);
q[p] = q[2 * p] * q[2 * p + 1];
}
};
build(1, 0, n);
function<void(int, int, int, const Poly &)> work = [&](int p, int l, int r, const Poly &num) {
if (r - l == 1) {
if (l < int(ans.size()))
ans[l] = num[0];
}
else {
int m = (l + r) / 2;
work(2 * p, l, m, num.mulT(q[2 * p + 1]).modxk(m - l));
work(2 * p + 1, m, r, num.mulT(q[2 * p]).modxk(r - m));
}
};
work(1, 0, n, mulT(q[1].inv(n)));
return ans;
}
};
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
int n;
cin >> n;
vector<int> f;
for (int i = 0; i < n; i++) {
int a;
cin >> a;
if (a >= int(f.size())) {
f.resize(a + 1);
}
f[a] = 1;
}
int mx = f.size() - 1;
auto g = f;
reverse(g.begin(), g.end()); //位置取反 //差的卷积
auto h = Poly(f) * Poly(g); //h 数组从mx开始为差值为0，后边以此类推
return 0;
}
FWT
OR
void fwt_or(vector<ll> &a, int len) {
for (int mid = 2; mid <= len; mid <<= 1) {
for (int i = 0; i < len; i += mid) {
for (int j = i; j < i + (mid >> 1); j++) {
a[j + (mid >> 1)] += a[j];
}
}
}
return ;
}
void ifwt_or(vector<ll> &a, int len) {
for (int mid = 2; mid <= len; mid <<= 1) {
for (int i = 0; i < len; i += mid) {
for (int j = i; j < i + (mid >> 1); j++) {
a[j + (mid >> 1)] -= a[j];
}
}
}
return ;
}
AND
void fwt_and(vector<ll> &a, int len) {
for (int mid = 2; mid <= len; mid <<= 1) {
for (int i = 0; i < len; i += mid) {
for (int j = i; j < i + (mid >> 1); j++) {
a[j] += a[j + (mid >> 1)];
}
}
}
return ;
}
void ifwt_and(vector<ll> &a, int len) {
for (int mid = 2; mid <= len; mid <<= 1) {
for (int i = 0; i < len; i += mid) {
for (int j = i; j < i + (mid >> 1); j++) {
a[j] -= a[j + (mid >> 1)];
}
}
}
return ;
}
XOR
void fwt_xor(vector<ll> &a, int len) {
for (int mid = 2; mid <= len; mid <<= 1) {
for (int i = 0; i < len; i += mid) {
for (int j = i; j < i + (mid >> 1); j++) {
ll x = a[j], y = a[j + (mid >> 1)];
a[j] = x + y, a[j + (mid >> 1)] = x - y;
}
}
}
return ;
}
void ifwt_xor(vector<ll> &a, int len) {
for (int mid = 2; mid <= len; mid <<= 1) {
for (int i = 0; i < len; i += mid) {
for (int j = i; j < i + (mid >> 1); j++) {
ll x = a[j], y = a[j + (mid >> 1)];
a[j] = (x + y) >> 1, a[j + (mid >> 1)] = (x - y) >> 1;
}
}
}
return ;
}
标签：return,数论,ll,Number,int,const,Mint,mod 来源： https://blog.csdn.net/u014711890/article/details/120841840