SDOI 2017 Round 1 大赏

2020-03-06

考虑到写详解有点花时间, 不写又容易忘 = =, 那就写个大概吧.

SD 两轮省选 (可能?) 会比较科学一点吧

不要因为 Round 1 而错怪 SDOI.

「SDOI2017」数字表格

题意简述

求

$\prod_{i=1}^n \prod_{j=1}^m f( \gcd(i, j) ) \bmod 10^9 + 7$

其中 $f(n)$ 为斐波那契数列数列第 $n$ 项, $1 \le n,m \le 10^6$, 多组测试数据, $1 \le T \le 1000$.

题目链接

https://loj.ac/problem/2000

看到 $\gcd$ 不如来反演一波. 不妨令 $n \le m$.

枚举约数并整理, 有

$\prod_{i=1}^n \prod_{j=1}^m \prod_{d = 1}^n f(d) [\gcd(i, j) = d]$ $\prod_{d=1}^n f(d) ^ {\sum\limits_{i=1}^n \sum\limits_{j=1}^m{\gcd(i, j) = d}}$

由莫比乌斯反演的套路, 得

$\prod_{d=1}^n f(d) ^ {\sum\limits_{T=1}^{\lfloor \frac{n}{d} \rfloor} \mu (T) \lfloor \frac{n}{Td} \rfloor \lfloor \frac{m}{Td} \rfloor}$ $\prod_{T=1}^n (\prod_{d\mid T} f(d) ^ {\mu (\frac{T}{d})}) ^ {\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor}$

设 $g(T) = \prod_{d\mid T} f(d) ^ {\mu (\frac{T}{d})}$, 则所求即为

$\prod_{T=1}^n g(T) ^ {\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor}$

考虑计算 $g(T)$. 可以发现, 由于 $\mu(n)$ 取值只有 $0, -1, 1$, 计算 $g(T)$ 可以通过枚举倍数简单地实现. 时间复杂度 $O(n \ln n)$.

配合前缀积和数论分块即可, 注意快速幂时间复杂度 $O(\log b)$.

时间复杂度 $O(n \log P + n \ln n + T \cdot n \sqrt n \log P)$.

代码实现

> 点击显示 / 隐藏

// LOJ #2000
// DeP
#include <cctype>
#include <cstdio>
#include <algorithm>
using namespace std;

namespace IO {
    const int MAXSIZE = 1 << 18 | 1;
    char buf[MAXSIZE], *p1, *p2;

    inline int Gc() {
        return p1 == p2 &&
            (p2 = (p1 = buf) + fread(buf, 1, MAXSIZE, stdin), p1 == p2)? EOF: *p1++;
    }
    template<typename T> inline void read(T& x) {
        x = 0; int f = 0, ch = Gc();
        while (!isdigit(ch)) f |= ch == '-', ch = Gc();
        while (isdigit(ch)) x = x * 10 + ch - '0', ch = Gc();
        if (f) x = -x;
    }
}
using IO::read;

const int MAXN = 1e6+5, P = 1e9 + 7;

int fpow(int base, int b) {
    int ret = 1;
    while (b > 0) {
        if (b & 1) ret = 1LL * ret * base % P;
        base = 1LL * base * base % P, b >>= 1;
    }
    return ret;
}

bool notPrime[MAXN];
int Prime[MAXN], tot;
int mu[MAXN], invf[MAXN], f[MAXN], g[MAXN];

void EulerSieve() {
    notPrime[1] = true;
    g[0] = g[1] = invf[1] = f[1] = mu[1] = 1;
    for (int i = 2; i < MAXN; ++i) {
        g[i] = 1;
        f[i] = (f[i-1] + f[i-2]) % P;
        invf[i] = fpow(f[i], P-2);
        if (!notPrime[i]) Prime[++tot] = i, mu[i] = -1;
        for (int j = 1; j <= tot && i*Prime[j] < MAXN; ++j) {
            notPrime[i*Prime[j]] = true;
            if (i % Prime[j] == 0) { mu[i*Prime[j]] = 0; break; }
            mu[i*Prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i < MAXN; ++i) if (mu[i] != 0)
        for (int j = i; j < MAXN; j += i)
            g[j] = 1LL * g[j] * (mu[i] == 1? f[j / i]: invf[j / i]) % P;
    for (int i = 2; i < MAXN; ++i) g[i] = 1LL * g[i] * g[i-1] % P;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("input.in", "r", stdin);
#endif
    EulerSieve();
    int T; read(T);
    while (T--) {
        static int n, m;
        read(n), read(m);
        if (n > m) swap(n, m);
        int ans = 1;
        for (int R, L = 1; L <= n; L = R + 1) {
            R = min(n / (n / L), m / (m / L));
            int s = 1LL * g[R] * fpow(g[L-1], P-2) % P;
            ans = 1LL * ans * fpow(s, 1LL * (n / L) * (m / L) % (P-1)) % P;
        }
        printf("%d\n", ans);
    }
    return 0;
}

「SDOI2017」树点涂色

题目链接

https://loj.ac/problem/2001

解题思路

观察操作 1, 以及初始颜色都不相同, 可以发现每次把一条从根节点开始的链染成新颜色的操作类似与 LCT 中的 access. 另加考虑可以发现, 某个节点 $u$ 到根节点的颜色数, 就是在 LCT 上经过虚边的个数 + 1. (假设初始状态中, 树上的边在 LCT 上都是虚边)

根据这个思路, 操作 2 可以简单计算, 也就是 dist(u) + dist(v) - 2 * dist(LCA) + 1. (此处用 dist 指代当前节点到根的颜色数, 因为 LCA 被减两次所以答案 + 1)

那么操作 3 呢? 回想 LCT 维护子树信息的过程, 体现在这道题中也就是在更改 “实儿子” 和 “虚儿子” 子树的信息, 换言之, 子树内加减, 然后维护最值就好了. 可以计算出 DFS 序后使用线段树维护.

实现中在建树的时候直接使用深度作为初始最值, 以及 access 修改的时候需要 findroot 一下找到实际子树的根再修改.

时间复杂度 $O(m \log ^2 n)$, ~~最近一群毒瘤卡树剖, 不敢再有 log^3 了~~

代码实现

> 点击显示 / 隐藏

// LOJ #2001
// DeP
#include <cctype>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

namespace IO {
    const int MAXSIZE = 1 << 18 | 1;
    char buf[MAXSIZE], *p1, *p2;

    inline int Gc() {
        return p1 == p2 &&
            (p2 = (p1 = buf) + fread(buf, 1, MAXSIZE, stdin), p1 == p2)? EOF: *p1++;
    }
    template<typename T> inline void read(T& x) {
        x = 0; int f = 0, ch = Gc();
        while (!isdigit(ch)) f |= ch == '-', ch = Gc();
        while (isdigit(ch)) x = x * 10 + ch - '0', ch = Gc();
        if (f) x = -x;
    }
}
using IO::read;

const int MAXN = 1e5+5;

int n, m;

namespace Graph {
    struct Edge { int nxt, to; } edges[MAXN << 1];
    int head[MAXN], eidx;

    inline void init() { memset(head, -1, sizeof head), eidx = 1; }
    inline void AddEdge(int from, int to) {
        edges[++eidx] = (Edge){ head[from], to };
        head[from] = eidx;
    }
}

namespace LCT {
    int ch[2][MAXN], pre[MAXN];

    inline int which(const int& u) { return pre[u]? ch[1][pre[u]] == u: 0; }
    inline bool nroot(const int& u) {
        return pre[u]? ch[0][pre[u]] == u || ch[1][pre[u]] == u: false;
    }

    inline void rotate(const int& u) {
        int fa = pre[u], w = which(u);
        pre[u] = pre[fa];
        if (nroot(fa)) ch[which(fa)][pre[fa]] = u;
        ch[w][fa] = ch[w^1][u];
        if (ch[w^1][u]) pre[ch[w^1][u]] = fa;
        ch[w^1][u] = fa, pre[fa] = u;
    }

    inline void splay(const int& u) {
        while (nroot(u)) {
            int fa = pre[u];
            if (nroot(fa)) which(fa) == which(u)? rotate(fa): rotate(u);
            rotate(u);
        }
    }
}

namespace HLD {
    using namespace Graph;
    int size[MAXN], son[MAXN], depth[MAXN], pre[MAXN];
    int topfa[MAXN], dfn[MAXN], rnk[MAXN], dfs_clock;

    void dfs1(int u, int fa) {
        depth[u] = depth[fa] + 1;
        size[u] = 1, son[u] = -1, LCT::pre[u] = pre[u] = fa;
        for (int v, i = head[u]; ~i; i = edges[i].nxt) {
            if ((v = edges[i].to) == fa) continue;
            dfs1(v, u), size[u] += size[v];
            if (son[u] == -1 || size[v] > size[son[u]]) son[u] = v;
        }
    }

    void dfs2(int u, int top) {
        topfa[u] = top;
        dfn[u] = ++dfs_clock, rnk[dfs_clock] = u;
        if (~son[u]) dfs2(son[u], top);
        for (int v, i = head[u]; ~i; i = edges[i].nxt)
            if ((v = edges[i].to) != pre[u] && v != son[u]) dfs2(v, v);
    }

    inline int LCA(int u, int v) {
        while (topfa[u] != topfa[v]) {
            if (depth[topfa[u]] < depth[topfa[v]]) swap(u, v);
            u = pre[topfa[u]];
        }
        return depth[u] > depth[v]? v: u;
    }

    inline void solve(int root = 1) { dfs1(root, 0), dfs2(root, root); }
}

namespace SGT {
#define lc (nd<<1)
#define rc (nd<<1|1)
    using namespace HLD;
    int datMax[MAXN << 2], tagAdd[MAXN << 2];

    inline void maintain(const int& nd) {
        datMax[nd] = max(datMax[lc], datMax[rc]);
    }

    inline void pushdown(const int& nd) {
        if (!tagAdd[nd]) return;
        tagAdd[lc] += tagAdd[nd], datMax[lc] += tagAdd[nd];
        tagAdd[rc] += tagAdd[nd], datMax[rc] += tagAdd[nd];
        tagAdd[nd] = 0;
    }

    void build(int nd, int L, int R) {
        if (L == R) return void( datMax[nd] = depth[rnk[L]] );
        int Mid = (L + R) / 2;
        build(lc, L, Mid), build(rc, Mid+1, R);
        maintain(nd);
    }

    void Mdy(int nd, int L, int R, const int& opL, const int& opR, const int& val) {
        if (opL <= L && R <= opR) return datMax[nd] += val, void( tagAdd[nd] += val );
        int Mid = (L + R) / 2;
        pushdown(nd);
        if (opL <= Mid) Mdy(lc, L, Mid, opL, opR, val);
        if (opR > Mid) Mdy(rc, Mid+1, R, opL, opR, val);
        maintain(nd);
    }
    inline void Mdy(const int& u, const int& val) {
        return Mdy(1, 1, n, dfn[u], dfn[u] + size[u] - 1, val);
    }

    int Qry(int nd, int L, int R, const int& opL, const int& opR) {
        if (opL <= L && R <= opR) return datMax[nd];
        pushdown(nd);
        int Mid = (L + R) / 2, ret = 0;
        if (opL <= Mid) ret = max(ret, Qry(lc, L, Mid, opL, opR));
        if (opR > Mid) ret = max(ret, Qry(rc, Mid+1, R, opL, opR));
        return ret;
    }

    inline int Qry(const int& u) { return Qry(1, 1, n, dfn[u], dfn[u]); }
    inline int Sub(const int& u) { return Qry(1, 1, n, dfn[u], dfn[u] + size[u] - 1); }
#undef lc
#undef rc
}

namespace LCT {
    inline int findroot(int u) {
        while (ch[0][u]) u = ch[0][u];
        return u;
    }

    void access(int u) {
        for (int v = 0; u; v = u, u = pre[u]) {
            splay(u);
            if (ch[1][u]) SGT::Mdy(findroot(ch[1][u]), 1);
            if (v) SGT::Mdy(findroot(v), -1);
            ch[1][u] = v;
        }
    }
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("input.in", "r", stdin);
#endif
    // init
    Graph::init();
    // input
    read(n), read(m);
    for (int u, v, i = 1; i < n; ++i)
        read(u), read(v), Graph::AddEdge(u, v), Graph::AddEdge(v, u);
    // solve
    HLD::solve(), SGT::build(1, 1, n);
    // output
    while (m--) {
        static int opt, x, y;
        read(opt), read(x);
        switch (opt) {
            case 1: LCT::access(x); break;
            case 2: read(y), printf("%d\n", SGT::Qry(x) + SGT::Qry(y) - 2 * SGT::Qry(HLD::LCA(x, y)) + 1); break;
            case 3: printf("%d\n", SGT::Sub(x)); break;
            default: fprintf(stderr, "ERR\n");
        }
    }
    return 0;
}

「SDOI2017」序列计数

高一刚学矩阵快速幂的时候整天写斐波那契数列的 n 倍经验 = =

然后某次学长测试, 出了一个矩阵快速幂, 差点就 A 了, 可惜矩阵写反了…

再然后就没独立写出过矩阵快速幂的题了, 凄惨 = =

题目链接

https://loj.ac/problem/2002

解题思路

要求和是 $p$ 的倍数, 那么显然是模 $p$ 意义下为 $0$ 了. 考虑到至少有一个数为质数的限制可以通过减法原理解决, 即计算一遍所有数的方案数, 再减去只用合数的答案.

容易想到一个暴力 DP:

设 $f(i, j)$ 表示选 $i$ 个数, 加起来模 $p$ 为 $j$ 的方案数. 则

$f(i, j) = \sum_k f(i-1, (j - k)\bmod p)$

每次选择满足限制的 $k$ 就好了.

观察到 DP 中每次转移都是一样的, 考虑用矩阵快速幂优化.

容易想到 $O(mp)$ 的构造矩阵方法, 考虑在模 $p$ 意义下, 矩阵中某一列 (或者是某一行) 是循环的, 于是可以在 $O(m + p^2)$ 的时间内构造.

具体实现的时候有一些细节, 还是参考代码实现吧.

时间复杂度 $O(m + p^2 + p^3 \log n)$.

代码实现

这道题中, 两个矩阵相乘可以通过某些方法优化到 $O(p^2)$, 研究不懂, $O(p^3)$ 养老…

> 点击显示 / 隐藏

// LOJ #2002
// DeP
#include <cstdio>
#include <cstring>

const int MAXM = 2e7+5, MAXP = 105;
const int MOD = 20170408;

int n, m, p;
int f[MAXP], g[MAXP];

bool notPrime[MAXM];
int Prime[MAXM], tot;

void EulerSieve() {
    notPrime[1] = true;
    for (int i = 2; i <= m; ++i) {
        if (!notPrime[i]) Prime[++tot] = i;
        for (int j = 1; j <= tot && i*Prime[j] <= m; ++j) {
            notPrime[i*Prime[j]] = true;
            if (i % Prime[j] == 0) break;
        }
    }
}

struct Matrix {
    int g[MAXP][MAXP];

    Matrix() { memset(g, 0, sizeof g); }

    Matrix operator * (const Matrix& rhs) const {
        Matrix ret;
        for (int i = 0; i < p; ++i)
            for (int k = 0; k < p; ++k)
                for (int j = 0; j < p; ++j)
                    ret.g[i][j] = (ret.g[i][j] + 1LL * g[i][k] * rhs.g[k][j] % MOD) % MOD;
        return ret;
    }

    inline void init() { for (int i = 0; i < p; ++i) g[i][i] = 1; }
};

Matrix fpow(Matrix base, int b) {
    Matrix ret; ret.init();
    while (b > 0) {
        if (b & 1) ret = ret * base;
        base = base * base, b >>= 1;
    }
    return ret;
}

int Part1() { // 所有数部分
    Matrix base;
    for (int i = 1; i <= m; ++i) ++f[i % p];
    for (int i = 1; i <= m; ++i) ++base.g[(-i % p + p) % p][0];
    for (int j = 1; j < p; ++j)
        for (int i = 0; i < p; ++i) base.g[i][j] = base.g[(i - 1 + p) % p][j-1];
    int ret = 0;
    Matrix Ans = fpow(base, n - 1);
    for (int i = 0; i < p; ++i)
        ret = (ret + 1LL * f[i] * Ans.g[i][0] % MOD) % MOD;
    return ret;
}

int Part2() { // 合数部分
    Matrix base;
    for (int i = 1; i <= m; ++i) if (notPrime[i]) ++g[i % p];
    for (int i = 1; i <= m; ++i) if (notPrime[i]) ++base.g[(-i % p + p) % p][0];
    for (int j = 1; j < p; ++j)
        for (int i = 0; i < p; ++i) base.g[i][j] = base.g[(i - 1 + p) % p][j-1];
    int ret = 0;
    Matrix Ans = fpow(base, n - 1);
    for (int i = 0; i < p; ++i)
        ret = (ret + 1LL * g[i] * Ans.g[i][0] % MOD) % MOD;
    return ret;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("input.in", "r", stdin);
#endif
    scanf("%d%d%d", &n, &m, &p);
    EulerSieve();
    printf("%d\n", (Part1() - Part2() + MOD) % MOD);
    return 0;
}

「SDOI2017」新生舞会

题目链接

https://loj.ac/problem/2003

解题思路

看到式子之后容易联想到 0/1 分数规划. 所以直接考虑二分答案.

假设当前答案 $\ge x$, 也就是

$C = \frac{ \sum a_i' }{ \sum b_i' } \ge x$

进一步转化,得

$\sum (a_i' - b_i' x) \ge 0$

然后就可以用二分图最大匹配解决了. 具体地说, 二分图内两两边权 $w_{i, j} = a_{i, j} - b_{i, j} x$, 如果最大权匹配 $\ge 0$ 则当前二分到的 $x$ 合法.

二分图最大权匹配部分采用 KM 算法, 时间复杂度 $O(n^3 \log w)$. 其中 $w$ 为权值上界.

突然有一种工程界很喜欢 KM 算法的错觉.

代码实现

听说写费用流有点卡常? 我这一个写 KM 的哪知道啊 = =

以及二分的时候不要玩弄 eps, 因为二分边界写挂调了好久.png

> 点击显示 / 隐藏

// LOJ #2003
// DeP
#include <cmath>
#include <cctype>
#include <cstdio>
#include <algorithm>
using namespace std;

namespace IO {
    const int MAXSIZE = 1 << 18 | 1;
    char buf[MAXSIZE], *p1, *p2;

    inline int Gc() {
        return p1 == p2 &&
            (p2 = (p1 = buf) + fread(buf, 1, MAXSIZE, stdin), p1 == p2)? EOF: *p1++;
    }
    template<typename T> inline void read(T& x) {
        x = 0; int f = 0, ch = Gc();
        while (!isdigit(ch)) f |= ch == '-', ch = Gc();
        while (isdigit(ch)) x = x * 10 + ch - '0', ch = Gc();
        if (f) x = -x;
    }
}
using IO::read;

const int MAXN = 105;
const double EPS = 1e-7, INFD = 1e18;

inline int dcmp(const double& p) {
    return (fabs(p) < EPS)? 0: (p < 0? -1: 1);
}

int n;
int A[MAXN][MAXN], B[MAXN][MAXN];
double W[MAXN][MAXN];

namespace KM {
    int S[MAXN], T[MAXN], Time;
    int left[MAXN];
    double Lx[MAXN], Ly[MAXN], slack[MAXN];

    bool dfs(int u) {
        S[u] = Time;
        for (int v = 1; v <= n; ++v) if (T[v] != Time) {
            double d = Lx[u] + Ly[v] - W[u][v];
            if (dcmp(d) == 0) {
                T[v] = Time;
                if (left[v] == -1 || dfs(left[v])) return left[v] = u, true;
            } else slack[v] = min(d, slack[v]);
        }
        return false;
    }

    inline void update() {
        double a = INFD;
        for (int j = 1; j <= n; ++j)
            if (T[j] != Time) a = min(a, slack[j]);
        for (int i = 1; i <= n; ++i) {
            if (S[i] == Time) Lx[i] -= a;
            if (T[i] == Time) Ly[i] += a; else slack[i] -= a;
            // 某篇讲 KM 算法的 CSDN 博客这里更新 slack 的时候写挂了 = =
        }
    }

    inline double KM() {
        for (int i = 1; i <= n; ++i)
            left[i] = -1, Lx[i] = *max_element(W[i]+1, W[i]+1+n), Ly[i] = 0;
        for (int u = 1; u <= n; ++u) {
            for (int i = 1; i <= n; ++i) slack[i] = INFD;
            while (true) {
                ++Time;
                if (!dfs(u)) update(); else break;
            }
        }
        double ret = 0;
        for (int u = 1; u <= n; ++u) if (~left[u]) ret += W[left[u]][u];
        return ret;
    }
}

inline bool check(const double& x) {
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j) W[i][j] = A[i][j] - x * B[i][j];
    return dcmp(KM::KM()) >= 0;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("input.in", "r", stdin);
#endif
    // input
    read(n);
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j) read(A[i][j]);
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j) read(B[i][j]);
    // solve
    double L = 0.0, R = 1e6;
    while (dcmp(R - L) > 0) {
        double Mid = (L + R) / 2.0;
        if (check(Mid)) L = Mid; else R = Mid;
    }
    // output
    printf("%.6lf\n", L);
    return 0;
}

「SDOI2017」硬币游戏

字符串题啊, 学提高的时候用哈希, 学省选的时候用 SAM 就够了.

不知道从哪里听来这句屁话 = =

题目链接

https://loj.ac/problem/2004

解题思路

~~还是 CQzhangyu 讲得好.~~

代码实现

0.5 ^ 300 = 4.909093465297727e-91

2.0 ^ 300 = 2.037035976334486e+90

你叫我用哪个…

> 点击显示 / 隐藏

// LOJ #2004
// DeP
#include <cmath>
#include <queue>
#include <cstdio>
#include <cassert>
#include <algorithm>
using namespace std;

const int MAXN = 305;
const double EPS = 1e-9;

inline int dcmp(const double& p) {
    return (fabs(p) < EPS)? 0: (p < 0? -1: 1);
}

int n, m;
char S[MAXN][MAXN];
double pow2[MAXN], A[MAXN][MAXN];

namespace AC {
    const int MAXN = ::MAXN * ::MAXN;
    int ch[2][MAXN], fail[MAXN], nidx;
    int len[MAXN], pos[MAXN];

    inline void init() { nidx = 1; }
    inline int val(const char& c) { return c == 'H'; }

    inline void insert(char* S, const int& idx) {
        int u = 1;
        for (int i = 1; i <= m; ++i) {
            int c = val(S[i]);
            if (!ch[c][u]) ch[c][u] = ++nidx;
            u = ch[c][u];
        }
        pos[idx] = u;
    }

    void getFail() {
        queue<int> Q;
        fail[1] = 1;
        for (int c = 0; c < 2; ++c) {
            int& v = ch[c][1];
            if (v) Q.push(v), fail[v] = 1; else v = 1;
        }
        while (!Q.empty()) {
            int u = Q.front(); Q.pop();
            for (int c = 0; c < 2; ++c) {
                int& v = ch[c][u];
                if (v) Q.push(v), fail[v] = ch[c][fail[u]];
                else v = ch[c][fail[u]];
            }
        }
    }

    inline void build() {
        pow2[0] = 1.0;
        for (int i = 1; i <= m; ++i) pow2[i] = pow2[i-1] * 0.5;
        for (int i = 1; i <= n; ++i) {
            for (int u = 1, j = 1; j <= m; ++j)
                u = ch[val(S[i][j])][u], len[u] = j;
            A[i][n + 1] = -pow2[m];
            for (int j = 1; j <= n; ++j)
                for (int u = pos[j]; u > 1; u = fail[u])
                    if (len[u]) A[i][j] += pow2[m - len[u]];
            for (int u = 1, j = 1; j <= m; ++j)
                u = ch[val(S[i][j])][u], len[u] = 0;
        }
        for (int i = 1; i <= n; ++i) A[n + 1][i] = 1;
        A[n + 1][n + 2] = 1;
    }
}

inline void Gauss() {
    for (int i = 1; i <= n + 1; ++i) {
        int r = i;
        for (int j = i+1; j <= n + 1; ++j)
            if (fabs(A[j][i]) > fabs(A[r][i])) r = j;
        // assert(dcmp(A[r][i]) != 0);
        if (r != i) swap(A[r], A[i]);
        for (int j = 1; j <= n + 1; ++j) if (i != j) {
            double d = A[j][i] / A[i][i];
            for (int k = i; k <= n + 2; ++k) A[j][k] -= d * A[i][k];
        }
    }
    for (int i = 1; i <= n + 1; ++i)
        if (dcmp(A[i][i]) != 0) A[i][n + 2] /= A[i][i];
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("input.in", "r", stdin);
#endif
    AC::init();
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= n; ++i)
        scanf("%s", S[i] + 1), AC::insert(S[i], i);
    AC::getFail(), AC::build();
    Gauss();
    for (int i = 1; i <= n; ++i) printf("%.7lf\n", A[i][n + 2]);
    return 0;
}

「SDOI2017」相关分析

之前听 lxl 讲洛谷网课的时候看他吐槽过 SDOI “那个无聊的东西”, 大概是指的这个吧 (

~~的确挺无聊的 = =~~

题目链接

https://loj.ac/problem/2005

解题思路

读过高中必修课本的 (忘了是哪本了 = =) 都知道, 询问的 $\hat a$ 可以展开, 当然不知道这个式子也是可推的…

$\hat a = \frac{\sum\limits_{i = L} ^ R (x_i - \bar{x})(y_i - \bar{y})}{\sum\limits_{i = L} ^ R (x_i - \bar{x}) ^ 2} = \frac{ \sum\limits_{i=L}^R x_i y_i - n \bar x \bar y }{\sum\limits_{i=L}^R x_i^2 - n \bar x^2}$

对于询问, 直接用线段树维护 $\sum {x_i}, \sum {y_i}, \sum x_i y_i, \sum x_i ^ 2$ 就可以回答了.

区间加… 把这几个东西加起来再展开就好了.

$\sum_{i=L}^R (x_i + S) = \sum_{i=L}^R x_i + (R-L+1) \cdot S$ $\sum_{i=L}^R (y_i + T) = \sum_{i=L}^R y_i + (R-L+1) \cdot T$ $\sum_{i=L}^R (x_i + S) (y_i + T) = \sum_{i=L}^R x_i y_i + T \sum_{i=L}^R x_i + S \sum_{i=L}^R y_i + (R-L+1) \cdot ST$ $\sum_{i=L}^R (x_i + S) ^ 2 = \sum_{i=L}^R x_i^2 + 2S \sum_{i=L}^R x_i + (R-L+1) \cdot S^2$

区间赋值… 注意所赋值具有特殊性, 也就是第 $i$ 个位置改成 $i$, 然后再跑一遍区间加的操作就好了.

~~看错条件多推了一倍的式子, 人没了~~

代码实现

线段树的细节无非是传标记的时候容易挂… 然后就喜闻乐见地挂了.

> 点击显示 / 隐藏

// LOJ #2005
// DeP
#include <cmath>
#include <cctype>
#include <cstdio>
using namespace std;

namespace IO {
    const int MAXSIZE = 1 << 18 | 1;
    char buf[MAXSIZE], *p1, *p2;

    inline int Gc() {
        return p1 == p2 &&
            (p2 = (p1 = buf) + fread(buf, 1, MAXSIZE, stdin), p1 == p2)? EOF: *p1++;
    }
    template<typename T> inline void read(T& x) {
        x = 0; int f = 0, ch = Gc();
        while (!isdigit(ch)) f |= ch == '-', ch = Gc();
        while (isdigit(ch)) x = x * 10 + ch - '0', ch = Gc();
        if (f) x = -x;
    }
}
using IO::read;

const int MAXN = 1e5 + 5;
const double EPS = 1e-9;

inline int dcmp(const double& p) {
    return (fabs(p) < EPS)? 0: (p < 0? -1: 1);
}

inline double s0(int L, int R) { return R - L + 1; }

inline double s1(int n) { return n / 2.0 * (n + 1); }
inline double s1(int L, int R) { return s1(R) - s1(L - 1); }

inline double s2(int n) { return n / 6.0 * (n + 1) * (2 * n + 1); }
inline double s2(int L, int R) { return s2(R) - s2(L - 1); }

int n, m;
int X[MAXN], Y[MAXN];

namespace SGT {
#define lc (nd<<1)
#define rc (nd<<1|1)
    struct Node {
        double x, x2, y, xy;
        Node() { x = y = xy = 0; }
        Node(double _x, double _x2, double _y, double _xy): x(_x), x2(_x2), y(_y), xy(_xy) { }
        Node operator + (const Node& rhs) const {
            return Node(x + rhs.x, x2 + rhs.x2, y + rhs.y, xy + rhs.xy);
        }
        inline double calc(const int& L, const int& R) {
            return (xy - x / s0(L, R) * y) / (x2 - x / s0(L, R) * x);
        }
    } dat[MAXN << 2];
    double tagAddX[MAXN << 2], tagAddY[MAXN << 2]; bool tagSet[MAXN << 2];

    inline void maintain(const int& nd) { dat[nd] = dat[lc] + dat[rc]; }

    inline void pushAdd(const int& nd, const int& L, const int& R, const double& S, const double& T) {
        dat[nd].x2 += 2.0 * S * dat[nd].x + s0(L, R) * S * S;
        dat[nd].xy += T * dat[nd].x + S * dat[nd].y + s0(L, R) * S * T;
        dat[nd].x += s0(L, R) * S;
        dat[nd].y += s0(L, R) * T;
        tagAddX[nd] += S, tagAddY[nd] += T;
    }

    inline void pushSet(const int& nd, const int& L, const int& R) {
        dat[nd].y = dat[nd].x = s1(L, R);
        dat[nd].x2 = dat[nd].xy = s2(L, R);
        tagSet[nd] = true, tagAddX[nd] = tagAddY[nd] = 0;
    }

    inline void pushdown(const int& nd, const int& L, const int& R) {
        int Mid = (L + R) / 2;
        if (tagSet[nd]) {
            pushSet(lc, L, Mid), pushSet(rc, Mid+1, R);
            tagSet[nd] = false;
        }
        if (dcmp(tagAddX[nd]) || dcmp(tagAddY[nd])) {
            pushAdd(lc, L, Mid, tagAddX[nd], tagAddY[nd]);
            pushAdd(rc, Mid+1, R, tagAddX[nd], tagAddY[nd]);
            tagAddX[nd] = tagAddY[nd] = 0;
        }
    }

    void build(int nd, int L, int R) {
        if (L == R) return void( dat[nd] = Node(X[L], 1.0 * X[L] * X[L], Y[L], 1.0 * X[L] * Y[L]) );
        int Mid = (L + R) / 2;
        build(lc, L, Mid), build(rc, Mid+1, R);
        maintain(nd);
    }

    void Add(int nd, int L, int R, const int& opL, const int& opR, const double& S, const double& T) {
        if (opL <= L && R <= opR) return pushAdd(nd, L, R, S, T);
        int Mid = (L + R) / 2;
        pushdown(nd, L, R);
        if (opL <= Mid) Add(lc, L, Mid, opL, opR, S, T);
        if (opR > Mid) Add(rc, Mid+1, R, opL, opR, S, T);
        maintain(nd);
    }

    void Mdy(int nd, int L, int R, const int& opL, const int& opR, const double& S, const double& T) {
        if (opL <= L && R <= opR) return pushSet(nd, L, R), pushAdd(nd, L, R, S, T);
        int Mid = (L + R) / 2;
        pushdown(nd, L, R);
        if (opL <= Mid) Mdy(lc, L, Mid, opL, opR, S, T);
        if (opR > Mid) Mdy(rc, Mid+1, R, opL, opR, S, T);
        maintain(nd);
    }

    Node Qry(int nd, int L, int R, const int& opL, const int& opR) {
        if (opL <= L && R <= opR) return dat[nd];
        pushdown(nd, L, R);
        int Mid = (L + R) / 2;
        if (opR <= Mid) return Qry(lc, L, Mid, opL, opR);
        if (opL > Mid) return Qry(rc, Mid+1, R, opL, opR);
        return Qry(lc, L, Mid, opL, opR) + Qry(rc, Mid+1, R, opL, opR);
    }
#undef lc
#undef rc
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("input.in", "r", stdin);
#endif
    // input
    read(n), read(m);
    for (int i = 1; i <= n; ++i) read(X[i]);
    for (int i = 1; i <= n; ++i) read(Y[i]);
    // solve
    SGT::build(1, 1, n);
    // output
    while (m--) {
        static int opt, L, R, S, T;
        static SGT::Node Ans;
        read(opt), read(L), read(R);
        switch (opt) {
            case 1: printf("%.8lf\n", SGT::Qry(1, 1, n, L, R).calc(L, R)); break;
            case 2: read(S), read(T), SGT::Add(1, 1, n, L, R, S, T); break;
            case 3: read(S), read(T), SGT::Mdy(1, 1, n, L, R, S, T); break;
            default: fprintf(stderr, "ERR\n");
        }
    }
    return 0;
}

DepletedPrism's Blog

SDOI 2017 Round 1 大赏

「SDOI2017」数字表格

题意简述

题目链接

代码实现

「SDOI2017」树点涂色

题目链接

解题思路

代码实现

「SDOI2017」序列计数

题目链接

解题思路

代码实现

「SDOI2017」新生舞会

题目链接

解题思路

代码实现

「SDOI2017」硬币游戏

题目链接

解题思路

代码实现

「SDOI2017」相关分析

题目链接

解题思路

代码实现