「解题报告」Sample

题目大意

从 \([1,n]\) 中随机选两个数，选中 \(i\) 的概率为 \(p_i\)，若两次选中的数相等，则权值为 \(0\)；否则，权值为两数之和。
求一种 \(p_i\) 序列满足 \(\sum p_i=1\)，使得最后得到的权值期望值最大。
\(n\le 10^6\)，保留六位小数，多测，\(T\le 10^5\)

思路

首先根据题意，可以写出期望为 \(E(p)=\sum_{i=1}^np_i(1-p_i)i=\sum_{i=1}^ni(p_i-p_i^2)\)。

---

首先先不考虑 \(p_i\ge 0\) 的条件，我们先看如何直接求这个函数的最值：

解法 1

考虑 拉格朗日乘数法：

拉格朗日乘数法

如果有多元函数 \(f(x_i)\)，要求在限制条件 \(\varphi(x_i)=0\) 下的极值，那么我们做拉格朗日函数 \(F(x_i,\lambda)=f(x_i)+\lambda \varphi(x_i)\)，令它对于 \(x_i,\lambda\) 的偏导数为 \(0\)，列出方程组，即：

\[\left\{\begin{aligned} F'_{x_1}&=f'_{x_1}(x_i)+\lambda \varphi'_{x_1}(x_i)=0\\ F'_{x_2}&=f'_{x_2}(x_i)+\lambda \varphi'_{x_2}(x_i)=0\\ \ \ \vdots\\ F'_{x_n}&=f'_{x_1}(x_i)+\lambda \varphi'_{x_n}(x_i)=0\\ F'_{\lambda}&=\varphi(x_i)=0\\ \end{aligned}\right.\]

那么解出的 \((x_i)\) 就是其极值点。

令 \(\varphi(p_i)=1-\sum_{i=1}^np_i\)，则

\[\left\{\begin{aligned} F'_{p_i}&=f'_{p_i}(p_i)+\lambda \varphi'_{p_i}(p_i)\\ &=(i(p_i-p_i^2))'+\lambda (-p_i')\\ &=i-2ip_i-\lambda=0\\ F'_{\lambda}&=\varphi(x_i)=0\\ \end{aligned}\right.\]

即：

\[\left\{\begin{aligned} i-2ip_i&=\lambda\\ \sum_{i=1}^np_i&=1\\ \end{aligned}\right.\]

解得：

\[\left\{\begin{aligned} p_i&=\frac{i-\lambda}{2i}\\ \lambda&=\frac{n-2}{H_n}\\ \end{aligned}\right.\]

解法 2

我们可以用另一个角度去考虑：设 \(f(p_i)=\sum_{i=1}^ng(p_i),g(p_i)=i(p_i^2-p_i)\)，那么如果 \(\exists i,j,g'(p_i)\ne g'(p_j)\)，那么令 \(p_i\) 和 \(p_j\) 进行一些调整，会使答案更优。那么当答案最大时，一定满足 \(\forall i,j,g'(p_i) = g'(p_j) = \lambda\)，实际上是和解法 1 等价的。

---

这样会发现，解出来的一些 \(p_i<0\)。这不好。所以我们考虑进行调整。

我们直接令一段前缀 \(p_i=0\)，然后再进行上面的计算，如果最后解出来的 \(p_i\) 都大于 \(0\)，那么就是合法的。

设 \(\forall i\in[1,b),p_i=0\)，那么对上面的式子进行一些调整，就可以得到：

\[\left\{\begin{aligned} i-2ip_i&=\lambda\\ \sum_{i=b}^np_i&=1\\ \end{aligned}\right.\]

注意下指标的更改。

解得：

\[\left\{\begin{aligned} p_i&=\frac{i-\lambda}{2i}\\ \lambda&=\frac{n-b-1}{H_n-H_{b-1}}\\ \end{aligned}\right.\]

考虑找到这个边界 \(b\)。

我们需要保证 \(p_i\ge0\)，也就是 \(\frac{i-\lambda}{2i}\ge0\)，即 \(i\ge \lambda\)，也就是要令：

\[b \ge \frac{n-b-1}{H_n-H_{b-1}} \]

可以通过 two-pointer 或者 二分 求出对于每一个 \(n\) 的边界 \(b\)。

令 \(h(x)=\frac{n-x-1}{H_n-H_{x-1}}-x\)，求导发现 \(h'(x)<0\)，所以可以二分。懒得导了

那么，可以求出答案：

\[\begin{aligned} E(p_i)&=\sum_{i=b}^np_i(1-p_i)i\\ &=\sum_{i=b}^ni(p_i-p_i^2)\\ &=\sum_{i=b}^ni\left(\left(\frac{1}{2}-\frac{n-b-1}{2i(H_n-H_{b-1})}\right)-\left(\frac{1}{2}-\frac{n-b-1}{2i(H_n-H_{b-1})}\right)^2\right)\\ &=\sum_{i=b}^ni\left(\frac{1}{2}-\frac{n-b-1}{2i(H_n-H_{b-1})}-\frac{1}{4}+\frac{n-b-1}{2i(H_n-H_{b-1})}-\frac{(n-b-1)^2}{4i^2(H_n-H_{b-1})^2}\right)\\ &=\sum_{i=b}^ni\left(\frac{1}{4}-\frac{(n-b-1)^2}{4i^2(H_n-H_{b-1})^2}\right)\\ &=\frac{(n+b)(n-b+1)}{8}-\frac{(n-b-1)^2}{4(H_n-H_{b-1})^2}\sum_{i=b}^n\frac{1}{i}\\ &=\frac{(n+b)(n-b+1)}{8}-\frac{(n-b-1)^2}{4(H_n-H_{b-1})}\\ \end{aligned}\]

代码

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1000005;
int T, n;
long double h[MAXN];
int b[MAXN];
int main() {
    freopen("sample.in", "r", stdin);
    freopen("sample.out", "w", stdout);
    b[0] = 1;
    for (int i = 1; i < MAXN; i++) {
        h[i] = h[i - 1] + 1.0l / i;
        b[i] = b[i - 1];
        while (b[i] < (i - b[i] - 1.0l) / (h[i] - h[b[i] - 1])) b[i]++;
    }
    scanf("%d", &T);
    while (T--) {
        scanf("%d", &n);
        printf("%.12Lf\n", 2.0l * (1.0l * (n + b[n]) * (n - b[n] + 1) / 8 - 1.0l * (n - b[n] - 1) * (n - b[n] - 1) / 4.0l / (h[n] - h[b[n] - 1])));
    }
    return 0;
}