作者: KSkun

[BZOJ3683]Falsita 题解

[BZOJ3683]Falsita 题解

题目地址:BZOJ:Problem 3683. — Falsita

题目描述

描述
到海边了呢……
如果没有那次选择,现在是不是会好些呢……
都过去了。
仰望着星空,迎面吹过一阵阵海风,倚靠着护栏,Fine 在海边静静地伫立着,在一个个无际的长夜后,Fine 终于放下了往事的痛楚,得到了治愈。
但是作为 Fine 的另一重人格的 Falsita 就没那么幸运了。她仍然被各种繁忙的事务困扰着。
虽然同在一副躯体中,Fine 与 Falsita 的精神世界却相差甚远,Fine 可以轻易地构造出幻梦时,Falsita 却只能停留在现实的痛楚中。
但是为了生活需要,她们还是需要经常达成共识。
让我们形式化的描述一下吧。
她们所在的精神世界是一棵以 1 号节点为根的树,每个树上的节点 u 都有一个权值Wu,她们每个人分别都在一个节点上,达成共识的方法就是两个人都到达一个共识节点(即到达它们的最近公共祖先)。
一个点 u 与另外一个点 v 之间想要达到共识需要花费的代价为Wu+Wv。
有时两人的精神有所触动时,有时点的权值会改变成某个数,有时以某个点的子树中的所有点的权值会加上某个数。
Falsita 和 Fine 经常需要达成共识,每一次询问,已知达成的共识节点,求她们花费的期望代价。

题意简述

有一棵有点权的树,以1为根。有三种操作:

  1. 单点点权加一个数
  2. 子树点权加一个数
  3. 求随机从$u$子树中选择两个点$(i, j)$且它们的LCA为$u$,$w_i+w_j$的期望

输入输出格式

输入格式:
输入共 m + 3 行。
第一行两个整数 n, m ,表示节点个数和操作数。
第二行 n – 1 个整数Pi,表示节点 i ( i = 2 . . . n ) 的父亲节点的编号。
第三行 n 个整数Wi。
接下来 m 行,每行表示一个操作。
1. S u delta 表示将节点 u 的权值加上 delta 。
2. M u delta 表示将以节点 u 为根的子树中的所有节点的权值加上 delta。
3. Q u 表示询问共识节点为 u 时的答案。
询问保证 u 不是叶子节点。

输出格式:
对于每组询问,输出答案,答案精确到小数点后 6 位。
你的程序输出的答案需要与标准答案之差不超过10^(-5)。

输入输出样例

输入样例#1:

4 6
1 2 2 
0 -6 3 0 
S 2 -5
M 3 8
S 1 -1
M 4 7
M 3 2
Q 1

输出样例#1:

2.000000

说明

前5个操作后,四个节点的权值分别为-1,-11,13,7。最近公共祖先为1的点对有(1,2),(1,3),(1,4),权值和分别为-12,12,6,故答案为(-12+12+6)/3=2。
对于 100% 的数据,1 ≤ n, m ≤ 300000, 0≤ |delta|, |Wi|≤ 20000。

题解

首先,期望可以表示为所有可能情况的和/情况数。对于确定的点$u$,它子树中所有可能情况的数量是可以预处理出来的。这个答案就是$u$子树中不同时在$u$的某个儿子的子树内的点对数,我们可以DFS树的使用利用子树大小来计算,假如计算$u$的方案数,$u$有儿子$v$,在遍历到$v$的时候,只需要对方案数加上$siz[v] \cdot (已经遍历过的子树大小和+1)$即可。
接下来考虑期望中分子的问题。如果不带修改,我们可以DFS的时候计算出所有节点的答案。如果用$sum[u]$表示$u$子树内所有节点的权值和,那么$u$有儿子$v$,$v$子树对$u$的答案贡献是$sum[v] \cdot (siz[u] – siz[v])$。
接下来我们考虑修改的问题,单点修改时,对该点答案的影响是$delta \cdot (siz[u] – 1)$,若该点有祖先$anc$,该点到该祖先路径上,祖先的儿子节点是$anc’$,那么对于该点的该祖先的影响是$delta \cdot (siz[anc] – siz[anc’])$。我们考虑进行树剖,在线段树上来维护这些个影响。显然,我们可以维护一个标记$delta \cdot (siz[u] – siz[v])$,其中$v$是$u$在重链上的儿子,这样,单点修改$u$就可以从$u$的父亲开始跳重链改了。但是当我们发现,一条重链最深的节点的重儿子可能不是修改点到根路径上该点的儿子,因此最深节点的情况要特别处理。我们直接把这个特别处理的修改改到该点处就好了。
对于子树加的情况,我们发现,只需要对子树里所有的节点的分子加上$delta \cdot 分母$就好了。而对于子树根的祖先们,其实跟单点修改产生的影响相同,只不过这里的影响是$delta \cdot siz[u]$。
总之,我们树剖以后在线段树上把这些个标记维护好就能实现正确的复杂度$O(n \log^2 n)$。

代码

// Code by KSkun, 2018/7
#include <cstdio>
#include <cctype>
#include <cstring>

#include <algorithm>
#include <vector>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2)
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1; register char c = fgc();
    for(; !isdigit(c); c = fgc()) if(c == '-') neg = -1;
    for(; isdigit(c); c = fgc()) res = (res << 1) + (res << 3) + c - '0';
    return res * neg;
}

inline bool isop(char c) {
    return c == 'S' || c == 'M' || c == 'Q';
}

inline char readop() {
    char c;
    while(!isop(c = fgc())) {}
    return c;
}

const int MAXN = 300005;

int n, m, w[MAXN];
std::vector<int> gra[MAXN];

int fa[MAXN], siz[MAXN], dep[MAXN], son[MAXN], dfn[MAXN], vn[MAXN], top[MAXN], clk;
LL ans1[MAXN], ans2[MAXN], sum[MAXN];

inline void dfs1(int u) {
    siz[u] = 1; sum[u] = w[u];
    for(int i = 0; i < gra[u].size(); i++) {
        int v = gra[u][i];
        if(v == fa[u]) continue;
        fa[v] = u; dep[v] = dep[u] + 1;
        dfs1(v);
        ans2[u] += 1ll * siz[v] * siz[u];
        siz[u] += siz[v];
        sum[u] += sum[v];
        if(siz[v] > siz[son[u]]) son[u] = v;
    }
}

inline void dfs2(int u, int tp) {
    dfn[u] = ++clk; vn[dfn[u]] = u; top[u] = tp;
    ans1[u] += 1ll * w[u] * (siz[u] - 1);
    if(son[u]) {
        ans1[u] += 1ll * (siz[u] - siz[son[u]]) * sum[son[u]];
        dfs2(son[u], tp);
    }
    for(int i = 0; i < gra[u].size(); i++) {
        int v = gra[u][i];
        if(v == fa[u] || v == son[u]) continue;
        ans1[u] += 1ll * (siz[u] - siz[v]) * sum[v];
        dfs2(v, v);
    }
}

#define lch o << 1
#define rch o << 1 | 1
#define mid ((l + r) >> 1)

LL tag1[MAXN << 2], tag2[MAXN << 2];
// 1: ptadd 2: subtadd

inline void pushdown(int o, int l, int r) {
    if(tag1[o]) {
        tag1[lch] += tag1[o]; tag1[rch] += tag1[o];
        if(l == mid) ans1[vn[l]] += tag1[o] * ans2[vn[l]];
        if(mid + 1 == r) ans1[vn[r]] += tag1[o] * ans2[vn[r]];
        tag1[o] = 0;
    }
    if(tag2[o]) {
        tag2[lch] += tag2[o]; tag2[rch] += tag2[o];
        if(l == mid) ans1[vn[l]] += (siz[vn[l]] - siz[son[vn[l]]]) * tag2[o];
        if(mid + 1 == r) ans1[vn[r]] += (siz[vn[r]] - siz[son[vn[r]]]) * tag2[o];
        tag2[o] = 0;
    }
}

inline void modify(int o, int l, int r, int ll, int rr, LL v, int type) {
    if(l >= ll && r <= rr) {
        if(type == 1) {
            if(l == r) ans1[vn[l]] += v * ans2[vn[l]];
            else tag1[o] += v;
        } else if(type == 2) {
            if(l == r) ans1[vn[l]] += (siz[vn[l]] - siz[son[vn[l]]]) * v;
            else tag2[o] += v;
        }
        return;
    }
    pushdown(o, l, r);
    if(ll <= mid) modify(lch, l, mid, ll, rr, v, type);
    if(rr > mid) modify(rch, mid + 1, r, ll, rr, v, type);
}

inline void update(int o, int l, int r, int x) {
    if(l == r) return;
    pushdown(o, l, r);
    if(x <= mid) update(lch, l, mid, x);
    else update(rch, mid + 1, r, x);
}

inline void modify(int u, int _lu, LL v) {
    int tu = top[u], lu = _lu;
    while(u) {
        if(dfn[tu] != dfn[u]) modify(1, 1, n, dfn[tu], dfn[u] - 1, v, 2);
        ans1[u] += v * (siz[u] - siz[lu]);
        lu = u; u = fa[tu]; tu = top[u];
    }
}

int main() {
    n = readint(); m = readint();
    for(int i = 2, p; i <= n; i++) {
        p = readint();
        gra[i].push_back(p); gra[p].push_back(i);
    }
    for(int i = 1; i <= n; i++) {
        w[i] = readint();
    }
    dfs1(1);
    dfs2(1, 1);
    char op; int u, delta;
    while(m--) {
        op = readop(); u = readint();
        if(op != 'Q') delta = readint();
        if(op == 'S') {
            ans1[u] += 1ll * delta * (siz[u] - 1);
            modify(fa[u], u, delta);
        } else if(op == 'M') {
            modify(1, 1, n, dfn[u], dfn[u] + siz[u] - 1, delta * 2ll, 1);
            modify(fa[u], u, 1ll * delta * siz[u]);
        } else {
            update(1, 1, n, dfn[u]);
            printf("%.6lf\n", double(ans1[u]) / ans2[u]);
        }
    }
    return 0;
}
[ZJOI2007]粒子运动 题解

[ZJOI2007]粒子运动 题解

题目地址:洛谷:【P4724】[ZJOI2007]粒子运动 – 洛谷、BZOJ:Problem 1094. — [ZJOI2007]粒子运动

题目描述

阿Q博士正在观察一个圆形器皿中的粒子运动。不妨建立一个平面直角坐标系,圆形器皿的圆心坐标为(x0, y0
),半径为R。器皿中有若干个粒子,假设第i个粒子在时刻0的位置为(xi, yi),速度为(vxi,vyi)(注:这是一个速度向量,若没有发生碰撞,t时刻的位置应该是(xi + t * vxi, yi + t * vyi) )。假设所有粒子的运动互不干扰;若某个粒子在某个时刻碰到了器皿壁,将发生完全弹性碰撞,即速度方向按照碰撞点的切线镜面反射,且速度大小不变(如图)。认为碰撞是瞬间完成的。
1094 - [ZJOI2007]粒子运动 题解
尽管碰撞不会影响粒子的速率,但是粒子却会受到一定的伤害,所以若某一个粒子碰撞了k次器皿壁,那么在第k次碰撞时它便会消亡。
出于研究的需要,阿Q博士希望知道从时刻0到所有粒子都消亡这段时间内,所有粒子之间的最近距离是什么。你能帮助他么?

题意简述

圆形容器中有$n$个粒子,各给定初始位置和速度,粒子会与容器壁碰撞,碰撞时发生反射,如图,速度会改变方向但不改变大小。每个粒子最多能和容器壁碰撞$k$次,否则会消失。求所有时刻所有粒子两两距离的最小值。

输入输出格式

输入格式:
第一行包含三个实数,分别为x0, y0, R,即圆形器皿的圆心坐标及半径。第二行包含两个正整数N, k,分别表示粒子的总数与消亡碰撞次数。接下来N行每行四个实数,分别为xi, yi, vxi , vyi,保证(xi, yi)都在圆内且(vxi, vyi)非零。

输出格式:
仅包含一个实数,即所有粒子的历史最近距离,精确到小数点后三位。

输入输出样例

输入样例#1:

0 0 10
2 10
0 -5 0 1
5 0 1 0

输出样例#1:

7.071

说明

对于所有的数据,2 ≤N ≤100。1≤k ≤100。
请注意实数精度问题。

题解

没用到计算几何算法的计算几何题。
我们考虑,在一段连续时间内,研究两个粒子的距离,假如在这段时间内不发生碰撞,则距离可以表示为一个二次函数,假设第一个粒子的位置可以表示为$P+tV$($P$是起始位置,$V$是速度向量),第二个可以表示为$P’+tV’$,则可以列出距离的表达式$f(t) = \mathrm{length}[(P+tV) – (P’+tV’)]$,展开后得到一个二次函数。既然是二次函数,我们可以在$O(1)$时间内求得区间最小值。
因此,本题的大体思路是,枚举每两个粒子,在每一个两个粒子都不发生碰撞的连续区间内求区间最小距离,来更新最后的答案。这个复杂度是$O(n^2k)$的。
我用了对称的方法解决反射的速度变化,先把射线的起点对称至碰撞点处的半径另一侧,然后在另一侧构造新的速度向量,并处理长度即可。过程如图所示。
particle - [ZJOI2007]粒子运动 题解

代码

// Code by KSkun, 2018/7
#include <cstdio>
#include <cctype>
#include <cmath>

#include <algorithm>

const int MAXN = 105;
const double EPS = 1e-8, PI = acos(-1);

inline int dcmp(double x) {
    if(fabs(x) < EPS) return 0; else return x < 0 ? -1 : 1;
}

struct Vector {
    double x, y;
    Vector(double x = 0, double y = 0) : x(x), y(y) {}
    inline Vector operator+(Vector rhs) {
        return Vector(x + rhs.x, y + rhs.y);
    }
    inline Vector operator-(Vector rhs) {
        return Vector(x - rhs.x, y - rhs.y);
    }
    inline Vector operator*(double rhs) {
        return Vector(x * rhs, y * rhs);
    }
    inline Vector operator/(double rhs) {
        return Vector(x / rhs, y / rhs);
    }
};

typedef Vector Point;

inline double dot(Vector u, Vector v) {
    return u.x * v.x + u.y * v.y;
}

inline double cross(Vector u, Vector v) {
    return u.x * v.y - u.y * v.x;
}

inline double length(Vector v) {
    return sqrt(dot(v, v));
}

inline double angle(Vector u, Vector v) {
    return acos(dot(u, v) / length(u) / length(v));
}

inline Vector rotate(Vector v, double rad) {
    return Vector(v.x * cos(rad) - v.y * sin(rad), v.x * sin(rad) + v.y * cos(rad));
}

struct Line {
    Point p; Vector v;
    Line(Point p = Point(), Vector v = Vector()) : p(p), v(v) {}
};

struct Ray : public Line {
    Ray(Point p = Point(), Vector v = Vector()) : Line(p, v) {}
};

struct Circle {
    Point c; double r;
    Circle(Point c = Point(), double r = 0) : c(c), r(r) {}
};

inline Point intersection(Circle o, Ray l) {
    double x0 = o.c.x, y0 = o.c.y, a = l.p.x, b = l.p.y, c = l.v.x, d = l.v.y;
    double A = c * c + d * d, B = 2 * (a * c - c * x0 + b * d - d * y0), 
        C = a * a - 2 * a * x0 + x0 * x0 + b * b - 2 * b * y0 + y0 * y0 - o.r * o.r;
    double delta = B * B - 4 * A * C;
    double t1 = (-B + sqrt(delta)) / 2 / A, t2 = (-B - sqrt(delta)) / 2 / A;
    if(dcmp(t1) > 0) return l.p + l.v * t1; else return l.p + l.v * t2;
}

inline Point projection(Line l, Point p) {
    return l.p + l.v * (dot(l.v, p - l.p) / dot(l.v, l.v));
}

inline double distance(Point p1, Point p2) {
    return length(p1 - p2);
}

inline double getmin(double A, double B, double C, double l, double r) {
    double v1 = A * l * l + B * l + C, v2 = A * r * r + B * r + C, v3 = 1e10;
    double m = -B / 2 / A;
    if(m >= l && m <= r) v3 = A * m * m + B * m + C;
    return std::min(std::min(v1, v2), v3);
}

Circle o;
int n, k;
Point p[MAXN];
Vector v[MAXN];

inline double sqr(double x) {
    return x * x;
}

inline double getmin(int a, int b) {
    int ak = 0, bk = 0; double res = 1e10;
    Ray ra(p[a], v[a]), rb(p[b], v[b]);
    while(ak < k && bk < k) {
        Point ai = intersection(o, ra), bi = intersection(o, rb);
        double t1 = distance(ra.p, ai) / length(ra.v), t2 = distance(rb.p, bi) / length(rb.v);
        double nt = std::min(t1, t2);
        double A = sqr(ra.v.y - rb.v.y) + sqr(ra.v.x - rb.v.x),
            B = 2 * ((ra.p.y - rb.p.y) * (ra.v.y - rb.v.y) 
                + (ra.p.x - rb.p.x) * (ra.v.x - rb.v.x)),
            C = sqr(ra.p.y - rb.p.y) + sqr(ra.p.x - rb.p.x);
        res = std::min(res, getmin(A, B, C, 0, nt));
        if(dcmp(nt - t1) == 0) {
            ak++; 
            Point proj = projection(Line(o.c, o.c - ai), ra.p);
            Point p1(proj.x * 2 - ra.p.x, proj.y * 2 - ra.p.y);
            Vector v1 = p1 - ai; v1 = v1 / length(v1) * length(ra.v);
            ra.p = ai; ra.v = v1;
        } else {
            ra.p = ra.p + ra.v * nt;
        }
        if(dcmp(nt - t2) == 0) {
            bk++; 
            Point proj = projection(Line(o.c, o.c - bi), rb.p);
            Point p1(proj.x * 2 - rb.p.x, proj.y * 2 - rb.p.y);
            Vector v1 = p1 - bi; v1 = v1 / length(v1) * length(rb.v);
            rb.p = bi; rb.v = v1;
        } else {
            rb.p = rb.p + rb.v * nt;
        }
    }
    return res;
}

int main() {
    double x0, y0, r;
    scanf("%lf%lf%lf", &x0, &y0, &r);
    o = Circle(Point(x0, y0), r);
    scanf("%d%d", &n, &k);
    for(int i = 1; i <= n; i++) {
        double xi, yi, vxi, vyi;
        scanf("%lf%lf%lf%lf", &xi, &yi, &vxi, &vyi);
        p[i] = Point(xi, yi); v[i] = Vector(vxi, vyi);
    }
    double ans = 1e10;
    for(int i = 1; i <= n; i++) {
        for(int j = i + 1; j <= n; j++) {
            ans = std::min(ans, getmin(i, j));
        }
    }
    printf("%.3lf", dcmp(ans) > 0 ? sqrt(ans) : 0);
    return 0;
}
[NOI2014]随机数生成器 题解

[NOI2014]随机数生成器 题解

题目地址:洛谷:【P2354】[NOI2014]随机数生成器 – 洛谷、BZOJ:Problem 3671. — [Noi2014]随机数生成器

题目描述

小 H 最近在研究随机算法。随机算法往往需要通过调用随机数生成函数(例如 Pascal 中的 random 和 C/C++中的 rand)来获得随机性。事实上,随机数生成函数也并不是真正的“随机”,其一般都是利用某个算法计算得来的。
比如,下面这个二次多项式递推算法就是一个常用算法:
算法选定非负整数 x0,a,b,c,d 作为随机种子,并采用如下递推公式进行计算。
对于任意 i ≥ 1, xi=(a*x[i-1]^2+b*x[i-1]+c)mod d 这样可以得到一个任意长度的非负整数数列{xi},i≥1,一般来说,我们认为这个数列是随机的。
利用随机序列{xi},i≥1,我们还可以采用如下算法来产生一个 1 到 K 的随机排列{Ti},i=1 to k:

  1. 初始设 T 为 1 到 K 的递增序列;
  2. 对 T 进行 K 次交换,第 i 次交换,交换 Ti 和 T[xi mod i + 1] 的值。

此外,小 H 在这 K 次交换的基础上,又额外进行了 Q 次交换操作,对于第i 次额外交换,小 H 会选定两个下标 ui 和 vi,并交换 T[ui] 和 T[vi] 的值。
为了检验这个随机排列生成算法的实用性,小 H 设计了如下问题:
小 H 有一个 N 行 M 列的棋盘,她首先按照上述过程,通过 N × M + Q 次交换操作,生成了一个 1~N × M 的随机排列 {Ti},i=1 to N*M,然后将这 N × M 个数逐行逐列依次填入这个棋盘:也就是第 i 行第 i 列的格子上所填入的数应为 T[(i-1)*M+uj]。
接着小 H 希望从棋盘的左上角,也就是第一行第一列的格子出发,每次向右走或者向下走,在不走出棋盘的前提下,走到棋盘的右下角,也就是第 N 行第M 列的格子。
小 H 把所经过格子上的数字都记录了下来,并从小到大排序,这样,对于任何一条合法的移动路径,小 H 都可以得到一个长度为 N + M − 1 的升序序列,我们称之为路径序列。
小 H 想知道,她可能得到的字典序最小的路径序列应该是怎样的呢?

题意简述

我们有一个随机序列生成算法,对于给定的$x_0, a, b, c, d$,有$x_i = (ax_{i-1}^2 + bx_{i-1} + c) \bmod d$。我们可以利用序列$x_i$来生成随机排列$T_i$,对于一个长为$k$的顺序排列$T_i$,我们进行$k$次交换,每次交换$T_i$和$T_{x_i \bmod i + 1}$。
我们有一个$n \times m$的矩阵,现在,我们按照上面的方法生成一个$n \times m$的随机排列,另外,给出$q$次额外交换,每次交换排列中$T_{u_i}$和$T_{v_i}$两个位置。然后,把这个排列放进矩阵中,即,矩阵$A$的元素$A_{i, j}$上填入$T_{(i-1)m+j}$。
现在,我们需要找到一条从矩阵左上角$(1, 1)$出发,每次向右或向下走一格,到右下角$(n, m)$的路径,这条路径上经过的格子的数字构成了一个长度$n+m-1$的序列,现在想求排序后字典序最小的这样的序列。

输入输出格式

输入格式:
第1行包含5个整数,依次为 x_0,a,b,c,d ,描述小H采用的随机数生成算法所需的随机种子。
第2行包含三个整数 N,M,Q ,表示小H希望生成一个1到 N×M 的排列来填入她 N 行 M 列的棋盘,并且小H在初始的 N×M 次交换操作后,又进行了 Q 次额外的交换操作。
接下来 Q 行,第 i 行包含两个整数 u_i,v_i,表示第 i 次额外交换操作将交换 T_(u_i )和 T_(v_i ) 的值。

输出格式:
输出一行,包含 N+M-1 个由空格隔开的正整数,表示可以得到的字典序最小的路径序列。

输入输出样例

输入样例#1:

1 3 5 1 71 
3 4 3 
1 7 
9 9 
4 9 

输出样例#1:

1 2 6 8 9 12 

输入样例#2:

654321 209 111 23 70000001 
10 10 0 

输出样例#2:

1 3 7 10 14 15 16 21 23 30 44 52 55 70 72 88 94 95 97

输入样例#3:

123456 137 701 101 10000007 
20 20 0 

输出样例#3:

1 10 12 14 16 26 32 38 44 46 61 81 84 101 126 128 135 140 152 156 201 206 237 242 243 253 259 269 278 279 291 298 338 345 347 352 354 383 395 

说明

对于样例 1,根据输入的随机种子,小 H 所得到的前 12 个随机数xi为:

9 5 30 11 64 42 36 22 1 9 5 30

根据这 12 个随机数,小 H 在进行初始的 12 次交换操作后得到的排列为:

6 9 1 4 5 11 12 2 7 10 3 8

在进行额外的 3 次交换操作之后,小 H 得到的最终的随机排列为:

12 9 1 7 5 11 6 2 4 10 3 8

棋盘为

12 9 1 7 
5 11 6 2 
4 10 3 8

最优路径依次经过的数字为 :12-9-1-6-28。
random noi14 - [NOI2014]随机数生成器 题解

题解

略有些无聊的卡空间阅读题。
首先,我们按照题目的要求处理出$T_i$,然后,发现开三个$5000 \times 5000$的int数组显然会爆炸,于是我们考虑重复利用最开始的$x_i$数组来存储每个数在$T_i$中的下标。
显然我们可以从小到大贪心地选择每个可行的数字,选择一个数字会对它所在行之前行可选择的列的右端点左移至该数所在列,且之后行可选择的列的左端点右移至该数所在列,我们开两个$5000$大小的int数组来存每一行能选择的列区间即可。
复杂度$O(nm)$,卡常卡空间,尽量控制吧。

代码

// Code by KSkun, 2018/7
#include <cstdio>
#include <cctype>
#include <cstring>

#include <algorithm>

typedef long long LL;

inline int min(int a, int b) {
    return a < b ? a : b;
}

inline int max(int a, int b) {
    return a > b ? a : b;
}

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2)
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1; register char c = fgc();
    for(; !isdigit(c); c = fgc()) if(c == '-') neg = -1;
    for(; isdigit(c); c = fgc()) res = (res << 1) + (res << 3) + c - '0';
    return res * neg;
}

const int MAXN = 5005;

LL a, b, c, d;
int n, m, q, x[MAXN * MAXN], T[MAXN * MAXN], l[MAXN], r[MAXN];

int main() {
    x[0] = readint(); a = readint(); b = readint(); c = readint(); d = readint();
    n = readint(); m = readint(); q = readint();
    for(int i = 1; i <= n * m; i++) {
        x[i] = (1ll * a * x[i - 1] * x[i - 1] + 1ll * b * x[i - 1] + c) % d;
    }
    for(int i = 1; i <= n * m; i++) {
        T[i] = i;
    }
    for(int i = 1; i <= n * m; i++) {
        std::swap(T[i], T[x[i] % i + 1]);
    }
    while(q--) {
        int u = readint(); int v = readint();
        std::swap(T[u], T[v]);
    }
    for(int i = 1; i <= n * m; i++) {
        x[T[i]] = i;
    }
    for(int i = 1; i <= n; i++) {
        r[i] = m + 1;
    }
    int cnt = 0;
    for(int i = 1; i <= n * m; i++) {
        int X = x[i] / m + 1, Y = x[i] % m;
        if(!Y) { X--; Y += m; }
        if(Y >= l[X] && Y <= r[X]) {
            printf("%d ", i); cnt++;
            for(int j = 1; j <= n; j++) {
                if(j < X) r[j] = min(r[j], Y);
                if(j > X) l[j] = max(l[j], Y);
            }
        }
        if(cnt >= n * m - 1) break;
    }
    return 0;
}