UESTC2022暑假前(?)集训 数学与几何

数学与几何解题报告

其实还有一个求圆上的整点个数的题，写完忘保存丢掉了。

目录

• 数学与几何解题报告
• B-新月之舞 推式子
  • 题意
  • 解题思路
  • 代码
• C-回转数 几何
  • 题意
  • 解题思路
  • 代码
• E-兔儿爷 几何
  • 题意
  • 解题思路
  • 代码
• F-多重和式 反演
  • 题意
  • 解题思路
  • 代码
• G-数三角形 几何
  • 题意
  • 解题思路
  • 代码
• H-洛水神女 水题
  • 题意
  • 解题思路
  • 代码
• N-黑灰游戏 改 博弈
  • 题意
  • 解题思路
  • 代码
• Q-挑战新高考 反演 杜教筛
  • 题意
  • 解题思路
  • 代码
• S- For D. E. Knuth 拓展欧拉定理
  • 题意
  • 解题思路
  • 代码
• W-奥日与风袭废墟 线性筛
  • 题意
  • 解题思路
  • 代码
• X-多项式乘法 FFT
  • 题意
  • 解题思路
  • 代码

B-新月之舞 推式子

题意

求$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}i \bmod j$

$1\le n\le 10^{12}$

解题思路

$原式=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{i}i \bmod j+\sum\limits_{i=1}^{n}\sum\limits_{j=i+1}^{n}i \bmod j$

$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{i}(i - \left \lfloor \frac{i}{j} \right \rfloor j )+\sum\limits_{i=1}^{n}i(n-i)$

$=\sum\limits_{i=1}^{n}i^2-\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{i} \left \lfloor \frac{i}{j} \right \rfloor j +\sum\limits_{i=1}^{n}i(n-i)$

$=\sum\limits_{i=1}^{{n}ni-\sum\limits_{i=1}}\sum\limits_{j=1}^{i} \left \lfloor \frac{i}{j} \right \rfloor j $

$=\frac{n^{2(n+1)}{2}-\sum\limits_{i=1}}\sum\limits_{j=1}^{i} \left \lfloor \frac{i}{j} \right \rfloor j $

结论：$\sum\limits_{j=1}^{i}\left \lfloor \frac{i}{j} \right \rfloor j =\sum\limits_{j=1}^{i} \sigma(j)$

证明：考虑$[1,i]$中的每个数出现在了几个$\sigma(j)$中，有$\sum\limits_{j=1}^{i} \sigma(j)=\sum\limits_{j=1}^{i} \sum\limits_{d|j}d=\sum\limits_{d=1}^{i}\left \lfloor \frac{i}{d} \right \rfloor d=\sum\limits_{j=1}^{i}\left \lfloor \frac{i}{j} \right \rfloor j$ $\Box$

因此$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{i} \left \lfloor \frac{i}{j} \right \rfloor j= \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{i} \sigma(j)$

$= \sum\limits_{i=1}^{n}\sigma(i)(n+1-i)$

$= (n+1)\sum\limits_{i=1}^{n}\sigma(i)-\sum\limits_{i=1}^{n}i\sigma(i)$

根据前面的结论，

$(n+1)\sum\limits_{i=1}^{n}\sigma(i)=(n+1)\sum\limits_{d=1}^{n}d\left \lfloor \frac{n}{d} \right \rfloor$

只需求出$\sum\limits_{i=1}^{n}i\sigma(i)$

类似地，我们考虑$id$对答案的贡献，可以得到

$\sum\limits_{i=1}^{n}i\sigma(i)=\sum\limits_{i=1}^{n}\sum\limits_{d|i}id=\sum\limits_{d=1}^{n}(d+2d+\dots+\left \lfloor \frac{n}{d} \right \rfloor d)d$

=$\sum\limits_{d=1}^{n}\frac{\left \lfloor \frac{n}{d} \right \rfloor(\left \lfloor \frac{n}{d} \right \rfloor+1)}{2}d^2$

综上所述，$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}i \bmod j=\frac{n^2(n+1)}{2}-(n+1)\sum\limits_{d=1}^{n}d\left \lfloor \frac{n}{d} \right \rfloor+\sum\limits_{d=1}^{n}\frac{\left \lfloor \frac{n}{d} \right \rfloor(\left \lfloor \frac{n}{d} \right \rfloor+1)}{2}d^2$

借助数论分块，可以$O(\sqrt n)$求出答案。

数论分块是这么一个东西:

出现$\left \lfloor \frac{n}{i} \right \rfloor$时，我们可以把$n$个数相加，转化成$O(\sqrt n)$个区间和相加，每一个区间里$\left \lfloor \frac{n}{i} \right \rfloor$是相等的。

例如最简单的$H(n)=\sum\limits_{i=1}^{n}\left \lfloor \frac{n}{i} \right \rfloor$

int H(int n) {
    int res = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res += (r - l + 1) * (n / l);
    }
    return res;
}

答案会达到$O(n^3)$，需要使用int128。

代码

//
// Created by vv123 on 2022/7/4.
//
#include <bits/stdc++.h>
#define int __int128
using namespace std;

inline int read() {
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9'){ if(ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9'){ x = x * 10 + ch - '0'; ch = getchar(); }
    return x * f;
}
inline void print(int x) {
    if(x < 0) { putchar('-'); x = -x; }
    if(x > 9) print(x / 10);
    putchar(x % 10 + '0');
}

int n, ans;


inline int A() {
    int res = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        res += (l + r) * (r - l + 1) * (n / l) / 2;
    }
    return res;
}

inline int sumsquare(int x) {
    return x * (x + 1) * (x + x + 1) / 6;
}

inline int B() {
    int res = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        int t = n / l;
        res += (sumsquare(r) - sumsquare(l - 1)) * t * (t + 1) / 2;
    }
    return res;
}

void solve() {
    n = read();
    ans =  n * n * (n + 1) / 2 - (n + 1) * A() + B();
    print(ans);putchar('\n');
}

signed main() {
    int T = read();
    while (T--) {
        solve();
    }
    return 0;
}

C-回转数 几何

题意

给出一段由$n$条有向线段组成的有向封闭折线，可能有交叉和重叠，并给出$m$个点，对每个点求折线围绕该点的回转数。

$3\le n\le 5000, 1\le m \le5000$

解题思路

• to-left测试：判断一个点P是否在一条有向直线AB左侧，只需要看$\overrightarrow{AB} \times\overrightarrow{AP}$是否大于0.

• 回转数：面内闭合曲线逆时针绕过该点的总次数。性质：当回转数为0，点在曲线外部。

• 朴素做法：是计算相邻两边夹角的和，需要使用浮点数。

• 无精度误差的解法：向右作一条射线，与朝上的边相交回转数+1，朝下的边相交回转数-1。

• 朝上的边与射线相交，当且仅当检测点纵坐标介于uv之间，且检测点在uv左侧。另一种情况同理。

• 细节：如果射线恰好经过顶点，不妨定义顶点在射线上方。（只统计$miny \le y_p\lt maxy$产生的答案）

代码

//
// Created by vv123 on 2022/7/18.
//
#include <bits/stdc++.h>
#define int long long
#define pii pair<int,int>
using namespace std;

const int N = 2e4 + 10;
struct point {
    int x, y;
    point operator-(const point &p) const { return {x - p.x, y - p.y}; }
    int operator*(const point &p) const { return x * p.x + y * p.y; }
    int operator^(const point &p) const { return x * p.y - p.x * y; }
    int quad() const {
        if (y > 0 || (y == 0 && x >= 0)) return 1;
        return 2;
    }//下平面定义为[-180,0),上平面定义为[0,180)
    int toleft(const point &p) { //p在这个向量的逆时针方向，则叉积为正
        int t = (*this) ^ p;
        if (t > 0) return 1;
        if (t == 0) return 0;
        if (t < 0) return -1;
    }
};

bool cmp(const point &a,const point &b) {
    if (a.quad() != b.quad()) return a.quad() < b.quad();
    else return (a ^ b) > 0;
}

struct line {
    point p, v;
    int toleft(const point &a) {
        return v.toleft(a - p);
    }
};

struct polygon {
    vector<point> vec;
    int nxt(int i) { return (i + 1) % vec.size(); }
    pii winding(const point &p) {
        int res = 0;
        for (int i = 0; i < vec.size(); i++) {
            point u = vec[i], v = vec[nxt(i)];
            if (((p - u) ^ (p - v)) == 0 && (p - u) * (p - v) <= 0) return {1, 0};
            if (u.y == v.y) continue;
            line uv = {u, v - u};
            if (u.y < v.y && uv.toleft(p) <= 0) continue;
            if (u.y > v.y && uv.toleft(p) >= 0) continue;
            if (u.y < p.y && v.y >= p.y) res++;
            if (u.y >= p.y && v.y < p.y) res--;
        }
        return {0, res};
    };
} poly;

void solve() {
    int n, m;
    cin >> n;
    for (int i = 1; i <= n; i++) {
        int x, y;
        cin >> x >> y;
        poly.vec.push_back({x, y});
    }
    cin >> m;
    while (m--) {
        int x, y;
        cin >> x >> y;
        point p = {x, y};
        pii res = poly.winding(p);
        if (res.first == 1) puts("EDGE");
        else cout << res.second << "\n";
    }
}
signed main() {
    int T = 1;
    //cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

E-兔儿爷 几何

题意

有三瓶药水，每瓶药水两种溶质的浓度分别为$a_i,b_i$。$n$次询问问能否配制出两种溶质浓度分别为$x,y$的药水。

$1\le n\le 10^5$

解题思路

容易发现"把两种药水按一定比例混合"，结果相当于"两个向量按一定比例相加"（系数非负，且和为1）。

我们可以将问题转化成，给定三个向量，是否存在一组和为1的非负系数，使得三个向量的线性组合为$(x,y)$。

或者说，$(x,y)$能否被表达为三个向量的凸组合。

根据凸包的定义，只需要检验$(x,y)$是否在三瓶药水构成的凸包（三角形）内即可。

这里直接采用了C题的回转数法。

代码

//
// Created by vv123 on 2022/7/18.
//
#include <bits/stdc++.h>
#define int long long
#define pii pair<int,int>
using namespace std;

const int N = 2e4 + 10;
struct point {
    int x, y;
    point operator-(const point &p) const { return {x - p.x, y - p.y}; }
    int operator*(const point &p) const { return x * p.x + y * p.y; }
    int operator^(const point &p) const { return x * p.y - p.x * y; }
    int quad() const {
        if (y > 0 || (y == 0 && x >= 0)) return 1;
        return 2;
    }//下平面定义为[-180,0),上平面定义为[0,180)
    int toleft(const point &p) { //p在这个向量的逆时针方向，则叉积为正
        int t = (*this) ^ p;
        if (t > 0) return 1;
        if (t == 0) return 0;
        if (t < 0) return -1;
    }
};

bool cmp(const point &a,const point &b) {
    if (a.quad() != b.quad()) return a.quad() < b.quad();
    else return (a ^ b) > 0;
}

struct line {
    point p, v;
    int toleft(const point &a) {
        return v.toleft(a - p);
    }
};

struct polygon {
    vector<point> vec;
    int nxt(int i) { return (i + 1) % vec.size(); }
    pii winding(const point &p) {
        int res = 0;
        for (int i = 0; i < vec.size(); i++) {
            point u = vec[i], v = vec[nxt(i)];
            if (((p - u) ^ (p - v)) == 0 && (p - u) * (p - v) <= 0) return {1, 0};
            if (u.y == v.y) continue;
            line uv = {u, v - u};
            if (u.y < v.y && uv.toleft(p) <= 0) continue;
            if (u.y > v.y && uv.toleft(p) >= 0) continue;
            if (u.y < p.y && v.y >= p.y) res++;
            if (u.y >= p.y && v.y < p.y) res--;
        }
        return {0, res};
    };
} poly;

void solve() {
    int n, m;
    //cin >> n;
    n = 3;
    for (int i = 1; i <= n; i++) {
        int x, y;
        cin >> x >> y;
        poly.vec.push_back({x, y});
    }
    cin >> m;
    while (m--) {
        int x, y;
        cin >> x >> y;
        point p = {x, y};
        pii res = poly.winding(p);
        if (res.first == 1 || res.second != 0) puts("YES");
        else puts("NO");
    }
}
signed main() {
    int T = 1;
    //cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

F-多重和式 反演

题意

给定$n, m, a(1\le n,m,a\le10^7)$, 求$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} lcm (i,j)[\gcd(i,j)\le a]$

解题思路

• 莫比乌斯函数：$\mu(n)=\left\{\begin{array}{cc} 1 & n=1 \\ (-1)^{|P|} & n=\prod_{i} p_{i} \\ 0 & \exists p_{i}, p_{i}^{2} \mid n \end{array}\right.$
• 反演结论：$\sum\limits_{d|n}\mu(d)=[n=1]$
• 莫比乌斯变换：
  • 形式1 $F(n)=\sum\limits_{d|n}f(d)\Rightarrow f(n)=\sum\limits_{d|n}\mu(d)F(\frac{n}{d})$
  • 形式2 $F(n)=\sum\limits_{n|d}f(d)\Rightarrow f(n)=\sum\limits_{n|d}\mu(\frac{d}{n})F(d)$
• 基础例题：$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} [\gcd(i,j)= 1]=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{d|gcd(i,j)}\mu(d)=\sum\limits_{d=1}^{n}\mu(d)\left \lfloor \frac{n}{d} \right \rfloor^2$(考虑枚举d)

考虑将lcm转化为gcd，有

$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} lcm (i,j)[\gcd(i,j)\le a]=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} \frac{ij}{\gcd(i,j)}[\gcd(i,j)\le a]$

$= \sum\limits_{d=1}^{\min(n, m, a)}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{ij}{d}[\gcd(i,j)=d]$

注意到两个数的gcd等于d，当且仅当它们除以d后互质。

$= \sum\limits_{d=1}^{\min(n, m, a)}\sum\limits_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum\limits_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}\frac{id\cdot jd}{d}[\gcd(i,j)=1]$

$= \sum\limits_{d=1}^{\min(n, m, a)}d \sum\limits_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum\limits_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}ij[\gcd(i,j)=1]$

令$g(n,m)= \sum\limits_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum\limits_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}ij[\gcd(i,j)=1]$,则原式$= \sum\limits_{d=1}^{\min(n, m, a)}d\cdot g(\left \lfloor \frac{n}{d} \right \rfloor,\left \lfloor \frac{m}{d} \right \rfloor)$

考虑如何求出$g(n, m)$.

令$f(d)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[\gcd(i,j)=d]$（n,m为已确定的常数）

那么$f(1)=g(n,m)$

令$F(x)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[x|\gcd(i,j)]$,那么$F(x)=\sum\limits_{x|d}f(d)$

根据形式2，$f(1)=\sum\limits_{d}^{min(n,m)}\mu(d)F(d)$

下面考虑如何快速求出$F(d)$。根据前面类似的手法，可以得到

$F(d)=\sum\limits_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum\limits_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}id\cdot jd[\gcd(i,j)=1]=d^2\frac{\left \lfloor \frac{n}{d} \right \rfloor(\left \lfloor \frac{n}{d} \right \rfloor+1)}{2}\frac{\left \lfloor \frac{m}{d} \right \rfloor(\left \lfloor \frac{m}{d} \right \rfloor+1)}{2}$

从而$g(n,m)=f(1)=\sum\limits_{d}^{min(n,m)}\mu(d)d^2\frac{\left \lfloor \frac{n}{d} \right \rfloor(\left \lfloor \frac{n}{d} \right \rfloor+1)}{2}\frac{\left \lfloor \frac{m}{d} \right \rfloor(\left \lfloor \frac{m}{d} \right \rfloor+1)}{2}$

代入原式$=\sum\limits_{d=1}^{\min(n, m, a)}d\cdot g(\left \lfloor \frac{n}{d} \right \rfloor,\left \lfloor \frac{m}{d} \right \rfloor)$求解即可

使用了两层整除分块，时间复杂度为$O(n)$

这里第一次接触二维的整除分块。一般只需要将$r =\left \lfloor \frac{n}{\left \lfloor \frac{n}{l} \right \rfloor} \right \rfloor $改成$r =\min(\left \lfloor \frac{n}{\left \lfloor \frac{n}{l} \right \rfloor} \right \rfloor,\left \lfloor \frac{m}{\left \lfloor \frac{m}{l} \right \rfloor} \right \rfloor)$即可。但要注意到这时r可能大于a，因此本题应取$r =\min(\left \lfloor \frac{n}{\left \lfloor \frac{n}{l} \right \rfloor} \right \rfloor,\left \lfloor \frac{m}{\left \lfloor \frac{m}{l} \right \rfloor} \right \rfloor， a)$

代码

//
// Created by vv123 on 2022/7/19.
//
#include <bits/stdc++.h>
#define int long long
using namespace std;


int n, m, a, ans;
const int mod = 1e9 + 7;
const int N = 1e7 + 10;

int pow(int a, int b, int p) {
    int res = 1;
    for (; b; b >>= 1) {
        if (b & 1) res = res * a % mod;
        a = a * a % mod;
    }
    return res;
}

int vis[N], prime[N], mu[N];

int s[N], sum[N];

void Shai(int n = N - 1) {
    vis[1] = mu[1] = 1;
    int cnt = 0;
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) {
            prime[++cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * prime[j] <= n; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0; //p1的指数大于1，莫比乌斯函数为0
                break;
            } else {
                mu[i * prime[j]] = -mu[i];
            }
        }
    }
    s[0] = sum[0] = 0;
    for (int i = 1; i <= n; i++) {
        s[i] = (s[i - 1] + i * i % mod * mu[i] % mod) % mod;
        sum[i] = (sum[i - 1] + i) % mod;
    }
}



inline int g(int n, int m) {
    int res = 0;
    for (int l = 1, r; l <= min(n, m); l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        int t = (n / l) * (n / l + 1) % mod * (m / l)  %  mod * (m / l + 1)  % mod * pow(4, mod - 2, mod) % mod;
        res = (res + (s[r] - s[l - 1] + mod) % mod * t % mod) % mod;
    }
    return res;
}

void solve() {
    Shai();
    cin >> n >> m >> a;
    int ans = 0;
    for (int l = 1, r; l <= min(min(n, m), a); l = r + 1) {
        r = min(a, min(n / (n / l), m / (m / l)));
        int t = g(n / l, m / l);
        //cout << l << " " << r << " " << (sum[r] - sum[l - 1] + mod) % mod * t % mod << "\n";
        ans = (ans + (sum[r] - sum[l - 1] + mod) % mod * t % mod) % mod;
    }
    cout << ans << "\n";
}

signed main() {
    int T = 1;
    //cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

G-数三角形 几何

题意

在二维平面上有 $n$个点，请问以这些点为顶点可以形成多少个不同的锐角三角形。（保证坐标两两不同）

$3\le n\le2000$

解题思路

我们将锐角定义为$[0,\frac{\pi}{2})$,非锐角定义为$[\frac{\pi}{2},\pi]$。那么任选三个点，只有两种情况：

1. 有一个非锐角和两个锐角。包括直角三角形、钝角三角形和三点共线。
2. 有三个锐角，为锐角三角形。

因此设锐角、非锐角的个数分别为$R,D$,显然锐角三角形的个数为$\frac{R-2D}{3}$

又因为总共有$n \binom{n-1}{2}$个角，只需统计出非锐角的个数。考虑枚举每一个点作为顶点，将其他点进行极角排序。

这里将极角按照$[-\pi \to\pi)$排序。考虑将第$i$条边作为始边，指针$l$指向逆时针方向的一个大于等于$\frac{\pi}2{}$（第一个不满足(v[i] ^ v[l]) >= 0 && v[i] * v[l] > 0)的位置，指针$r$指向逆时针方的第一个大于$\frac{3\pi}{2}$（第一个满足(v[i] ^ v[l]) < 0 && v[i] * v[l] > 0）的位置，则以第$i$条边作为始边构成的非锐角个数为$r - l$；逆时针枚举每一条边，指针$l, r$是单调的，可以$O(n)$统计出非锐角个数。

时间复杂度的瓶颈来自$n$次极角排序，为$O(n^2\log n)$。

代码

//
// Created by vv123 on 2022/7/13.
//
#include <bits/stdc++.h>
#define int long long
using namespace std;

const int N = 2e4 + 10;
struct point {
    int x, y;
    point operator-(const point &p) const { return {x - p.x, y - p.y}; }
    int operator*(const point &p) const { return x * p.x + y * p.y; }
    int operator^(const point &p) const { return x * p.y - p.x * y; }
    int quad() const {
        if (y > 0 || (y == 0 && x >= 0)) return 1;
        return 2;
    }//下平面定义为[-180,0),上平面定义为[0,180)
} p[N], v[N];

bool cmp(const point &a,const point &b) {
    if (a.quad() != b.quad()) return a.quad() < b.quad();
    else return (a ^ b) > 0;
}

int n, R, D;

void work(int id) {
    point O = p[id];
    int m = 0;
    for (int i = 1; i <= n; i++) {
        if (i == id) continue;
        v[++m] = p[i] - O;
    }
    sort(v + 1, v + 1 + m, cmp);
    //for (int i = 1; i <= m; i++) cout << v[i].x << " " << v[i].y << "\n";

    int l = 1, r = 1;
    for (int i = 1; i <= m; i++) {
        while (l <= m && (v[i] ^ v[l]) >= 0 && v[i] * v[l] > 0) l++;
        while (r <= m && ((v[i] ^ v[r]) >= 0 || v[i] * v[r] <= 0)) r++;
        //cout << i << " " << l << " " << r << "\n";
        D += r - l;
    }
}

void solve() {
    D = 0;
    cin >> n;
    for (int i = 1; i <= n; i++) {
        cin >> p[i].x >> p[i].y;
    }
    for (int i = 1; i <= n; i++) {
        work(i);
    }
    //cout << D << endl;
    R = n * (n - 1) * (n - 2) / 2 - D;
    cout << (R - D * 2) / 3 << "\n";
}
signed main() {
    int T = 1;
    //cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

H-洛水神女 水题

题意

有$n$只猫，每只猫在生日之后有$P_i$的概率死亡，求他们的期望寿命之和 对$10^9 + 7$取模.

$n \le 5\times 10^5$

解题思路

设一只猫每次生日后的死亡率为$p$,则它的寿命期望为

$\sum \limits_{i=1}^{{\infty}i(1-p)}p =-p\sum \limits_{i=1}^{{\infty}[(1-p)}i]'=-p[\sum \limits_{i=1}^{{\infty}(1-p)}i]'=-p(\frac{1}{p}-1)'=\frac{1}{p} $

因此只需求出$\sum\frac{1}{P_i}=\sum a_i^{-1}b_i \pmod{10^9+7}$

根据费马小定理，$a^{-1} \equiv a^{p-2}\pmod p, 这里p = 10^9 + 7$

代码

//
// Created by vv123 on 2022/7/12.
//
#include <bits/stdc++.h>
#define int long long
using namespace std;

const int mod = 1e9 + 7;

typedef long long ll;

pair<signed, signed> nextPair(signed &x) {
    signed a, b;
    x ^= x << 13; x ^= x >> 17; x ^= x << 5; x %= 10000000; a = x;
    x ^= x << 13; x ^= x >> 17; x ^= x << 5; x %= 10000000; b = x;
    b = max(b, 1); a = max(a % b, 1);
    //cout << a << " " << b << "\n";
    return make_pair(a, b);
}

int pow(int a, int b, int p) {
    int res = 1;
    for (; b; b >>= 1) {
        if (b & 1) res = res * a % p;
        a = a * a % p;
    }
    return res;
}



void solve() {
    int n; signed x;
    cin >> n >> x;
    ll ans = 0;
    for (int i = 1; i <= n; i++) {
        pair<signed, signed> t = nextPair(x);
        ans = (ans + 1ll * t.second % mod * 1ll * pow(t.first, mod - 2, mod)) % mod;
    }
    cout << ans << "\n";
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    int T = 1;
    //cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

N-黑灰游戏 改 博弈

题意

有$n$个石子堆，每堆有$a_i$个石子。

奇数轮由Kate指定一堆石子，Emilico从中拿走至少一个石子。

偶数轮由Emilico指定一堆石子，Kate从中拿走至少一个石子。

拿走最后一个石子的一方获胜。问Emilico是否有必胜策略。

解题思路

我们将石子堆分为两类：

• A类：只有一个石子

• B类：多于一个石子

考虑最简单的情况：如果所有石子堆都只有一个石子(都是A类堆)，显然当堆数为奇数时Emilico获胜。

考虑只有一个B类堆。

如果上一种情况的胜者指定了B类堆，上一种情况的败者一定会将它变成A类堆，改变胜负。

如果上一种情况的败者指定了B类堆，上一种情况的胜者一定会把它拿完，依然获胜。

不管哪种情况，都是指定了B类堆的一方失败。

因此，双方一定不会主动指定B类堆，不得不指定B类堆的一方必败。

容易发现：A类堆个数为偶数时，B类堆一定会被Kate指定，从而Emilico获胜。否则Kate获胜。

下面使用数学归纳法证明：存在K(K ≥ 1)个B类堆时，结论都是"A类堆个数为偶数数时Emilico获胜"

首先K=1的情况已被证明。

假设当B类堆个数为K时，结论是"A类堆个数为偶数时Emilico获胜"，则当B类堆个数为K+1时：

考虑哪一方指定了第一个B类堆。

如果B类堆个数为K时的胜者指定了第一个B类堆，B类堆个数为K时的败者一定会将它变成A类堆，使胜负与B类堆个数为K的情况相反，得以取胜。

如果B类堆个数为K时的败者指定了第一个B类堆，B类堆个数为K时的胜者者一定会将它取完，使胜负与B类堆个数为K的情况相同，得以取胜。

因此不管哪种情况，都是不得不指定第一个B类堆的一方失败，由Kate还是Emilico被迫指定第一个B类堆，只与A类堆的数量奇偶性有关。结论仍为"A类堆个数为偶数时Emilico获胜"。因此结论对任意正整数K都成立。

综上所述，只有奇数个A类堆，或者存在B类堆且A类堆个数为为偶数时,Emilico获胜。

代码

//
// Created by vv123 on 2022/7/20.
//
#include <bits/stdc++.h>
#define int long long
#define pii pair<int,int>
#define seed 131
using namespace std;

const int mod = 1e9 + 7;
const int N = 2e6 + 10;
const int inf = 0x3f3f3f3f;


void solve() {
    int n, t, cntA = 0, cntB = 0;
    cin >> n;
    for (int i = 1; i <= n; i++) {
        cin >> t;
        if (t == 1) cntA++;
        else cntB++;
    }
    if ((cntB == 0) ^ (cntA % 2 == 0)) puts("Win");
    else puts("Lose");
}


signed main() {
    int T = 1;
    cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

Q-挑战新高考 反演 杜教筛

题意

从L至R中的R-L+1个整数中随机取2个不同的数。求这2个数互质的概率。

$1\le T\le10, 1\le L, R\le 10^9$

解题思路

我们在课上学过

• $\mu(n)=\left\{\begin{array}{cc} 1 & n=1 \\ (-1)^{|P|} & n=\prod_{i} p_{i} \\ 0 & \exists p_{i}, p_{i}^{2} \mid n \end{array}\right.$

• $\sum\limits_{d|n}\mu(d)=[n=1]$

• $\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n} [\gcd(i,j)= 1]=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{d|gcd(i,j)}\mu(d)=\sum\limits_{d=1}^{n}\mu(d)\left \lfloor \frac{n}{d} \right \rfloor^2$

同样地，我们考虑$d$的倍数在$[L,R]$中的出现次数，有$\sum\limits_{i=L}^{R}\sum\limits_{j=L}^{R} [\gcd(i,j)= 1]=\sum\limits_{d=1}^{R} \mu(d)\left(\left\lfloor\frac{R}{d}\right\rfloor-\left\lfloor\frac{L-1}{d}\right\rfloor\right)^{2}$

（ppt上的结论有误）

注意到互质的$(i, j)$除$(1,1)$外都可以互换，因此如果$L=1$，上述和式的结果是符合要求的排列数+1,否则恰好为所有可能的排列个数。

因此此题答案为$\frac{\sum\limits_{i=L}^{R}\sum\limits_{j=L}^{R} [\gcd(i,j)= 1]-[L=1]}{(R-L+1)(R-L)}=\frac{\sum\limits_{d=1}^{R} \mu(d)\left(\left\lfloor\frac{R}{d}\right\rfloor-\left\lfloor\frac{L-1}{d}\right\rfloor\right)^{2}-[L=1]}{(R-L+1)(R-L)}$

用数论分块进行计算的过程中需要用到莫比乌斯函数的$10^9$前缀和，可以使用杜教筛。

杜教筛被用来处理数论函数的前缀和问题。对于求解一个前缀和，杜教筛可以在低于线性时间的复杂度内求解。

对于数论函数 $f$，要求我们计算 $S(n)=\sum_{i=1}^{n}f(i)$.

我们想办法构造一个 $S(n)$ 关于 $S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)$ 的递推式

对于任意一个数论函数 $g$，必满足

$$\begin{aligned} \sum_{i=1}^{n}\sum_{d \mid i}g(d)f\left(\frac{i}{d}\right)&=\sum_{i=1}^{n}g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\\ \iff \sum_{i=1}^{n}(f\ast g)(i)&=\sum_{i=1}^{n}g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \end{aligned}$$

略证：

$g(d)f\left(\frac{i}{d}\right)$ 就是对所有 $i\leq n$ 的做贡献，因此变换枚举顺序，枚举 $d,\frac{i}{d}$（分别对应新的 $i,j$）

$$\begin{split} &\sum_{i=1}^n\sum_{d \mid i}g(d)f\left(\frac{i}{d}\right)\\ =&\sum_{i=1}^n\sum_{j=1}^{\left\lfloor\frac{n}{i}\right\rfloor}g(i)f(j)\\ =&\sum_{i=1}^ng(i)\sum_{j=1}^{\left\lfloor\frac{n}{i}\right\rfloor}f(j)\\ =&\sum_{i=1}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \end{split}$$

那么可以得到递推式

$$g(1)S(n)=\sum_{i=1}^n(f\ast g)(i)-\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)$$

那么假如我们可以快速对 $\sum_{i=1}^n(f \ast g)(i)$ 求和，并用数论分块求解 $\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)$ 就可以在较短时间内求得 $g(1)S(n)$.

由狄利克雷卷积，我们知道：

$\because \epsilon =\mu \ast 1$（$\epsilon(n)=~[n=1]$）

$\therefore \epsilon (n)=\sum_{d \mid n} \mu(d)$

$S_1(n)=\sum_{i=1}^n \epsilon (i)-\sum_{i=2}^n S_1(\lfloor \frac n i \rfloor)$

$= 1-\sum_{i=2}^n S_1(\lfloor \frac n i \rfloor)$

观察到 $\lfloor \frac n i \rfloor$ 最多只有 $O(\sqrt n)$ 种取值，我们就可以应用数论分块来计算每一项的值了。

直接计算的时间复杂度为 $O(n^{\frac 3 4})$。考虑先线性筛预处理出前 $n^{\frac 2 3}$ 项，剩余部分的时间复杂度为

$O(\int_{0}^{n^{\frac 1 3}} \sqrt{\frac{n}{x}} ~ dx)=O(n^{\frac 2 3})$

对于较大的值，可以用 unordered_map 存下其对应的值，方便以后使用时直接使用之前计算的结果。

代码

//
// Created by vv123 on 2022/7/29.
//
#include <bits/stdc++.h>
#define int long long
using namespace std;


int n, m, a, ans;
const int N = 1e7 + 10;

int min(int a, int b, int c) {
    return min(min(a, b), c);
}
int pow(int a, int b, int p) {
    int res = 1;
    for (; b; b >>= 1) {
        if (b & 1) res = res * a % p;
        a = a * a % p;
    }
    return res;
}

int gcd(int a, int b) {
    return b ? gcd(b, a % b) : a;
}

int vis[N], prime[N], mu[N], sum[N];

void Shai(int n = N - 1) {
    vis[1] = mu[1] = 1;
    int cnt = 0;
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) {
            prime[++cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * prime[j] <= n; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0; 
                break;
            } else {
                mu[i * prime[j]] = -mu[i];
            }
        }
    }
    for (int i = 1; i <= n; i++) {
        sum[i] = sum[i - 1] + mu[i];
    }
}

unordered_map<int, int> sum_mu;
inline int get_mu(int x) {
    if (x < N) return sum[x];
    if (sum_mu[x]) return sum_mu[x];
    int res = 1;
    for (int l = 2, r; l <= x; l = r + 1) {
        r = x / (x / l);
        res -= (r - l + 1) * get_mu(x / l);
    }
    return sum_mu[x] = res;
}


int f(int L, int R) {
    int res = 0;
    for (int l = 1, r; l <= L; l = r + 1) {
        r = min(L / (L / l), R / (R / l));
        int t = ((R / l) - (L / l)) * ((R / l) - (L / l));
        res += (get_mu(r) - get_mu(l - 1)) * t;
    }
    for (int l = L + 1, r; l <= R; l = r + 1) {
        r = R / (R / l);
        int t = (R / l) * (R / l);
        res += (get_mu(r) - get_mu(l - 1)) * t;
    }
    return res;
}

int bf(int L, int R) {
    int res = 0;
    for (int i = L; i <= R; i++)
        for (int j = L; j <= R; j++)
            if (gcd(i, j) == 1) res++;
    return res;
}

void solve() {
    int L, R;
    cin >> L >> R;
    int a = f(L-1, R) - (L == 1), b = (R - L + 1) * (R - L), d = gcd(a, b);
    cout << a / d << "/" << b / d << "\n";
}

signed main() {
    Shai();
    int T = 1;
    cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

S- For D. E. Knuth 拓展欧拉定理

题意

求$a\uparrow \uparrow b\pmod p$

$1\le a,b,p\le 10^9$

解题思路

前置知识：拓展欧拉定理

$$a^b\equiv \begin{cases} a^{b\bmod\varphi(p)},\,&\gcd(a,\,p)=1\\ a^b,&\gcd(a,\,p)\ne1,\,b<\varphi(p)\\ a^{b\bmod\varphi(p)+\varphi(p)},&\gcd(a,\,p)\ne1,\,b\ge\varphi(p) \end{cases} \pmod p$$

由题意，$a\uparrow\uparrow1=a, a\uparrow\uparrow b=a^{a\uparrow\uparrow (b-1)}$

于是当$a>1$且$b$充分大时，有

$a\uparrow\uparrow b\pmod p=a^{a\uparrow\uparrow (b-1)}\pmod p$

$=a^{a\uparrow\uparrow (b-1)\bmod \varphi(p)+\varphi(p)}\pmod p$

$=a^{a^{a\uparrow\uparrow(b-2)}\bmod \varphi(p)+\varphi(p)}\pmod p$

$=a^{a^{a\uparrow\uparrow(b-2)\bmod \varphi(\varphi(p))+\varphi(\varphi(p))}\bmod \varphi(p)+\varphi(p)}\pmod p$

$=\dots$

注意到欧拉函数的多次嵌套很快会收敛到$\varphi(1)=1$,因此可以递归计算。

我们设$f(a,b,p)=(a\uparrow \uparrow b)\pmod p$

则有$f(a,b,p)\equiv\left\{\begin{array}{cc}1&a=1，p>1\\0 & p=1 \\a^{f(a,b-1,\varphi(p))+\varphi(p)} &a>1,p>1,b\gt5\\暴力计算&b\le5 \end{array}\right.$

......

我们发现$2\uparrow\uparrow5=2^{65536}$是一个很大的数字

这可以确保$a\uparrow\uparrow (b-1)\ge\varphi(p)$在$a>1,b\gt5$时一定成立.

并且由于$a^{\varphi(p)}\equiv1\pmod p,\gcd(a,p)=1$,因此拓展欧拉定理的第三条对$\gcd(a,p)=1$也适用

因此采取这种方法计算，结果一定是正确的

代码

//
// Created by vv123 on 2022/7/10.
//
#include <bits/stdc++.h>
#define int long long
using namespace std;

int a, b, p;

inline int phi(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 gcd(int a, int b) {
    return b ? gcd(b, a % b) : a;
}

int pow(int a, int b, int p) {
    if (gcd(a, p) == 1) b %= phi(p);
    if (gcd(a, p) != 1 && b >= phi(p)) b = b % phi(p) + phi(p);
    int res = 1;
    for (; b; b >>= 1) {
        if (b & 1) res = res * a % p;
        a = a * a % p;
    }
    return res;
}

int f(int a, int b, int p) {
    if (p == 1) return 0;
    else if (b > 5){
        int t = pow(a, f(a, b - 1, phi(p)) + phi(p), p);
        //printf("f(%d,%d,%d)=%d\n", a, b, p, t);
        return t;
    } else {
        int res = 1;
        for (int i = 1; i <= b; i++) {
            res = pow(a, res, p);
        }
        return res;
    }
}

void solve() {
    cin >> a >> b >> p;
    if (a == 1) {
        if (p == 1) puts("0");
        else cout << 1 << "\n";
        return;
    }
    cout << f(a, b, p)  << "\n";
}

signed main() {
    int T = 1;
    //cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

W-奥日与风袭废墟 线性筛

题意

求$\varphi ,\mu,\sigma_0, \sigma_1, (\mu \cdot \varphi)\ast \sigma_1$的前$500000$项

解题思路

以上积性函数在质数的取值通常可以由定义得到。

而对于合数，线性筛保证了每次由最小的质因子筛掉。我们设合数$n = p \times n'$，$p$是它的最小质因子。则有以下两种情况

1. $n' \bmod p \ne 0$,则$p$与$n'$互质，根据积性函数的性质，$f(n) = f(p)f(n')$

2. $n' \bmod p = 0$,则$p$也是$n'$的最小质因子。

   设$n=p_1^{k_1}p_2^{k_2}...$。其中$p_1$是它的最小质因子。我们求一个辅助函数$g(n)=f(p_1^{k_1})$

   那么$f(n) = \frac{f(n')}{g(n')}g(n)$

而$g(n)$是很好求的:

1. 当$n' \bmod p \ne 0$,$n$被唯一的最小质因子$p$筛去，$g(n)=f(p)$。

2. 当$n' \bmod p = 0$，只需在$g(n')$的基础上添加一个$p$的贡献即可。

代码

//
// Created by vv123 on 2022/7/8.
//
#include <bits/stdc++.h>
#define int long long
using namespace std;

const int N = 500000 + 10;
int n;
int vis[N], prime[N], phi[N], mu[N], c[N], d[N], g[N], f[N], F[N], G[N], gg[N];

void Shai() {
    vis[1] = phi[1] = mu[1] = d[1] = g[1] = f[1] = G[1] = gg[1] = F[1] = 1;
    int cnt = 0;
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) {
            prime[++cnt] = i;
            phi[i] = i - 1;
            mu[i] = -1;
            c[i] = 1; //质数的最小质因子指数为1
            d[i] = c[i] + 1;   //质数的约数个数为2。d[n] = PI(c_i + 1)
            g[i] = i + 1;//g[i]为i的最小质因子0~k次方和,约数和f[i]为所有k次质因子g[i]的积
            gg[i] = i;
            f[i] = g[i];
            G[i] = mu[1] * phi[1] * f[i] + mu[i] * phi[i] * f[1];
            F[i] = G[i];
        }
        for (int j = 1; j <= cnt && i * prime[j] <= n; j++) {
            //printf("i=%d prime[j]=%d\n", i, prime[j]);
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                //printf("%d mod %d == 0\n", i, prime[j]);
                phi[i * prime[j]] = prime[j] * phi[i];  //n'包含了n的所有质因子，与n的欧拉函数只差p1
                mu[i * prime[j]] = 0; //p1的指数大于1，莫比乌斯函数为0
                c[i * prime[j]] = c[i] + 1; //注意到i和i * prime[j]的最小质因子都是prime[j],后面都很显然了
                d[i * prime[j]] = d[i] / (c[i] + 1) * (c[i * prime[j]] + 1);
                g[i * prime[j]] = g[i] * prime[j] + 1;
                gg[i * prime[j]] = gg[i] * prime[j];
                f[i * prime[j]] = f[i] / g[i] * g[i * prime[j]];
                G[i * prime[j]] = G[i] + gg[i];         //
                F[i * prime[j]] = F[i] / G[i] * G[i * prime[j]];

                break;
            } else {
                //printf("%d mod %d == %d\n", i, prime[j], i % prime[j]);
                phi[i * prime[j]] = phi[prime[j]] * phi[i]; //最小质因子不整除其他部分，n'和p1互质
                mu[i * prime[j]] = mu[prime[j]] * mu[i];
                c[i * prime[j]] = 1;
                d[i * prime[j]] = d[prime[j]] * d[i];
                g[i * prime[j]] = 1 + prime[j];
                gg[i * prime[j]] = prime[j];
                f[i * prime[j]] = f[prime[j]] * f[i];
                G[i * prime[j]] = G[prime[j]];
                F[i * prime[j]] = F[prime[j]] * F[i];
            }
        }
    }
}
 
void solve() {
    cin >> n;
    Shai();
    for (int i = 1; i <= n; i++) printf("%d ", phi[i]); puts("");
    for (int i = 1; i <= n; i++) printf("%d ", mu[i]); puts("");
    for (int i = 1; i <= n; i++) printf("%d ", d[i]); puts("");
    for (int i = 1; i <= n; i++) printf("%d ", f[i]); puts("");
    //for (int i = 1; i <= n; i++) printf("%d ", gg[i]); puts("");
    //for (int i = 1; i <= n; i++) printf("%d ", G[i]); puts("");
    for (int i = 1; i <= n; i++) printf("%d ", F[i]); puts("");
}

signed main() {
    int T = 1;
    //cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

X-多项式乘法 FFT

题意

给出两个多项式的系数，求两个多项式乘积的系数

$0\le n,m\le10^5,系数是0到9的整数$

解题思路

FFT是这样一种东西：通过由分治实现的DFT，我们可以$O(n\log n)$将一个多项式从系数表示法转化成点值表示法。已知两个多项式的点值表示，则它们的乘积的点值表示就是点值的乘积。我们得到乘积的点值表示后，可以经过 IDFT（傅里叶反变换）的作用转换成系数形式 。

代码

//
// Created by vv123 on 2022/7/11.
//
#include <bits/stdc++.h>
#define int long long
#define pi M_PI
using namespace std;
typedef complex<double> cp;
const cp I(0, 1);
const int N = 5e6 + 10;

cp t[N], f[N], g[N], h[N];
int a[N], b[N];

void fft(cp *f, int n, int rev) {
    if (n == 1) return;
    for (int i = 0; i < n; i++) t[i] = f[i];
    for (int i = 0; i < n; i++) {
        if (i % 2) f[n / 2 + i / 2] = t[i];
        else f[i / 2] = t[i];
    }
    cp *g = f, *h = f + n / 2;
    fft(g, n / 2, rev); fft(h, n / 2, rev);
    cp cur(1, 0), step = exp(I * (2 * pi / n * rev));
    for (int k = 0; k < n / 2; k++) {
        t[k] = g[k] + cur * h[k];
        t[k + n / 2] = g[k] - cur * h[k];
        cur *= step;
    }
    for (int i = 0; i < n; i++) f[i] = t[i];
}



void solve() {
    int n, m;
    cin >> n >> m;
    for (int i = 0; i <= n; i++) {
        int x;
        cin >> x;
        f[i].real(x);
    }
    for (int i = 0; i <= m; i++) {
        int x;
        cin >> x;
        g[i].real(x);
    }
    int k = 1;
    while (k <= n + m) k <<= 1;
    fft(f, k, 1); fft(g, k, 1);
    for (int i = 0; i <= k; i++)
        h[i] = f[i] * g[i];
    fft(h, k, -1);
    for (int i = 0; i <= n + m; i++)
        cout << (int) ((h[i].real() / k) + 0.5)  << " ";
}



signed main() {
    int T = 1;
    //cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}