數學
逆元
離線法
class C_INV
{
private:
int Pow(int a, lnt b, int P)
{
int res = 1;
while (b)
{
if (b & 1)
{
res = (1ll * res * a) % P;
}
b >>= 1, a = (1ll * a * a) % P;
}
return res;
}
public:
vector<int> operator()(vector<int>& arr, int P)
{
vector<int>sum = arr;
for (int i = 1; i < sum.size(); ++i)
{
sum[i] = 1ll * sum[i - 1] * arr[i] % P;
}
vector<int>inv = sum;
inv.back() = Pow(inv.back(), P - 2, P);
for (int i = inv.size() - 2; i >= 0; --i)
{
inv[i] = 1ll * inv[i + 1] * arr[i + 1] % P;
}
for (int i = 1; i < inv.size(); ++i)
{
inv[i] = 1ll * sum[i - 1] * inv[i] % P;
}
return inv;
}
};
快速冪法
class C_INV
{
private:
int P;
/*====================*/
int Pow(int a, lnt b, int P)
{
int res = 1;
while (b)
{
if (b & 1)
{
res = (1ll * res * a) % P;
}
b >>= 1, a = (1ll * a * a) % P;
}
return res;
}
public:
void Init(int P)
{
this->P = P;
}
int operator[](int x)
{
return Pow(x % P, P - 2, P);
}
};
Farey序列
class C_INV
{
private:
struct Frac
{
int up, dw;
Frac(int _up = 0, int _dw = 0)
{
up = _up, dw = _dw;
}
};
/*====================*/
int P, B1, B2;
vector<Frac>s;
vector<int>inv;
vector<int>pre, suf;
/*====================*/
int Calc(int x, Frac frac)
{
lnt pos = 1ll * x * frac.dw;
if (abs(pos - 1ll * P * frac.up) > P / B1) return -1;
if ((pos %= P) <= P / B1) return 1ll * frac.dw * inv[pos] % P;
return P - 1ll * frac.dw * inv[P - pos] % P;
}
public:
void Init(int P)
{
this->P = P;
B1 = pow(P, 1.0 / 3); B2 = B1 * B1;
s.assign(B2 + 1, Frac());
inv.assign(P / B1 + 1, 0);
pre.assign(B2 + 1, 0); suf.assign(B2 + 1, 0);
inv[1] = 1;
for (int i = 2; i <= P / B1; ++i)
{
inv[i] = 1ll * (P - P / i) * inv[P % i] % P;
}
s[0] = Frac(0, 1); s[B2] = Frac(1, 1);
pre[B2] = suf[B2] = B2;
for (int i = 2; i <= B1; ++i)
{
for (int j = 1; j < i; ++j)
{
int pos = 1ll * j * B2 / i;
if (pre[pos]) continue;
pre[pos] = suf[pos] = pos;
s[pos] = Frac(j, i);
}
}
for (int i = 1; i <= B2; ++i)if (!pre[i])pre[i] = pre[i - 1];
for (int i = B2; i >= 0; --i)if (!suf[i])suf[i] = suf[i + 1];
}
int operator[](int x)
{
if (x <= P / B1) return inv[x];
int pos = 1ll * x * B2 / P, res = Calc(x, s[pre[pos]]);
if (res == -1) res = Calc(x, s[suf[pos]]);
return res;
}
};
線性遞推法
class C_INV
{
private:
vector<int>inv;
public:
void Init(int n, int P)
{
inv.assign(n + 1, 0); inv[1] = 1;
for (int i = 2; i <= n; ++i)
{
inv[i] = 1ll * (P - P / i) * inv[P % i] % P;
}
}
int operator[](int x)
{
return inv[x];
}
};
擴充套件歐幾里得法
class C_INV
{
private:
int P;
/*====================*/
void exgcd(int a, int b, int& x, int& y)
{
if (b == 0)
{
x = 1, y = 0;
}
else
{
exgcd(b, a % b, y, x);
y -= a / b * x;
}
}
public:
void Init(int P)
{
this->P = P;
}
int operator[](int x)
{
int a, b;
exgcd(x, P, a, b);
return (a % P + P) % P;
}
};
質數
判定
六倍試除法
bool IsPrime(int x)
{
if (x <= 1)return false;
if (x == 2 || x == 3 || x == 5)return true;
if (x % 2 == 0 || x % 3 == 0)return false;
for (int i = 5; i * i <= x; i += 6)
{
if (x % i == 0 || x % (i + 2) == 0)
{
return false;
}
}
return true;
}
Miller-Rabin
class C_Prime
{
private:
lnt Mul(lnt a, lnt b, lnt p)
{
return (__int128)a * b % p;
}
lnt Pow(lnt a, lnt b, lnt p)
{
lnt res = 1;
while (b)
{
if (b % 2 == 1)
{
res = Mul(res, a, p);
}
b /= 2, a = Mul(a, a, p);
}
return res;
}
public:
bool operator()(lnt x)
{
if (x < 3 || x % 2 == 0) return x == 2;
lnt a = x - 1, b = 0;
while (a % 2 == 0) a /= 2, ++b;
//lnt lib[] = { 2,7,61 };
lnt lib[] = { 2,325,9375,28178,450775,9780504,1795265022 };
for (lnt r : lib)
{
lnt v = Pow(r, a, x);
if (v == 1 || v == x - 1 || v == 0) continue;
for (lnt j = 1; j <= b; j++)
{
v = Mul(v, v, x);
if (v == x - 1 && j != b) { v = 1; break; }
if (v == 1) return false;
}
if (v != 1) return false;
}
return true;
}
};
篩法
尤拉篩
template<int MAX_RANGE>
class C_PrimeSieve
{
private:
vector<int>prime;
bitset<MAX_RANGE + 1>notprime;
public:
C_PrimeSieve(void)
{
notprime[0] = notprime[1] = true;
for (int i = 2; i <= MAX_RANGE; ++i)
{
if (!notprime[i])prime.push_back(i);
for (auto p : prime)
{
if (i * p > MAX_RANGE)break;
notprime[i * p] = true;
if (i % p == 0)break;
}
}
}
int Size(void)
{
return prime.size();
}
bool operator()(int x)
{
return !notprime[x];
}
int operator[](int x)
{
return prime[x - 1];
}
};
埃拉託斯特尼篩法
template<int MAX_RANGE>
class C_PrimeSieve
{
private:
vector<int>prime;
bitset<MAX_RANGE + 1>notprime;
public:
C_PrimeSieve(void)
{
notprime[0] = notprime[1] = true;
for (int p = 2; p <= MAX_RANGE; ++p)
{
if (!notprime[p])
{
prime.push_back(p);
for (lnt k = p; k * p <= MAX_RANGE; ++k)
{
notprime[k * p] = true;
}
}
}
}
int Size(void)
{
return prime.size();
}
bool operator()(int x)
{
return !notprime[x];
}
int operator[](int x)
{
return prime[x - 1];
}
};
Wheel Factorization
namespace SIEVE {
using size_type = uint32_t;
using mask_type = unsigned char;
struct SieveNode {
size_type m_val, m_pos[8];
};
static constexpr size_type remainder_30[8] = { 1, 7, 11, 13, 17, 19, 23, 29 };
static constexpr size_type block_size = 7 * 11 * 13 * 17;
constexpr size_type get_estimated_ln(size_type x) {
return x <= 7 ? 1
: x <= 32 ? 2
: x <= 119 ? 3
: x <= 359 ? 4
: x <= 1133 ? 5
: x <= 3093 ? 6
: x <= 8471 ? 7
: x <= 24299 ? 8
: x <= 64719 ? 9
: x <= 175196 ? 10
: x <= 481451 ? 11
: x <= 1304718 ? 12
: x <= 3524653 ? 13
: x <= 9560099 ? 14
: x <= 25874783 ? 15
: x <= 70119984 ? 16
: x <= 189969353 ? 17
: x <= 514278262 ? 18
: x <= 1394199299 ? 19
: 20;
}
constexpr size_type get_estimated_Pi(size_type x) { return x / get_estimated_ln(x); }
constexpr size_type get_estimated_sqrt_Pi(size_type x) { return 5 << size_type(get_estimated_ln(x) * 0.55); }
template <size_type MAX_RANGE>
struct FastSieve {
static constexpr size_type max_r = (get_estimated_sqrt_Pi(MAX_RANGE) + block_size * 2 - 1) / (block_size * 2) * block_size * 2, max_pi = get_estimated_Pi(MAX_RANGE);
static size_type s_primes[max_pi], s_prime_cnt;
static mask_type s_masks[block_size], s_buffer[block_size];
static SieveNode s_nodes[max_r];
template <typename Callback>
static void _init_sieve(size_type range, Callback&& call) {
if (range < 19) return;
std::vector<bool> vis(range + 1);
for (size_type i = 3; i * i <= range; i += 2)
if (!vis[i])
for (size_type j = i * i; j <= range; j += i << 1) vis[j] = true;
for (size_type i = 19; i <= range; i += 2)
if (!vis[i]) call(i);
}
static void _add_prime(size_type p) { s_primes[s_prime_cnt++] = p; }
FastSieve(size_type range = MAX_RANGE) {
for (size_type p : {2, 3, 5, 7, 11, 13, 17})
if (p <= range) _add_prime(p);
if (range <= block_size) {
_init_sieve(range, _add_prime);
return;
}
std::fill_n(s_masks, block_size, -1);
for (size_type p : {7, 11, 13, 17})
for (size_type i = 0; i != 8; i++) {
size_type j = p;
while (j % 30 != remainder_30[i]) j += p << 1;
for (j /= 30; j < block_size; j += p) s_masks[j] &= ~(1 << i);
}
SieveNode* end = s_nodes;
auto add_node = [&](size_type p) {
for (size_type i = 0; i != 8; i++) {
size_type j = p * p;
while (j % 30 != remainder_30[i]) j += p << 1;
end->m_pos[i] = j / 30;
}
end++->m_val = p;
};
_init_sieve(sqrt(range), add_node);
for (size_type first = 0, tot = (range - 1) / 30 + 1; first < tot; first += block_size) {
size_type last = std::min(tot, first + block_size);
std::copy_n(s_masks, block_size, s_buffer);
if (!first) s_buffer[0] &= 0xfe;
for (auto it = s_nodes; it != end; ++it)
for (size_type p = it->m_val, i = 0; i != 8; i++) {
size_type j = it->m_pos[i];
for (; j < last; j += p) s_buffer[j - first] &= ~(1 << i);
it->m_pos[i] = j;
}
for (size_type i = first; i != last; i++)
for (mask_type mask = s_buffer[i - first]; mask;) {
size_type x = std::countr_zero(mask);
_add_prime(i * 30 + remainder_30[x]), mask -= size_type(1) << x;
}
}
while (s_primes[s_prime_cnt - 1] > range) s_prime_cnt--;
}
std::bitset<MAX_RANGE + 1> to_bitset() const {
std::bitset<MAX_RANGE + 1> res;
for (size_type i = 0; i != s_prime_cnt; i++) res.set(s_primes[i]);
return res;
}
size_type query_kth_prime(size_type k) const { return s_primes[k]; }
size_type count() const { return s_prime_cnt; }
};
template <size_type MAX_RANGE>
size_type FastSieve<MAX_RANGE>::s_primes[FastSieve<MAX_RANGE>::max_pi];
template <size_type MAX_RANGE>
size_type FastSieve<MAX_RANGE>::s_prime_cnt;
template <size_type MAX_RANGE>
mask_type FastSieve<MAX_RANGE>::s_masks[block_size];
template <size_type MAX_RANGE>
mask_type FastSieve<MAX_RANGE>::s_buffer[block_size];
template <size_type MAX_RANGE>
SieveNode FastSieve<MAX_RANGE>::s_nodes[FastSieve<MAX_RANGE>::max_r];
}
子集
列舉全部真子集
vector<int> AllSubset(int s)
{
vector<int>res;
for (int i = s; i; i = (i - 1) & s)
{
res.push_back(i);
}
return res;
}
列舉k大小真子集
vector<int> GospersHack(int k, int n)
{
vector<int>res;
int cur = (1 << k) - 1, limit = (1 << n);
while (cur < limit)
{
res.push_back(cur);
int lb = cur & -cur, r = cur + lb;
cur = (((r ^ cur) >> 2) / lb) | r;
}
return res;
}
快速冪
lnt Pow(lnt a, lnt b, const lnt P)
{
lnt res = 1;
while (b)
{
if (b & 1)
{
res = (res * a) % P;
}
b >>= 1, a = (a * a) % P;
}
return res;
}
龜速乘
int Mul(int a, int b, int p)
{
int res = 0;
while (b)
{
if (b & 1)
{
res = (res + a) % p;
}
b >>= 1, a = (a * 2) % p;
}
return res;
}
lnt Mul(lnt a, lnt b, lnt p)
{
return (a * b - (lnt)(a / (long double)p * b + 1e-3) * p + p) % p;
}
lnt Mul(lnt a, lnt b, lnt p)
{
return (__int128)a * b % p;
}
逆序對
template<class Type>
class _CountInversions
{
public:
lnt operator()(int n, Type arr[])
{
res = 0;
a = new Type[n + 1];
temp = new Type[n + 1];
for (int i = 1; i <= n; ++i)
{
a[i] = arr[i];
}
MergeSort(1, n);
delete[] a; delete[] temp;
return res;
}
private:
Type* a = NULL;
lnt res = 0;
Type* temp = NULL;
/*====================*/
void MergeSort(int l, int r)
{
if (l < r)
{
int i = l;
int mid = (l + r) >> 1;
int p = l, q = mid + 1;
MergeSort(l, mid + 0);
MergeSort(mid + 1, r);
while (p <= mid && q <= r)
{
if (a[p] <= a[q])
{
temp[i++] = a[p++];
}
else
{
temp[i++] = a[q++];
res += (lnt)(mid - p + 1);
}
}
while (p <= mid)temp[i++] = a[p++];
while (q <= r) temp[i++] = a[q++];
for (i = l; i <= r; i++)a[i] = temp[i];
}
}
};
多項式
快速傅立葉變換
typedef complex<double> Comp;
const ll SIZE = 4000000 + 10;
const double PI = acos(-1.0);
Comp temp[SIZE];
void FFT(Comp* p, ll len, ll inv)//inv=+1 FFT;inv=-1 IFFT
{
if (len == 1) return;
const int E = 0, O = len / 2;
for (ll i = 0; i < len; ++i) temp[i] = p[i];
for (ll i = 0; i < len; ++i)
{
if ((i & 1) == 1)p[i / 2 + O] = temp[i];
if ((i & 1) == 0)p[i / 2 + E] = temp[i];
}
Comp* pe = p + E; FFT(pe, len / 2, inv);
Comp* po = p + O; FFT(po, len / 2, inv);
Comp omega(1, 0);
const double Angle = 2 * PI / len;
const Comp step(cos(Angle), sin(inv * Angle));
for (ll k = 0; k < len / 2; ++k, omega *= step)
{
temp[k + E] = pe[k] + omega * po[k];
temp[k + O] = pe[k] - omega * po[k];
}
for (ll i = 0; i < len; ++i) p[i] = temp[i];
}
離散對數
ll BSGS(ll a, ll b, ll m)
{
static unordered_map<ll, ll> hs;
hs.clear();
ll cur = 1, t = sqrt(m) + 1;
for (int B = 1; B <= t; ++B)
{
(cur *= a) %= m;
hs[b * cur % m] = B; // 雜湊表中存B的值
}
ll now = cur; // 此時cur = a^t
for (int A = 1; A <= t; ++A)
{
auto it = hs.find(now);
if (it != hs.end())
return A * t - it->second;
(now *= cur) %= m;
}
return -1; // 沒有找到，無解
}
基數排序
template <int B = 8, class Type>
void RadixSort(vector<Type>& a)
{
const int mask = (1 << B) - 1, n = a.size();
vector<Type> b(n); vector<int> cnt(1 << B);
Type maxV = *max_element(a.begin(), a.end());
for (int i = 0; maxV; i += B, maxV >>= B)
{
fill(cnt.begin(), cnt.end(), 0);
for (int j = 0; j < n; j++)cnt[a[j] >> i & mask] += 1;
for (int j = 1; j < (1 << B); j++)cnt[j] += cnt[j - 1];
for (int j = n - 1; j >= 0; j--)b[--cnt[a[j] >> i & mask]] = a[j];
swap(a, b);
}
}
尤拉函式
試除法
class PHI
{
public:
int operator[](int x)
{
return GetPhi(x);
}
private:
int GetPhi(int x)
{
int res = x;
for (int i = 2; i * i <= x; ++i)
{
if (x % i == 0)
{
res = res / i * (i - 1);
while (x % i == 0) x /= i;
}
}
if (x > 1) res = res / x * (x - 1);
return res;
}
};
尤拉篩法
class PHI
{
public:
~PHI(void)
{
delete[] phi;
}
void init(int n)
{
GetPhi(n);
}
int operator[](int x)
{
return phi[x];
}
private:
int* phi = NULL;
void GetPhi(int n)
{
phi = new int[n + 1];
bool* vis = new bool[n + 1];
int* table = new int[n + 1];
for (int i = 0; i <= n; ++i)
{
vis[i] = true;
}
int cnt = 0; phi[1] = 1;
for (int i = 2; i <= n; ++i)
{
if (vis[i])
{
phi[i] = i - 1;
table[++cnt] = i;
}
for (int j = 1; j <= cnt; ++j)
{
if (i * table[j] > n)break;
vis[i * table[j]] = false;
if (i % table[j] == 0)
{
phi[i * table[j]] = phi[i] * table[j]; break;
}
else
{
phi[i * table[j]] = phi[i] * (table[j] - 1);
}
}
}
delete[] vis; delete[] table;
}
};
尤拉降冪
class EX_Euler
{
public:
int operator()(int a, string s, int p)
{
int b = 0;
bool flag = false;
int phi = GetPhi(p);
for (auto c : s)
{
b = (b * 10 + c - '0');
if (b >= phi)flag = true, b %= phi;
}
if (flag)b += phi; return Pow(a % p, b, p);
}
private:
int GetPhi(int x)
{
int res = x;
for (int i = 2; i * i <= x; ++i)
{
if (x % i == 0)
{
res = res / i * (i - 1);
while (x % i == 0) x /= i;
}
}
if (x > 1) res = res / x * (x - 1);
return res;
}
int Pow(int a, int b, int p)
{
int res = 1;
while (b != 0)
{
if (b % 2 == 1)
{
res = (res * a) % p;
}
b /= 2, a = (a * a) % p;
}
return res;
}
};
歐幾里得
最大公因數
int gcd(int a, int b)
{
return b == 0 ? a : gcd(b, a % b);
}
最小公倍數
int lcm(int a, int b)
{
return a / gcd(a, b) * b;
}
擴充套件歐幾里得
void exgcd(int a, int b, int& x, int& y)
{
if (b == 0)
{
x = 1, y = 0;
}
else
{
exgcd(b, a % b, y, x);
y -= a / b * x;
}
}
分解質因數
Pollard-Rho
class C_PrimeFactorization
{
private:
lnt Mul(lnt a, lnt b, lnt p)
{
return (__int128)a * b % p;
}
lnt Pow(lnt a, lnt b, lnt p)
{
lnt res = 1;
while (b)
{
if (b % 2 == 1)
{
res = Mul(res, a, p);
}
b /= 2, a = Mul(a, a, p);
}
return res;
}
bool MillerRabin(lnt x)
{
if (x < 3 || x % 2 == 0) return x == 2;
lnt a = x - 1, b = 0;
while (a % 2 == 0) a /= 2, ++b;
//lnt lib[] = { 2,7,61 };
lnt lib[] = { 2,325,9375,28178,450775,9780504,1795265022 };
for (lnt r : lib)
{
lnt v = Pow(r, a, x);
if (v == 1 || v == x - 1 || v == 0) continue;
for (lnt j = 1; j <= b; j++)
{
v = Mul(v, v, x);
if (v == x - 1 && j != b) { v = 1; break; }
if (v == 1) return false;
}
if (v != 1) return false;
}
return true;
}
/*====================*/
lnt GCD(lnt a, lnt b)
{
return b == 0 ? a : GCD(b, a % b);
}
lnt GetFactor(lnt X)
{
if (X == 4)return 2;
if (MillerRabin(X))return X;
mt19937 rng(time(NULL));
while (1)
{
lnt c = rng() % (X - 1) + 1;
auto f = [=](lnt x) { return ((__int128)x * x + c) % X; };
lnt t = 0, r = 0, p = 1, q;
do
{
for (int i = 0; i < 128; ++i)
{
t = f(t), r = f(f(r));
if (t == r || (q = (__int128)p * abs(t - r) % X) == 0)break;
p = q;
}
lnt d = GCD(p, X); if (d > 1)return d;
} while (t != r);
}
}
void GetAllFactor(lnt X, vector<lnt>& lib)
{
lnt fac = GetFactor(X);
if (fac == X)lib.push_back(fac);
else GetAllFactor(fac, lib), GetAllFactor(X / fac, lib);
}
public:
void operator()(lnt x, vector<lnt>& num, vector<lnt>& cnt)
{
num.clear(), cnt.clear();
GetAllFactor(x, num);
sort(num.begin(), num.end());
for (int i = 0; i < num.size(); ++i)
{
if (i == 0 || num[i] != num[i - 1])
{
cnt.push_back(0);
}
cnt.back()++;
}
num.erase(unique(num.begin(), num.end()), num.end());
}
};
預處理質數表法
class C_PrimeFactorization
{
private:
vector<int>prime;
public:
void Init(int n)
{
vector<bool>isprime;
isprime.assign(n + 1, true);
isprime[0] = isprime[1] = false;
for (int i = 2; i <= n; ++i)
{
if (isprime[i])prime.push_back(i);
for (int j = 0; j < prime.size(); ++j)
{
if (i * prime[j] > n)break;
isprime[i * prime[j]] = false;
if (i % prime[j] == 0)break;
}
}
}
void operator()(int x, vector<int>& num, vector<int>& cnt)
{
num.clear(), cnt.clear();
for (auto p : prime)
{
if (p * p > x)break;
if (x % p == 0)
{
num.push_back(p), cnt.push_back(0);
while (x % p == 0)cnt.back()++, x /= p;
}
}
if (x != 1)
{
num.push_back(x), cnt.push_back(1);
}
}
};
預處理最小質因子法
class C_PrimeFactorization
{
private:
vector<int>mpf;
public:
void Init(int n)
{
vector<int>prime;
vector<bool>isprime;
mpf.assign(n + 1, 0);
isprime.assign(n + 1, true);
isprime[0] = isprime[1] = false;
for (int i = 2; i <= n; ++i)
{
if (isprime[i])prime.push_back(i), mpf[i] = i;
for (int j = 0; j < prime.size(); ++j)
{
if (i * prime[j] > n)break;
isprime[i * prime[j]] = false;
mpf[i * prime[j]] = prime[j];
if (i % prime[j] == 0)break;
}
}
}
void operator()(int x, vector<int>& num, vector<int>& cnt)
{
num.clear(), cnt.clear();
while (x > 1)
{
if (num.empty() || num.back() != mpf[x])
{
num.push_back(mpf[x]); cnt.push_back(1);
}
else
{
cnt.back()++;
}
x /= mpf[x];
}
}
};
獲得全部因數
void GetDivisor(int x, vector<int>& divisor)
{
vector<int>num, cnt;
PFF(x, num, cnt);
divisor.push_back(1);
for (int i = 0; i < num.size(); ++i)
{
int val = 1;
int lim = divisor.size();
for (int j = 1; j <= cnt[i]; ++j)
{
val *= num[i];
for (int k = 0; k < lim; ++k)
{
divisor.push_back(divisor[k] * val);
}
}
}
}
隨機數生成器
mt19937 Rand(random_device{}());