莫反好题,不仅仅是莫反,还有很多思维含量。
由于推式子过程太过于漫长了,所以我仅仅讲下大概。
题目是给你一个长度为 $n$ 的数组,请求出 $\sum\limits_{i=1}^n\sum\limits_{j=i+1}^n \operatorname{lcm}(A_i, A_j)$
莫反通常是对于值域考虑,直接推是不可行的,所以开一个桶 $b$,$b_k$ 存储 $A_i=k$ 的个数。
对值域分类的话就只能把忽略大小关系的 $i$ $j$ 对算出来后再容斥得到答案,即:
$\sum\limits_{i=1}^V\sum\limits_{j=1}^V \operatorname{lcm}(i,j)\cdot b_i \cdot b_j$。
多说一句,我刚开始没忽略大小关系,是 $\sum\limits_{i=1}^V\sum\limits_{j=i+1}^V \operatorname{lcm}(i,j)\cdot b_i \cdot b_j$,但这个实际上是错的,我当时没注意,推到最后发现没办法推了才想到这个办法的。。。
推式子,大概就是先把 $\operatorname{lcm}$ 变成相乘除以最大公约数之后枚举最大公约数,然后再同时除以最大公约数,然后再用莫比乌斯反演之后再枚举。
此时有一个整除的条件,枚举的时候枚举倍数就行了,最后把所有 $\sum$ 内的能拿出来的都拿出来,用乘法分配律拿。
推完后是:$\sum\limits_{d=1}^V d\cdot \sum\limits_{x=1}^{\lfloor\frac{V}{d}\rfloor}\cdot \mu(x)\cdot x^2\cdot\sum\limits_{i=1}^{\lfloor\frac{V}{dx}\rfloor}i\cdot b_{i\cdot d\cdot x}\cdot \sum\limits_{j=1}^{\lfloor\frac{V}{dx}\rfloor}j\cdot b_{j\cdot d\cdot x}$。
此时令 $f(x)=\sum\limits_{i=1}^{\lfloor\frac{V}{x}\rfloor}i\cdot b_{i\cdot x}$,然后就可以化为 $\sum\limits_{d=1}^V d\cdot \sum\limits_{x=1}^{\lfloor\frac{V}{d}\rfloor} \mu(x)\cdot x^2\cdot f^2(dx)$。
对 $f$ 进行预处理,然后调和级数 $n\log n$ 就可求得,注意最后要除以二,用逆元,还有减去 $\sum_{i=1}^n a_i$,因为我们多算了要容斥掉。
代码:
#include <bits/stdc++.h> #define int long long #define mp make_pair #define pi pair <double, double> #define For(i, a, b) for (int i = (a); i <= (b); i ++) #define foR(i, a, b) for (int i = (a); i >= (b); i --) using namespace std; int n, cnt1, cnt2; struct Node {int x, y, id;}a[105]; bool cmp (Node n1, Node n2) {return n1.x < n2.x || n1.x == n2.x && n1.y < n2.y;} pi up[205], dn[205]; int id1[205], id2[205]; double ans[205]; const double inf = 1e18, PI = 3.141592653; double k (pi a, pi b) { if (a.first == b.first) return inf; return (a.second - b.second) * 1.0 / (a.first - b.first); } double sq (int x) {return x * x;} double sqdis (pi a, pi b) {return sq (a.first - b.first) + sq (a.second - b.second);} double dis (pi a, pi b) {return sqrt (sqdis (a, b) );} void solve () { cin >> n; For (i, 1, n) { cin >> a[i].x >> a[i].y; a[i].id = i; } if (n == 1) { cout << 1; return; } if (n == 2) { cout << 0.5 << '\n' << 0.5 << '\n'; return; } sort (a + 1, a + n + 1, cmp); For (i, 1, n) { while (cnt1 >= 2) { if (k (up[cnt1], up[cnt1 - 1]) <= k (up[cnt1], mp (a[i].x, a[i].y) ) ) -- cnt1; else break; } up[++ cnt1] = mp (a[i].x, a[i].y); id1[cnt1] = a[i].id; while (cnt2 >= 2) { if (k (dn[cnt2], dn[cnt2 - 1]) >= k (dn[cnt2], mp (a[i].x, a[i].y) ) ) -- cnt2; else break; } dn[++ cnt2] = mp (a[i].x, a[i].y); id2[cnt2] = a[i].id; } foR (i, cnt2 - 1, 1) { up[++ cnt1] = dn[i]; id1[cnt1] = id2[i]; } up[0] = dn[2]; id1[0] = id2[2]; For (i, 1, cnt1 - 1) { ans[id1[i] ] = (PI - acos ( (sqdis (up[i], up[i - 1]) + sqdis (up[i], up[i + 1]) - sqdis (up[i - 1], up[i + 1]) ) / 2.0 / dis (up[i], up[i - 1]) / dis (up[i], up[i + 1]) ) ) / 2.0 / PI; if (up[i - 1] == up[i + 1]) ans[id1[i] ] = 0.5; } For (i, 1, n) cout << fixed << setprecision (20) << ans[i] << '\n'; } signed main () { int _ = 1; // cin >> _; while (_ --) { solve (); cout << '\n'; } return 0; }