筛质数
埃式筛,复杂度
class Solution {
public:
int countPrimes(int n) {
vector<int> isPrime(n, 1);
int ans = 0;
for (int i = 2; i < n; ++i) {
if (isPrime[i]) {
ans += 1;
if ((long long)i * i < n) {
for (int j = i * i; j < n; j += i) {
isPrime[j] = 0;
}
}
}
}
return ans;
}
};
线性筛,复杂度
class Solution {
public:
int countPrimes(int n) {
vector<int> primes;
vector<int> isPrime(n, 1);
for (int i = 2; i < n; ++i) {
if (isPrime[i]) {
primes.push_back(i);
}
for (int j = 0; j < primes.size() && i * primes[j] < n; ++j) {
isPrime[i * primes[j]] = 0;
if (i % primes[j] == 0) {
break;
}
}
}
return primes.size();
}
};
预处理质因数个数
int f[N + 10];
bool inited = false;
// 预处理每个数有几个质因数
void init() {
if (inited) return;
inited = true;
memset(f, 0, sizeof(f));
for (int i = 2; i <= N; i++)
if (!f[i])
for (int j = i; j <= N; j += i) f[j]++;
}
预处理每个数的质因数
bool inited = false;
vector<int> fac[N + 1];
void init() {
if (inited) return;
inited = true;
for (int i = 2; i <= N; i++)
if (fac[i].empty())
for (int j = i; j <= N; j += i) fac[j].push_back(i);
}
void divide(int x) {
for (int i = 2; i <= x / i; i++) // i <= x / i:防止越界,速度大于 i < sqrt(x)
if (x % i == 0) { // i为底数
int s = 0; // s为指数
while (x % i == 0) x /= i, s++;
cout << i << ' ' << s << endl; //输出
}
if (x > 1) cout << x << ' ' << 1 << endl; //如果x还有剩余,单独处理
cout << endl;
}
结论:预处理的复杂度是 的
平均来看,每个自然数有 个真因子
判断公因子
枚举 和 的 GCD 的因子
组合数
对(3)的推导
之前的模板
int c[N][N];
void init() { //组合数公式
for (int i = 0; i < N; ++i) {
for (int j = 0; j <= i; ++j)
if (j == 0)
c[i][j] = 1;
else
c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
}
}
//用求逆元的方式求出分母阶乘的逆元 预处理阶乘
int fact[N], infact[N];
LL qmi(int a, int k, int p) {
LL res = 1;
while (k) {
if (k & 1) res = res * a % p;
a = (LL)a * a % p;
k >>= 1;
}
return res;
}
void init_1() {
fact[0] = infact[0] = 1;
for (int i = 1; i < N; ++i) {
fact[i] = (LL)fact[i - 1] % mod;
infact[i] = (LL)infact[i - 1] * qmi(i, mod - 2, mod);
}
}
//数据爆大 Lucus定理
// C(a, b) 模p同余 C(a mod p, b mod p) * C(a / p, b / p)
int p;
int C(int a, int b) {
int res = 1;
for (int i = 1, j = a; i <= b; ++i, --j) {
res = (LL)res * j % p;
res = (LL)res * qmi(i, p - 2, p) % p;
}
return res;
}
int lucus(LL a, LL b) {
if (a < p && b < p) return C(a, b);
return (LL)C(a % p, b % p) * lucus(a / p, b / p) % p;
}
新模板,包含逆元
LL qmi(LL x, int n)
{
LL res = 1;
for (; n; n /= 2)
{
if (n % 2) res = res * x % MOD;
x = x * x % MOD;
}
return res;
}
LL fac[N], inv[N];
auto init = [] {
fac[0] = 1;
for (int i = 1; i < N; i++)
fac[i] = fac[i - 1] * i % MOD;
inv[N - 1] = qmi(fac[N - 1], MOD - 2);
for (int i = N - 1; i; i--)
inv[i - 1] = inv[i] * i % MOD;
return 0;
}();
LL comb(int n, int k)
{
return (k < 0 || n < k) ? 0 : fac[n] * inv[k] % MOD * inv[n - k] % MOD;
}
关于逆元,当 为质数时,可以用快速幂来求,否则就只能用扩展欧几里得法来求
如何求解 个数的逆元?
首先求前缀积,然后求逆元,因为 sv[n]
是 个数积的逆元,把它乘上 a[n]
,就会和 a[n]
的逆元抵消,得到 sv[n-1]
,以此类推
s[0] = 1;
for (int i = 1;i <= n; i++) s[i] = s[i - 1] * a[i] % p;
sv[n] = qmi(s[n], p - 2);
for (int i = n; i >= 1; i--) sv[i - 1] = sv[i] * a[i] % p;
for (int i = 1; i <= n; i++) inv[i] = sv[i] * s[i - 1] % p;
不用快速幂来求逆元的模板
LL fac[N + 10], ifac[N + 10];
auto init = [] {
fac[0] = 1;
for (int i = 1; i <= N; i++) {
fac[i] = fac[i - 1] * i % MOD;
}
ifac[0] = ifac[1] = 1;
for (int i = 2; i <= N; i++) {
ifac[i] = (MOD - MOD / i) * ifac[MOD % i] % MOD;
}
for (int i = 2; i <= N; i++) {
ifac[i] = ifac[i - 1] * ifac[i] % MOD;
}
return 0;
}();
LL C(int a, int b) {
if (a < b) {
return 0;
}
return fac[a] * ifac[b] % MOD * ifac[a - b] % MOD;
}
组合数前缀和的递推公式
由 C[i][j]=C[i-1][j]+C[i-1][j-1]
可知,上一行的所有元素都会转移给下一行的两个元素,因此设 ,有
例题可见 LC3428
滚动数组求 C(n,*)
如果只用到第 层的 ,可以直接用数组递推
c[0] = 1;
for (int i = 1; i <= n; i++) {
for (int j = i; j > 0; j--) {
c[j] = (c[j] + c[j - 1]) % mod;
}
}
模数不是质数的情况
LC3463 中出现了需要求 的情况,不能直接套用逆元的模板来求
方法一:求出每个数中因子 和因子 的个数,单独计算,即 可以写成 ,其中 、 和 与 互质,可以通过欧拉定理计算它们在模 下的逆元;对于 和 的次幂直接计算即可
int comb(int n, int k) {
// 由于每项都 < 10,所以无需中途取模
return f[n] * inv_f[k] * inv_f[n - k] *
qpow(2, p2[n] - p2[k] - p2[n - k]) *
qpow(5, p5[n] - p5[k] - p5[n - k]) % MOD;
}
预处理的代码为
auto init = []() {
f[0] = 1;
for (int i = 1; i <= MX; i++) {
int x = i;
// 计算 2 的幂次
int e2 = countr_zero((unsigned) x);
x >>= e2;
// 计算 5 的幂次
int e5 = 0;
while (x % 5 == 0) {
e5++;
x /= 5;
}
f[i] = f[i - 1] * x % MOD;
p2[i] = p2[i - 1] + e2;
p5[i] = p5[i - 1] + e5;
}
// 欧拉定理 a^(phi(m)) = 1 (mod m)
// 这里 m = 10
// phi(m) 为小于等于 10 的与 10 互质的数的个数
// 逆元就是 a^(phi(m)-1) = a^3
inv_f[MX] = qpow(f[MX], 3);
for (int i = MX; i > 0; i--) {
int x = i;
x >>= countr_zero((unsigned) x);
while (x % 5 == 0) {
x /= 5;
}
inv_f[i - 1] = inv_f[i] * x % MOD;
}
return 0;
}();
方法二:利用 Lucas 定理和中国剩余定理(CRT)
先通过 Lucas 定理计算出 和 的结果
根据 CRT 的计算过程,现在有两条方程: 和 ,这里的 就是 。我们要求解这个 ,计算出 ,,那么 ,就可以求出答案
// 计算 C(n, k) % p,要求 p 是质数
int lucas(int n, int k, int p) {
if (k == 0) {
return 1;
}
return c[n % p][k % p] * lucas(n / p, k / p, p) % p;
};
int comb(int n, int k) {
// 结果至多为 5 + 4 * 6 = 29,无需中途取模
return lucas(n, k, 2) * 5 + lucas(n, k, 5) * 6;
}
循环计算简单组合数
参考题解,通过循环枚举分子的每一个乘数,然后把当前乘积除以 计算部分结果。由于连续 个数中必有 的倍数,计算过程一定是整数
long long comb(int n, int k) {
k = min(k, n - k);
long long res = 1;
for (int i = 1; i <= k; i++) {
res = res * (n + 1 - i) / i;
}
return res;
}
全排列
康托展开用于计算给定全排列的 rank,逆康托展开用于给定 rank 反推全排列
基于这样一条公式:
其中
原理是计算有多少个排列的字典序小于给定排列:固定前 个数相同,第 个数要小于 ,只能从后面找。这样枚举是不重不漏的
// 来自 oi-wiki
class BIT {
int n;
std::vector<int> su;
public:
BIT(int n) : n(n), su(n + 1) {}
// Add v to the x-th number.
void add(int x, int v) {
for (; x <= n; x += x & (-x)) {
su[x] += v;
}
}
// Get the cumulative sum till the x-th number.
int query(int x) {
int res = 0;
for (; x; x &= x - 1) {
res += su[x];
}
return res;
}
};
// Get the rank of a permutation of 1~n.
long long find_rank(const std::vector<int>& nums) {
int n = nums.size();
BIT bit(n);
long long fac = 1;
long long res = 0;
// Reverse iteration.
for (int i = n - 1; i >= 0; --i) {
// Count the number of elements smaller than the current one.
res += bit.query(nums[i] - 1) * fac;
// Insert the current element into the BIT.
bit.add(nums[i], 1);
// Update the factorial.
fac *= n - i;
}
return res + 1;
}
那怎么反求排列呢?需要利用 的公式反求出 数组。这里的计算原理是倒着求,每次对 取模,前面的都消掉,剩下的就是 ,然后需要找到当前未使用的数中第 小的数,用到了一个树状数组的技巧
// 来自 oi-wiki
class BIT {
int n;
std::vector<int> su;
public:
BIT(int n) : n(n), su(n + 1) {}
// Fill the BIT with one.
// 等同于求每个点的管辖范围有多少个数
void fill() {
for (int x = 1; x <= n; ++x) {
su[x] += x & (-x);
}
}
// Add v to the x-th number.
void add(int x, int v) {
for (; x <= n; x += x & (-x)) {
su[x] += v;
}
}
// Get the k-th smallest element.
// 这个技巧很关键
int find_kth(int k) {
int ps = 0, x = 0;
for (int i = log2(n); i >= 0; --i) {
x += 1 << i;
if (x >= n || ps + su[x] >= k) {
x -= 1 << i;
} else {
ps += su[x];
}
}
return x + 1;
}
};
// Find the k-th permutation of 1~n.
std::vector<int> find_permutation(int n, long long k) {
--k;
std::vector<int> A(n);
for (int i = 1; i <= n; ++i) {
// 这里的计算过程要结合样例理解
A[n - i] = k % i;
k /= i;
}
BIT bit(n);
// Set all values in BIT to one.
bit.fill();
std::vector<int> res(n);
for (int i = 0; i < n; ++i) {
// Find the A[i]-th smallest unused element.
res[i] = bit.find_kth(A[i] + 1);
// Remove it from the BIT.
bit.add(res[i], -1);
}
return res;
}
计算限制下的第 k 个全排列
题源 LC3470,要求相邻位置的奇偶性不同
枚举每一位填什么,难点是可能会爆数据范围,技巧是当剩下的排列数大于 时,直接找当前未用的最小数即可
using LL = long long;
vector<LL> f{1};
auto init = []() {
for (int i = 1; f.back() < 1e15; i++) {
f.push_back(f.back() * i);
}
return 0;
}();
class Solution {
public:
vector<int> permute(int n, long long k) {
vector<bool> st(n + 1);
// now:现在该填偶数还是奇数,-1 表示还不确定(n 为偶数时,第一位可以填任意数)
int now = (n & 1 ? 1 : -1);
vector<int> res;
for (int i = 1; i <= n; i++) {
int r = n - i;
// 计算剩下 n - i 个数任意组合,能产生几种方案
LL t;
int x;
if (r / 2 >= f.size() || f[r / 2] > k / f[r / 2]) {
x = 1;
} else {
if (r & 1) {
t = f[(r + 1) / 2] * f[r / 2];
} else {
t = f[r / 2] * f[r / 2];
}
// 我们需要填第 x 小的数
x = (k + t - 1) / t;
// 第 i 位填第 x 小的数,也就是说我们跳过了 t(x - 1) 种方案
// 那么要算出剩下的 n - i 位中,第 k - t(x - 1) 大的方案
k -= t * (x - 1);
}
// 通过枚举找出第 x 小的能填的数
for (int j = 1; j <= n; j++) {
if (!st[j] && (now == -1 || j % 2 == now)) {
x --;
if (x == 0) {
st[j] = true;
now = 1 - (j % 2);
res.push_back(j);
break;
}
}
}
if (x > 0) {
return {};
}
}
return res;
}
};