Chengyuan Ma's Blog

幽居默处而观万物之变，尽其自然之理，而断之于中

蒟蒻的数论笔记——素数

Posted on 2018-09-01

定义：只能被和自身整除的数。

素数的判定

试除法

如果不是素数那么它必定有一个不大于 $\sqrt n$ 的因子，因此我们枚举 $[1, \sqrt n]$ 中的所有数挨个试除即可，时间复杂度 $O(\sqrt n)$ ，虽然慢但是保证正确。

费马判定法

逆用费马小定理，即若为素数则 $a^{n - 1} \equiv 1 \pmod n$ ，对于一个可能的我们尝试不同的并判断是否符合定理，如果不是那么一定不是素数，如果是，那么对于单个的我们正确的概率约为 $1 \over 4$ ，因此我们可以多试几次（比如次），这样基本上就可以保证正确性了。时间复杂度就是快速幂的复杂度 $O(\log n)$ 。

缺点：会被卡迈克尔数干死。如果合数是卡迈克尔数，那么对于所有和互素的，也有 $a^{n - 1} \equiv 1 \pmod n$ ，简而言之就是费马小定理逆定理的反例。最小的卡迈克尔数是 $561 = 3 \times 11 \times 17$ 。虽然这种数很稀少但是，在 $[1, 1 \times 10^8]$ 之间总共有 255 个卡迈克尔数，出题人想干死你还是很简单。

Miller-Rabin 判定法

引理：模素数意义下的平方根模为 $\pm 1$ 。

证明：设此时的平方根为：

$\begin{aligned} x^2 &\equiv 1 &\pmod n\\ x^2 - 1 &\equiv 0 &\pmod n\\ (x + 1)(x - 1) &\equiv 0 &\pmod n\\ \end{aligned}$

因此 $n\mid (x + 1)$ 或 $n\mid (x - 1)$ ，即 $x \equiv \pm1 \pmod n$ ，得证。

对于一个奇素数，偶数 n - 1 一定可以写成 2^sd 的形式，其中是奇数。我们断言对于所有和某个以下两个不等式必有一个成立：

$\begin{aligned} a^d &\equiv 1 &\pmod n \\ a^{2^rd} &\equiv -1 &\pmod n \end{aligned}$

其中 $0 < a < n, 0 \le r < s$ 。

证明：由费马小定理可得 $a^{2^sd} \equiv 1 \pmod n$ ，我们不断地对这个式子两边开根，由引理，如果开出，那么下式成立，如果开了次没开出，则上式成立。

Miller-Rabin 判别法基于以上命题的逆否命题，即如果存在一个，使得 $a^d \not \equiv 1 \pmod n$ 且对于所有的 $0 \le r < s$ 都有 $a^{2^rd} \not\equiv -1 \pmod n$ ，那一定不是素数。反之，就很有可能是素数。

因此，Miller Rabin的工作流程如下：

bool millerRabin(int n) { // n > 3
    int s = 0, d = n - 1;
    while (!(d & 1)) s++, d >>= 1; // 分解s, d
    for (int i = 1; i <= K; i++) { // K 越大越准
        int a = rand() % (n - 3) + 2; // 随机一个[2, n - 2]的a
        int x = qpow(a, d, n); // a^d mod n, qpow快速幂
        if (x == 1 || x == n - 1) continue; // 不满足逆否命题, n - 1是r = 0时的下式
        bool flag = false;
        for (int j = 1; j <= s - 1; j++) {
            x = x * x % n;
            if (x == n - 1) {
                flag = true; // 不满足逆否命题
                break;
            }
        }
        if (flag) continue;
        return false;
    }
    return true;
}

要我说还是蛮复杂的，但是这个测试就已经不会被卡迈克尔数干死了。而且在 K = 1 时就已经很准了。

时间复杂度为 $O(n\log n)$ 。

自己的判定法？

可以通过研究素数的更多性质，例如高次剩余的特性等等发明自己的判别法，如果在复杂度允许的情况下能够在int范围内没问题也可以用。

计算的所有素数

如何快速获得 [1, n] 的所有素数呢？如果对于个数都进行以上的素数判定肯定是非常慢的，因此我们需要筛法，顾名思义，就是把不是素数的全部筛掉，那么剩下的就都是素数了。

埃拉托色尼筛法

最朴素的筛法，基于的思想是如果我把每一个素数的倍数都标记成合数的话那么剩下没有标记的就是素数了，代码非常简洁，如下：

bool v[N]; // 不是质数
int prime[P], tot = 0;
for (int i = 2; i <= n; i++) {
    if (v[i]) continue;
    prime[++tot] = i;
    for (int j = 2; j <= n / i; j++)
        v[i * j] = true;
}

其复杂度应该为 $\sum_{素数p \le n} \left\lfloor \frac{n}{p} \right\rfloor$ 。由麦腾斯定理，这个式子在 $n \to \infty$ 时逼近 $n\ln\ln n$ ，因此埃拉托色尼筛法的时间复杂度为 $O(n\ln\ln n)$ ，一般可以达到 10^7 。这个筛法有一个小优化我也在之前的笔记中有提到。

埃拉托色尼筛法的精髓在于它利用了约数和倍数的关系是成对的这一点，将素数判断这一个枚举约数的问题转化为了枚举倍数的问题并加快了速度，这一思想十分重要，在计算狄利克雷卷积和反演的时候也会用到（把对因数求和转化为对倍数算贡献）。

欧拉筛

虽然 $\ln\ln n$ 这个因子已经很小了，但是面对 n = 10^8 的时候还是会跪掉。因此我们需要一种筛法能够在 O(n) 的时间内筛出所有的素数。

我们发现埃氏筛效率欠佳的一点是：一个合数会被每一个它的素因数筛掉共计好几次。我们能不能让一个合数只被筛一次呢？

bool v[N]; // 不是质数
int prime[P], tot = 0;
for (int i = 2; i <= n; i++) {
    if (!v[i]) prime[++tot] = i;
    for (int j = 1; i * prime[j] <= n; j++) {
        v[i * prime[j]] = true;
        if (i % prime[j] == 0) break;
    }
}

欧拉筛的精髓在于那句if (i % prime[j] == 0) break;。

结论：在欧拉筛中每个合数会且仅会被它的最小素因子筛去。

证明：设合数的最小素因子为 p_o ，因为 $\frac{N}{p_0}$ 不可能被 p_0 更小的素数所整除，因此在循环到 $i = \frac{N}{p_0}$ 时，内循环一定会在break之前筛到 $\frac{N}{p_0} \cdot p_0 = N$ 。而假设有有另外一个素因子 p_1 > p_0 ，则显然 $p_0 \mid \frac{N}{p_1}$ ，因此在循环到 $i = \frac{N}{p_1}$ 时，内循环一定会在筛到（prime[j]为） p_0 之后就立刻break，因此一定不会再筛到 $\frac{N}{p_1} \cdot p_1$ 。

推论：因为合数的数量少于，因此欧拉筛的复杂度为 O(n) 。

欧拉筛原名线性筛，冠以欧拉之名纯粹是因为我们不仅用它筛素数，还可以顺便计算以欧拉函数为首的一票积性函数。

素因数分解

试除法

如果为合数那么一定有一个不大于 $\sqrt n$ 的素因子，因此我们枚举 $[1, \sqrt n]$ 中的所有数挨个试除即可，时间复杂度 $O(\sqrt n)$ ，是最常用的方法。

Pollard's-Rho方法

引理（生日悖论）：在 [1, N] 当中随机选出若干个数，平均需要抽出 $O(\sqrt N)$ 个数才会抽到两个相同的数。

简介：抽前 k - 1 个数两两不同的概率为：

$\frac{P(N, k - 1)}{N^{k - 1}} = \frac{N!}{(N - k+1)!N^{k-1}}$

如果你抽第

个数，与之前

个数相同的概率为 $\frac{k - 1}{N}$ ，因此你恰巧在刚刚抽完第

个数后得到两个相同的数的概率为：

$\frac{N!}{(N - k+1)!N^{k-1}} \cdot \frac{k - 1}{N} = \frac{(k - 1)N}{(N - k +1)!N^k}$

因此为了得到两个相同的数你需要抽出的数的期望是：

$\sum_{k = 2}^{N+1} \frac{k(k - 1)N}{(N - k +1)!N^k}$

它也等价于TAOCP中的：

$1 + \sum_{k = 1}^N \frac{N!}{(N-k)!N^k}$

神奇的数学家拉马努金曾对该函数进行过研究，并指出它渐进逼近 $\sqrt {\frac{\pi N}{2}}$ 。

好的，现在假设我们需要分解出 n = pq 的一个非平凡因子，假设我们有一个多项式同余函数（例如 $g(x) = (x^2 + 1) \bmod n$ ）。我们使用来迭代生成一个伪随机序列 $\{x_n\}$ ，它和其模的序列 $\{x_n \bmod p\}$ 相关联。由之前的生日悖论引理，前者平均在 $O(\sqrt n)$ 次迭代后产生循环，而后者只需要 $O(\sqrt p) < O(\sqrt n)$ 次。我们可以在迭代中检测这个循环，从而完成对于素因数的分解。

我们知道 $\{x_n \bmod p\}$ 产生循环当且仅当存在 x_i, x_j ，有 $x_i \equiv x_j \pmod p$ ，此时 x_i 和 x_j 之间就是这个序列的一个循环节，但是我们预先是不知道的，那怎么判断以上条件是否成立呢？我们发现以上条件等价于 $p \mid |x_i - x_j|$ ，因为同时 $p \mid n$ ，所以此时最大可以为 $\gcd(n, |x_i - x_j|)$ ，因此如果我们发现 $\gcd(n, |x_i - x_j|) \neq 1$ 我们就已经找到了一个 $p = \gcd(n, |x_i - x_j|)$ ，非常巧妙！

在算法实现中，我们使用弗洛伊德判圈算法，其精髓就是环上的追及问题，我们维护两个指针 x_i, x_j ，从出发点开始 x_i 一个一个地跳， x_j 两个两个地跳，他们的速度差为，只要出现了了一个环那么在某一时刻 x_j 总是可以套 x_i 一圈！每一轮过后我们检测$(n, |x_i - x_j|) $，如果不为那说明我们找到一个因子，但是如果为的话就是 $\{x_n\}$ 和 $\{x_n \bmod p\}$ 一起循环了，我们的算法就GG了，需要重新选定一个出发点（给的种子）。

这个启发式算法的复杂度取决于 $O(\sqrt p)$ ，但是因为一般来说 $p \le \sqrt n$ ，因此算法的复杂度平均下来是 $O(n^{1 \over 4})$ 的，实测在int范围内快得飞起。

整除分块

Posted on 2018-09-01

一维情况

在莫比乌斯反演中，我们经常需要计算类似如下的和式：

$\sum_{i = 1}^n \left\lfloor \frac{n}{i} \right\rfloor$

如果朴素计算的话时间复杂度是

的，这个代价对于很多多询问的题目来说是在是太高了，因此我们需要更快的方法。

结论： $\left\lfloor \frac{n}{i} \right\rfloor$ 总共只有 $O(\sqrt n)$ 种不同的取值。

证明：我们首先证明：对于所有的 $k \le \lfloor \sqrt n \rfloor$ 我们都可以找到一个使得 $\lfloor n / i \rfloor = k$ ：

令 $i = \left\lfloor n / k\right\rfloor \ge k$ ：

$ik \le n < (i + 1)k = ik + k \le ik + i = i(k + 1) \\ \Rightarrow k \le n / i < k + 1 \\ \Rightarrow \left\lfloor n / i \right\rfloor = k$

因此，通过适当地选择

我们已经可以让 $\lfloor n / i \rfloor$ 取到 $\lfloor \sqrt n\rfloor$ 个值，因为这些 $i \ge \lfloor \sqrt n\rfloor$ ，因此我们进一步推得：对于所有 $i \ge \lfloor \sqrt n\rfloor$ ， $\lfloor n / i \rfloor \le \lfloor \sqrt n\rfloor$ 共有 $\lfloor \sqrt n\rfloor$ 种取值。

对于所有的 $i \le \lfloor \sqrt n \rfloor$ ，我们证明对于这些， $\lfloor n / i \rfloor$ 的值各不相同，因为假设存在一对相邻的 i - 1, i ，使得 $\lfloor {n \over i - 1} \rfloor = \lfloor {n \over i} \rfloor = k$ ，则：

$k \le {n \over i - 1} < k + 1 \\ k \le {n \over i} < k + 1 \\ \Rightarrow {n \over i - 1} - {n \over i} = {n \over i(i - 1)} < 1$

而事实上因为 $i \le \lfloor \sqrt n\rfloor$ ， $i(i - 1) < n \Rightarrow {n \over i(i - 1)} > 1$ ，矛盾。

因此对于 $i \le \lfloor \sqrt n\rfloor$ 的 $\lfloor \sqrt n\rfloor$ 个， $\lfloor n / i \rfloor \ge \lfloor \sqrt n\rfloor$ 此时也有 $\lfloor \sqrt n\rfloor$ 种取值。

因此总共大约有 $2\lfloor \sqrt n \rfloor$ 种取值，证毕。

这样子我们就可以一段一段地计算了，时间复杂度骤降为 $O(\sqrt n)$ ，但是一段到底到哪里结束呢？

结论：值为的一段结束于 $i = \lfloor n / k\rfloor$ 。

证明：随着的变大， $n / i - \lfloor n / i\rfloor$ 的值必定单调递减，临界点就是 $n / i = \lfloor n / i\rfloor = k$ ，此时 i = n / k 。考虑到此时可能不是整数，因此结束位置为 $\lfloor n / k\rfloor$ 。

推论：所在的一段结束于 $\lfloor n / \lfloor n / i \rfloor \rfloor$ 。

有这两条结论，我们就可以很快地计算一开头和式的值了：

for (int l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    ans += (r - l + 1) * (n / l);
}

二维情况

假设 n < m ，我们有时候还会碰到类似如下的和式：

$\sum_{i = 1}^n \left\lfloor \frac{n}{i} \right\rfloor \left\lfloor \frac{m}{i} \right\rfloor$

结论： $\left\lfloor \frac{n}{i} \right\rfloor \left\lfloor \frac{m}{i} \right\rfloor$ 共有 $O(\sqrt n + \sqrt m)$ 种取值

证明：我们可以这样想象：假设有一根数轴表示的不同取值， $\left\lfloor \frac{n}{i} \right\rfloor$ 用 $O(\sqrt n)$ 个分割点把数轴分成了 $O(\sqrt n)$ 段（姑且叫它们段）， $\left\lfloor \frac{m}{i} \right\rfloor$ 用 $O(\sqrt m)$ 个分割点在不同的位置把数轴也分成了 $O(\sqrt m)$ 段（姑且叫它们段），现在把这些段overlay一下，你会发现现在总共有 $O(\sqrt n + \sqrt m)$ 个分割点把数轴分成了 $O(\sqrt n + \sqrt m)$ 段，每一段的取值可能不同（想一下，两个点只有不完全在相同的段和段时计算出来的 $\left\lfloor \frac{n}{i} \right\rfloor \left\lfloor \frac{m}{i} \right\rfloor$ 才可能不一样！）

结论：所在的一段结束于 $\min\{\lfloor n / \lfloor n / i \rfloor \rfloor, \lfloor m / \lfloor m / i \rfloor \rfloor\}$ 。

证明：类似一维情况。

HDU 2204 / 容斥原理

Posted on 2018-09-01

题面

大意：求在 [1, n] 中可以表示为 m^k 形式的数的个数（ $m, k \in \mathbb{Z}^+, k > 1,n \le 10^{18}$ ）

这道题不像是有 O(1) 的解法，但是哪怕是 $O(\sqrt{n})$ 的算法也会超时，因此直觉上我们应该寻找的是 $O(\log n)$ 的算法。

首先，因为 m = 1 的情形是trivial的，所以考虑 $m \ge 2$ 此时 $k \le \log_m n \le \log_2 n \le 60$ ，这也是的一个上界。

接下来我们发现：如果一个数可以表示为 m^4 的形式，也一定可以表示为 (m^2)^2 的形式，类似地，如果 k = pq 是合数的话，那么 m^k 一定可以表示为 (m^p)^q 和 (m^q)^p 。

令 P(i) 表示可以表示为 m^i 形式的数的个数，那么题目显然要我们求的是：

$\left| \bigcup_{素数p\le\log_2n}P(p) \right|$

接下来就用容斥原理很快就能解决这道题了，代码非常的简洁，注意精度：

#include <cstdio>
#include <cmath>

using namespace std;
typedef long long ll;
const double EPS = 1e-8;
const int p[] = { 0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59 };
ll n;

ll calc(ll d, ll a) {
    if (a > 60) return 0;
    if (d == 0) 
        return a == 1 ? 0 : (ll)(pow(n, 1.0 / a) + EPS) - 1;
    return calc(d - 1, a) - calc(d - 1, a * p[d]);
}

int main() {
    while (scanf("%lld", &n) != EOF)
        printf("%lld\n", -calc(17, 1) + 1);
    return 0;
}

NOIP模拟赛基站选址

Posted on 2018-09-01

题意

给定一个个顶点条边的无向图，第个节点有可靠度 r_i ，要求在节点中选取四个节点建设基站，其中有两个主基站，两个副基站。要求两个主基站必须和分别剩下的三个基站相连边，而两个副基站必须分别和两个主基站相连边。一个由两个主基站 a, b 和两个副基站 c, d 组成的基站的总可靠度定义为 (r_a + 1)(r_b + b) + r_cr_d ，求最大总可靠度。

$n \le 5 \times 10^4, m \le 2 \times 10^5$ 。

解法

题目要求我们寻找的四个点其实是三元环的并，主基站就是公共边的两个端点，而副基站就是公共边两侧的另外两个节点，因此，如果我们能够找出所有的三元环，并且对于每一条边维护他所在三元环相对顶点可靠度的最大和次大值，我们就可以通过枚举每一条边的方式完成最大总可靠度的计算。

因此，这个题目的核心就是寻找三元环，朴素的算法是 O(n^3) 的，明显不足以通过此题，因此我们要寻找复杂度更优的做法。

我们把顶点分为两类——大点和小点，分类的依据是点的度数，如果一个顶点的度数大于 $\sqrt m$ 就分在大点里，反之就落在小点里。我们对大点和小点分别处理。

对于大点，我们暴力枚举每个点：

因为对于无向图来说，有

$\sum_{i = 1}^n \deg(i) = 2m$

因此我们可以断言：大点的数量是 $O(\sqrt m)$ 的，因此枚举大点的时间复杂度是 $O((\sqrt m)^3) = O(m\sqrt m)$ 。

在这里我们处理了全部由大点构成的三元环。

对于小点，我们暴力枚举它的两条出边并判断对面两个点是否相连：

因为我们在枚举一个顶点的出边的时候，显然每条边都被枚举了 $\deg(u) - 1$ 次，所以图中小点发出的每一条边都被枚举了 $O(\sqrt m)$ 次，总共有 O(m) 条边，因此处理小点的时间复杂度还是 $O(m \sqrt m)$ 。

在这里我们处理了至少包含一个小点的三元环。

合并在一起，这个算法可以找到所有的三元环，总的时间复杂度为 $O(m \sqrt m)$ ，虽然对于稠密图而言因为 m = O(n^2) 因此算法的复杂度还是接近 O(n^3) ，但对于题目中的图大小还是非常优秀的。

代码实现也不是非常复杂，使用链式前向星存图，然后额外维护数组grt和snd表示一条边所在三元环中所对点可靠值最大和次大的点的编号，再首先一个哈希表用于快速查询点和点之间的相连关系就可以了，注意常数。

#include <cstdio>
#include <cmath>

using namespace std;
typedef long long ll;
const int N = 50100;
const int M = 500100;
const int HASH = 1000007;
const int INF = 0x3f3f3f3f;
int n, m;
int r[N];
int tot, fst[N], nxt[M], to[M];
int deg[N], bc, sc;
int big[N], small[N];
int grt[M], snd[M];

#define max(a, b) ((a) > (b) ? (a) : (b))
#define hashEdge(a, b) (((a) << 12) | (b))

int hfst[HASH], hnxt[M], hval[M], htot;
int hl[M], hr[M];

void link(int a, int b) {
    deg[a]++;
    nxt[++tot] = fst[a];
    fst[a] = tot; to[tot] = b;
    hl[++htot] = a; hr[htot] = b;
    hval[htot] = tot;
    int ent = hashEdge(a, b) % HASH;
    hnxt[htot] = hfst[ent];
    hfst[ent] = htot;
} 

inline int findEdge(int a, int b) { 
    int ent = hashEdge(a, b) % HASH;
    for (int i = hfst[ent]; i; i = hnxt[i])
        if (hl[i] == a && hr[i] == b) return hval[i];
    return -1; 
}

inline void updateEdge(int e, int c) {
    if (c == grt[e] || c == snd[e]) return;
    if (r[c] >= r[grt[e]]) snd[e] = grt[e], grt[e] = c;
    else if (r[c] >= r[snd[e]]) snd[e] = c;
}

inline void updateTriplet(int a, int b, int c, int ab, int ac, int bc) {
    updateEdge(ab, c);
    updateEdge(ab ^ 1, c);
    updateEdge(ac, b);
    updateEdge(ac ^ 1, b);
    updateEdge(bc, a);
    updateEdge(bc ^ 1, a);
}

int main() {    
    scanf("%d%d", &n, &m); tot = htot = 1; r[0] = -INF;
    for (int i = 1; i <= n; i++) scanf("%d", &r[i]);
    for (int i = 1; i <= m; i++) {
        int a, b;
        scanf("%d%d", &a, &b);
        link(a, b);
        link(b, a);
    }
    int s = sqrt(m);
    sc = bc = 0;
    for (int i = 1; i <= n; i++) {
        if (deg[i] < s) small[++sc] = i;
        else big[++bc] = i;
    }
    for (int i = 1; i <= sc; i++) {
        int w = small[i];
        for (int wu = fst[w]; wu; wu = nxt[wu]) {
            int u = to[wu];
            for (int wv = fst[w]; wv != wu; wv = nxt[wv]) {
                int v = to[wv];
                int uv = findEdge(u, v);
                if (uv > 0) updateTriplet(w, u, v, wu, wv, uv);
            }
        }
    }
    for (int i = 1; i <= bc; i++) {
        int w = big[i];
        for (int j = 1; j < i; j++) {
            int u = big[j];
            int wu = findEdge(w, u);
            if (wu < 0) continue;
            for (int k = 1; k < j; k++) {
                int v = big[k];
                int wv = findEdge(w, v), uv = findEdge(u, v);
                if (wv > 0 && uv > 0)
                    updateTriplet(w, u, v, wu, wv, uv);
            }
        }
    }
    int ans = 0; 
    for (int i = 2; i <= tot; i++) {
        int a = to[i], b = to[i ^ 1];
        int c = grt[i], d = snd[i];
        if (c == 0 || d == 0) continue;
        ans = max(ans, (r[a] + 1) * (r[b] + 1) + r[c] * r[d]);
    }
    printf("%d\n", ans);
    return 0;
}

NOIP模拟赛密码锁

Posted on 2018-09-01

题意

有把串级连接的锁（前一把的输入连接到后一把的输出），每一把锁将输入和一个操作数进行与，或，异或之一的二元位运算的结果作为输出，给定初始时每把锁的参数和运算类型。接下来个操作共两种：

查询某个数通过这一串锁的结果
修改某个锁的操作数和运算

$n ^5, m ^5 $，保证查询的数和操作数不超过 1000 。

解法

从数据范围可知正解的复杂度是 $O(n\log n)$ 级别的，考虑其操作中又有查询又有修改，直觉上我们断定这是一道数据结构题。

那维护什么呢？位运算？肯定是不行的，原因就是位运算是不结合的，比如说(a | b) & c就不一定等于a | (b & c)，而无论是线段树还是树状数组都只能维护有结合律的代数结构，因此直接维护位运算是行不通的。

但是因为二元位运算的一个操作数已经给出，而且每个位的运算有互不干扰，我们可以把给定的一个二元位运算转化为若干个依赖于输入某一位的一元位运算或函数。

那么这样做有什么好处呢？我们看看这样的一元函数有哪些：

：将输入位直接输出。
：将输入位取反输出。
：输出恒为0。
：输出恒为1。

定义它们之间的二元运算为串联： $(f \cdot g)(x)=g(f(x))$ ，则有如下性质：

结合律： $(f \cdot g) \cdot h = f \cdot (g \cdot h)$
单位元：
封闭性：任意两个运算的复合都是这四个运算之一（自己列张表）

可以看到这构成了一个幺半群，因此可以用线段树进行维护。接下来就很简单了，开棵线段树分别维护每一位，然后单点修改全局查询即可。

代码：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
#include <queue>
#define rg register

using namespace std;
typedef long long ll;
typedef unsigned int u32;
#define lch (rt << 1)
#define rch (lch | 1)
#define larg lch, l, mid
#define rarg rch, mid + 1, r
#define bit(a, b) (((a) >> (b)) & 1)
const int INF = 0x3f3f3f3f;
const int BITS = 10;
const int N = 200010;
const int T = N * 4;
const int r[4][4] = {
    { 0, 1, 1, 0 }, // 0
    { 0, 1, 0, 1 }, // 1
    { 0, 1, 3, 2 }, // N
    { 0, 1, 2, 3 } // E 
};
const int p[4][2] = {
    { 0, 0 },
    { 1, 1 },
    { 1, 0 },
    { 0, 1 }
};
const int x[3][2] = {
    { 0, 3 }, // AND
    { 3, 1 }, // OR
    { 3, 2 } // XOR
};

char t[BITS][T], a[BITS][N];
char s[5];
int n, m;

inline int read() {
    char ch = getchar(), s = 1;
    while (!isdigit(ch)) {
        if (ch == '-') s = -1;
        ch = getchar();
    }
    int ans = 0;
    while (isdigit(ch)) {
        ans = ans * 10 + ch - 48;
        ch = getchar();
    }
    return ans * s;
}

inline void pushup(int z, int rt) {
    t[z][rt] = r[t[z][lch]][t[z][rch]];
}

void build(int z, int rt, int l, int r) {
    if (l == r) {
        t[z][rt] = a[z][l];
        return;
    }
    int mid = (l + r) >> 1;
    build(z, larg);
    build(z, rarg);
    pushup(z, rt);
}

void update(int z, int rt, int l, int r, int pos, int d) {
    if (l == r) {
        t[z][rt] = d;
        return;
    }
    int mid = (l + r) >> 1;
    if (pos <= mid) update(z, larg, pos, d);
    else update(z, rarg, pos, d);
    pushup(z, rt);
}

int main() {
    freopen("cipher.in", "r", stdin);
    freopen("cipher.out", "w", stdout);
    scanf("%d%d", &n, &m);
    for (rg int i = 1; i <= n; ++i) {
        int tmp, op;
        scanf("%s%d", s, &tmp);
        if (s[0] == 'A') op = 0;
        else if (s[0] == 'O') op = 1;
        else op = 2;
        for (rg int j = 0; j < BITS; ++j)
            a[j][i] = x[op][bit(tmp, j)];
    }
    for (rg int i = 0; i < BITS; ++i)
        build(i, 1, 1, n);
    while (m--) {
        int op; scanf("%d", &op);
        if (op == 1) {
            int tmp; scanf("%d", &tmp);
            int res = 0;
            for (rg int i = 0; i < BITS; ++i)
                res |= p[t[i][1]][bit(tmp, i)] << i;
            printf("%d\n", res);
        } else {
            int pos, tmp;
            scanf("%d%s%d", &pos, s, &tmp);
            if (s[0] == 'A') op = 0;
            else if (s[0] == 'O') op = 1;
            else op = 2;
            for (rg int i = 0; i < BITS; ++i)
                update(i, 1, 1, n, pos, x[op][bit(tmp, i)]);
        }
    }
    fclose(stdin);
    fclose(stdout);
    return 0;
}

NOIP模拟赛城市规划

Posted on 2018-09-01

题意

有座高度为正整数的建筑物，定义不美观度为相邻建筑物高度差绝对值的最大值，你可以任意改动至多座建筑物的高度，求能够达到的最小不美观度。

$k \le n \le 2000, h \le 2 \times 10^9$

解法

从数据范围看出这应该是一道DP题，所以很容易陷入的一条歪路就是，定义 f[i, j] 为前座建筑物修改了次的最小不美观度，然后会发现没办法消除后效性~~当场去世~~。

但是如果我们把题目paraphrase一下，变成：在修改至多座建筑物高度的情况下，最大的相邻高度差最小是多少？首先我们发现随着修改数增加，不美观度肯定是可以单调递减的，其次我们从这熟悉的句型就发现这其实是一道二分答案的题目，问题就变成了如何判定一个不美观度能不能达到了，我们还是需要使用DP来完成。

定义 f[i] 为修改前 i - 1 座建筑物但是不修改的情况下为了达到我们二分的答案所需要的最小修改数，我们把它们全部初始化为 i - 1 ，即最坏情况下前面的我们都需要修改。

接下来我们思考状态转移，设我们此时二分的答案为，我们发现如果对于 j < i ，如果 $|h_j - h_i| \le (i - j)x$ ，那么就说明我们总可以通过调整和之间（不包括 i, j ）的 i - j - 1 座建筑物达到我们需要二分的答案。因此状态转移方程为：

$f[i] = \min_{1 \le j < i, |h_j - h_i| \le (i - j)x} \{ f[j] + i - j - 1\}$

那么最后总体上为了达到二分的答案所需要的最小修改数是多少呢？是 $\min_{1 \le i \le n} \{f[i]\}$ 吗？其实并不是，因为对于

，最坏情况下我们还需要修改

后面的

座建筑物，因此总的最小修改数为：

$\min_{1 \le i \le n} \{f[i] + n - i\}$

接下来就可以愉快的二分答案了，总的复杂度为 $O(n^2\log h_{max})$ ，还有要注意的是如果你采用int mid = (l + r) / 2的话是会爆int的，解决方法是使用int mid = l + (r - l) / 2：

#include <cstdio>

using namespace std;
const int INF = 0x7fffffff;
const int N = 2010;
int f[N], h[N];
int n, k;

#define abs(a) ((a) > 0 ? (a) : -(a))
#define min(a, b) ((a) < (b) ? (a) : (b))

int valid(int x) {
    for (int i = 1; i <= n; i++)
        f[i] = i - 1;
    for (int i = 2; i <= n; i++) {
        for (int j = 1; j < i; j++)
            if (abs(h[i] - h[j]) <= 1ll * x * (i - j))
                f[i] = min(f[i], f[j] + i - j - 1);
    }
    int ans = INF;
    for (int i = 1; i <= n; i++)
        ans = min(ans, f[i] + n - i);
    return ans;
} 

int main() {
    scanf("%d%d", &n, &k);
    for (int i = 1; i <= n; i++) scanf("%d", &h[i]);
    int l = 0, r = INF;
    while (l < r) {
        int mid = l + (r - l) / 2;
        if (valid(mid) <= k) r = mid;
        else l = mid + 1;
    }
    printf("%d\n", l);
    return 0;
}

最大多段连续子段和

Posted on 2018-09-01

问题概要

给定一个个元素的环形序列（号与号元素相邻），元素有正有负，求段连续的子段，使子段和最大，保证至少有个正数。

$n \le 3 \times 10^5, m \le 10^5$ 。

暴力

破环成链，类似序列上的最大连续子段和，定义 dp[i][j] 表示以为右端点选了段的最大子段和，口胡状态转移：

$dp[i][j] = a[i] + \max \begin{cases} dp[i - 1][j] \\ \max_{1 \le k < i} dp[k][j - 1] \end{cases}$

后面的 $\max$ 随便搞搞搞成

的，然后整个状态转移就是

的，加上破环成链的代价，总的时间复杂度为

，实测拿50分。

$O(n\log n)$ 正解

我们发现：既然我们要选一个正数，不如把这个正数所在的一整段正数全部选了；类似地，既然我们要不选一个负数，不如把这个负数所在的一整段负数排除掉。

于是乎我们就可以把整个环缩成若干段正负交替的“段”，假设我们缩出来是段正的夹着段负的，那么分情况讨论：

$p \le m$ ，直接选所有的正数就行了，注意我们可以把这里的一段正数断成两段或更多来选，因此无论如何都能选够段，此时明显最优。
，我们不能每一段都选了，但是我们可以假设我们先选了段，然后通过或消除，或合并的方式进行“调整”。

考虑调整的方式：

选一段负数，这等价于把两边的正数段合并了，总段数，对答案有负贡献。
不选一段正数，这等价于把两边的负数段合并了，总段数，对答案也有负贡献。

注意：我们不能先选一段负数，然后再不选其相邻的正数（或反过来），因为这样不仅血亏，对总段数也没有影响。

因为调整一次总段数减少1，以上的“调整”操作要进行 p - m 次。

既然对答案的贡献都是负的，我们把每一段的“值”变成该段数和的绝对值，然后问题就被转化为了：选择个不相邻的数，使它们的和（对答案的负贡献）最小。

这又是另一个经典问题了，参见CTSC 2007 数据备份，我们维护一个小根堆，每次取出堆顶元素，并在答案上减去其的贡献，然后执行“合并”，即如果这个数的前驱是，本身是，右边是，我们删除，然后插入一个新的权值为的等效节点。这样的话，如果我们下次再选到这个节点，对答案的影响等价于我们选了而不选。这在思想上类似于网络流的退流。为了维护前驱/后继关系，我们同时需要引入一个环形链表，总的复杂度约为 $O(n\log n)$ 。

我们发现这样做环形和非环形的差别就很小了，甚至是环形的更好写（不需要处理边界情况）。

代码：（注意因为人懒所以用了priority_queue~~饮鸩止渴~~，所以用的是懒惰删除法，还有就是代码中是大根堆维护负绝对值，本质上是一样的）

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
#include <queue>
#define register

using namespace std;
typedef long long ll;
typedef unsigned int u32;
const int INF = 0x3f3f3f3f;
const int N = 600100;
int n, m, ans, cnt;
int a[N], s[N], tot = 0;
int pre[N], nxt[N];
bool inv[N];

struct cmp {
    bool operator ()(int a, int b) {
        return s[a] < s[b];
    }
};
priority_queue<int, vector<int>, cmp> q;

#define max(a, b) ((a) > (b) ? (a) : (b))
#define abs(a) ((a) > 0 ? (a) : -(a))
#define sgn(a) ((a) >= 0 ? 1 : -1)

void merge(int x) {
    ans += s[x];
    s[++cnt] = s[pre[x]] + s[nxt[x]] - s[x];
    inv[pre[x]] = inv[nxt[x]] = true;
    pre[cnt] = pre[pre[x]];
    nxt[pre[cnt]] = cnt;
    nxt[cnt] = nxt[nxt[x]];
    pre[nxt[cnt]] = cnt;
    q.push(cnt);
}

void solve() {
    ans = 0;
    for (int i = 1; i <= n; ) {
        s[++tot] = 0;
        if (a[i] >= 0) while (i <= n && a[i] >= 0) s[tot] += a[i++];
        else while (i <= n && a[i] < 0) s[tot] += a[i++];
        if (s[i] >= 0) ans += s[i]; 
    }
    if (sgn(s[1]) == sgn(s[tot])) s[1] += s[tot--];
    if (tot / 2 <= m) return;
    for (int i = 1; i <= tot; i++) {
        pre[i] = i == 1 ? tot : i - 1;
        nxt[i] = i == tot ? 1 : i + 1;
        s[i] = -abs(s[i]);
        q.push(i);
    }
    cnt = tot;
    int left = tot / 2 - m;
    while (left--) {
        int now = q.top(); q.pop();
        if (inv[now]) left++;
        else merge(now);
    }
}

int main() {
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= n; ++i) scanf("%d", &a[i]);
    solve();
    printf("%d\n", ans);
    return 0;
}

线性筛卡常技巧

Posted on 2018-09-01

啊？卡常，没有啊？标程只是加了个小优化啊？

于是某人的欧拉筛就在 n = 10^8 的情况下被卡常死于非命。

素数的 $6k \pm 1$ 判定方法

任何正整数都可以表示为 6k + r 的形式，其中 k,r 为自然数且 $r \in [0, 6)$ ，我们讨论每一种：

：显然有，为合数。
：可能是素数。
：，在 $k \ge 1$ 时为合数。
：，在 $k \ge 1$ 时为合数。
：，为合数。
：可能是素数。

从上述分类讨论可知：除了 2, 3 以外，所有素数都可以表示为 $6k \pm 1$ 的形式，反之，如果 $n \bmod 6 \neq \pm 1$ ，我们就知道一定不是素数。

常数优化

在刚刚发现的 $6k\pm1$ 判定的指导下我们就可以对欧拉筛进行常数优化了：

void sieve() {
    tot = 0;
    for (int i = 1; i <= n; i++) v[i] = false;
    for (int i = 2; i <= 5; i++) {
        if (!v[i]) prime[++tot] = i;
        for (int j = 1; i * prime[j] <= n; j++) {
            v[i * prime[j]] = true;
            if (i % prime[j] == 0) break;
        }
    }
    for (int k = 6; k <= n; k += 6) {
        int i = k + 1;
        if (!v[i]) prime[++tot] = i;
        for (int j = 1; i * prime[j] <= n; j++) {
            v[i * prime[j]] = true;
            if (i % prime[j] == 0) break;
        }
        i += 4;
        if (!v[i]) prime[++tot] = i;
        for (int j = 1; i * prime[j] <= n; j++) {
            v[i * prime[j]] = true;
            if (i % prime[j] == 0) break;
        }
    }
}

可以看到我们手动处理了 [2,5] 的素数，然后使用了 $6k \pm 1$ 的结论，直接六个六个筛，但是这么筛会在一定程度上漏筛很多形式的合数，但因为这些形式的合数都是的倍数，在后面筛的循环中不受影响。

这个优化仅限于我们求取所有素数的情况，如果我们依赖于筛的其他副产物，如最小素因子，欧拉函数等，则该算法不正确。

实测中加了这个优化代码能快上2倍！

多重集组合

Posted on 2018-09-01

无限元素多重集组合

对于一个含有种元素的多重集 $S = \{a_1 \cdot \infty, a_2 \cdot \infty, \cdots, a_k \cdot \infty \}$ ，从当中取出个元素进行组合，共有

$\binom{r + k - 1}{k - 1} = \binom{r + k - 1}{r}$

种取法，这个结论用隔板法易证：

$\overbrace{\underbrace{\times \cdots \times \times | \times \times \cdots | \times \times}_{r个元素\times}}^{k -1个隔板|}$

有限元素多重集组合

对于一个含有种元素共计个的多重集 $S = \{a_1 \cdot n_1, a_2 \cdot n_2, \cdots, a_k \cdot n_k \}$ ，从当中取出个元素进行组合，方案数为：

$\sum_{X \subseteq \{1, 2, 3, \cdots, k\}} (-1)^{|X|} \binom{r - \sum_{i \in X}(n_i + 1) + k - 1}{k - 1}$

证明：

如果没有数量限制那么答案显然是 $\binom{r + k - 1}{k - 1}$ ，但是由于数量有限，我们要排除 a_i 多于 n_i 的方案，令 P_i 表示第种元素选了至少 n_i + 1 个的方案集合。要构造这样的方案，我们只需要先选出 n_i + 1 个该种元素，再在剩下的集合中选出 r - n_i - 1 个即可，于是我们有

$|P_i| = \binom{r - (n_i + 1) + k - 1}{k - 1}$

所有所有不合法的方案可以表示为

$\bigcup_{1 \le i \le k} P_i$

接下来结合容斥原理易证。

例子：CF451E

就是以上公式的应用，对于容斥原理使用递归的写法，求取组合数时使用了Lucas定理加以优化，代码很短：

#include <cstdio>

using namespace std;
typedef long long ll;
const ll P = 1000000007;
const int N = 22;
ll n, m, inv[N], a[N];

ll binom(ll a, ll b) {
    if (a < 0 || b < 0 || a < b) return 0;
    a %= P;
    if (a == 0 || b == 0) return 1;
    ll ret = 1;
    for (ll i = a; i > a - b; i--)
        ret = ret * i % P;
    for (ll i = 1; i <= b; i++)
        ret = ret * inv[i] % P;
    return ret;
}

ll calc(int d, ll x) {
    if (d == 0) return x == 0 ? 0 : binom(n + m - 1 - x, n - 1);
    else return (calc(d - 1, x) - calc(d - 1, x + a[d] + 1)) % P;
}

void setup() {
    inv[1] = 1;
    for (int i = 2; i <= n; i++)
        inv[i] = (-(P / i) * inv[P % i] % P + P) % P;
}

int main() {
    scanf("%lld%lld", &n, &m);
    setup();
    for (int i = 1; i <= n; i++)
        scanf("%lld", &a[i]);
    ll ans = binom(n + m - 1, n - 1) + calc(n, 0);
    printf("%lld\n", (ans % P + P) % P);
    return 0;
}

优化埃氏筛的时间复杂度分析

Posted on 2018-09-01 Edited on 2022-06-21

对于优化后的埃氏筛法：

bool v[N];

void sieve(int n) {
    memset(v, 0, sizeof(v));
    for (int i = 2; i <= n; i++) {
        if (v[i]) continue;
        for (int j = i; j <= n / i; j++)
            v[i * j] = true;
    }
}

简单分析一下它的流程，我们发现当i不是素数的时候内层循环才会执行，而内层循环执行的次数约为 $\frac{n}{i} - i$ 次，暗示了其实对于 $i > \sqrt{n}$ ，这个循环也不会被执行，这已经是埃氏筛的一个小优化：即一个数只会被它前半部分的素数筛去，那么总体的时间复杂度就可以表示为：

$\begin{aligned} \sum_{素数p \le \sqrt{n}} \left(\frac{n}{p} - p\right) &= \sum_{素数p \le \sqrt{n}} \frac{n}{p} - \sum_{素数p \le \sqrt{n}} p \\ &= n \sum_{素数p \le \sqrt{n}} \frac{1}{p} - \sum_{素数p \le \sqrt{n}} p \end{aligned}$

方便起见我们假设 $\sqrt{n}$ 为整数，这样就可以避免写 $\left\lfloor \sqrt{n} \right\rfloor$ ，我们先着手估计后半部分，即 $\sum_{素数p \le \sqrt{n}} p$ 的值，通过计算贡献的办法我们得到：

$\sum_{素数p \le \sqrt{n}} p = \sqrt{n}\pi\left(\sqrt{n}\right) - \sum_{i = 2}^{\sqrt{n} - 1}\pi(i)$

根据素数定理，我们有 $\pi(x) \approx \frac{x}{\ln x}$ ，那么我们可以得到：

$\begin{aligned} \sqrt{n}\pi\left(\sqrt{n}\right) - \sum_{i = 2}^{\sqrt{n} - 1}\pi(i) &\approx \frac{n}{\ln \sqrt{n}} - \sum_{i = 2}^{\sqrt{n} - 1}\frac{i}{\ln i} \\ &= \frac{2n}{\ln n} - \sum_{i = 2}^{\sqrt{n} - 1}\frac{i}{\ln i} \\ &\approx \frac{2n}{\ln n} - \int_2^{\sqrt{n}}\frac{x}{\ln x}dx \\ &< \frac{2n}{\ln n} - \frac{n}{\ln n} \\ &= O\left(\frac{n}{\ln n}\right) \end{aligned}$

对于前半部分，我们有Mertens Theorem指出：

$\lim_{n \to \infty} \left( \sum_{素数p \le n}\frac{1}{p} - \ln\ln n - M \right) = 0$

其中 $M \approx 0.2614972128$ 为Meissel-Mertens’ Constant。

这告诉我们前半部分的复杂度是 $O\left(n\ln\ln\sqrt{n}\right) = O (n\ln\ln n)$ 的，因此合并起来，整个算法的复杂度为 $O\left(n\ln\ln n - \frac{n}{\ln n}\right) = O(n \ln\ln n)$ ，和优化前的复杂度一致。

素数的判定

试除法

费马判定法

Miller-Rabin 判定法

自己的判定法？

计算的所有素数

埃拉托色尼筛法

欧拉筛

素因数分解

试除法

Pollard's-Rho方法

一维情况

二维情况

题意

解法

题意

解法

题意

解法

问题概要

暴力

正解

素数的判定方法

常数优化

无限元素多重集组合

有限元素多重集组合

例子：CF451E

$O(n\log n)$ 正解

素数的 $6k \pm 1$ 判定方法