Cow Camp (USACO Gold 2022 February)

官方解答中给的是从前向后递推的思路，这里给一个从后往前的：

考虑Bessie的最后一次尝试，这时得分就是二项分布，所以期望得分是 $T/2$ （sample case不考虑，所以开始把T减1，最后输出答案时再加1）。
然后我们考虑倒数第二次，比较容易想出来，如果Bessie这一次得分小于最后一次的期望得分 $T/2$ ，那么她就就应该继续尝试，如果大于这个，那她就应该停止。这样她的总的期望得分最高。
通过这个分析，我们尝试建立递推式。设 $p(x)$ = 二项分布值小于 $x$ 的概率， $q(x)$ = 二项分布值大于 $x$ 的部分的期望。则第 $i$ 次之后总期望为：

E_i = (1-p(E_{i+1})) \cdot q(E_{i+1}) + p(E_{i+1}) \cdot E_{i+1}

初值 $E_n = T/2$ ，最后就是要求 $E_1$ .

直接这样递推计算就可以得到一半的分。后面TLE原因是因为 $K$ 很大，所以需要一次跳多格，这里我们观察到中间 $p(x)$ 和 $q(x)$ 会有不变的情况，跳变就是在 $E$ 整数部分变化的时候。简单推下公式发现 $q-E$ 每次迭代都变成原来的 $p$ 倍，所以就是等比数列关系，二分查找，或者直接解方程就可以找到下一个点（见官方解答或者下面的代码）。

另一个要解决的问题是二项概率的求解， $1000 \choose k$ 是很大值，怎么放是个问题。浴谷上有一个解答是用long double，是可以的。这里给出一个不同的办法，就是使用杨辉三角型变型来求。

代码中还给了用log来计算的方法，但误差有点大，过不了测试点，undef掉了，也作为参考。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int,int> pi;
 
// 两种二项概率的求法：log法，以及pascal's triangle
// #define USELOG
 
int T, K;   // T < 1e3, K <= 1e9
double P[1005][1005];
double p[1005], cp[1005];       // prob and cumulated prob
double ex[1005];                // ex[x]: prob * x for >= x
 
#ifdef USELOG
double logfact[1005];
void fact() {
    for (int i = 1; i <= T; i++)
        logfact[i] = logfact[i-1] + log((double)i);
}
#endif
 
int main() {
    scanf("%d %d", &T, &K);
    T--;        // remove sample question
 
#ifdef USELOG
    fact();     // 用log法求应该也是对的，但通不过USACO的cases，应该是因为误差积累
    for (int i = 0; i <= T; i++) 
        p[i] = exp(logfact[T]-logfact[i]-logfact[T-i] - T*log(2));
#endif
    P[0][0] = 1.0;
    for (int i = 1; i <= T; i++) {  // Pascal's triangle for binomial probabilities
        P[i][0] = P[i-1][0] / 2.0;
        for (int j = 1; j <= T; j++)
            P[i][j] = (P[i-1][j-1] + P[i-1][j]) / 2.0;
        // printf("p[%d] = %lf\n", i, p[i]);
    }
#endif
 
    for (int i = 0; i <= T; i++) {
        p[i] = P[T][i];
        cp[i] = i ? cp[i-1] + p[i] : p[i];
    }
    for (int i = T; i >= 0; i--)
        ex[i] = i == T ? p[i] * i : ex[i+1] + p[i] * i;
 
    double e = (double)T/2;
    for (int j = 1; j < K; ) {
        double ee = cp[(int)e] * e + ex[(int)e+1];
        if ((int)ee != (int)e) {
            e = ee;
            j++;
        } else {
            // jump multiple points
            double P = cp[(int)e], Q = ex[(int)e+1];
            // e[x] = Q/(1-P) - (Q/(1-P) - e) * P^x
            double e2 = floor(e) + 1;
            int step = 1e8;
            if (Q/(1-P)-e2 > 0) 
                step = log((Q/(1-P)-e2) / (Q/(1-P)-e)) / log(P);
            double _e2 = Q/(1-P) - (Q/(1-P)-e)*pow(P, step);
            if (_e2 < e2 || step == 0)
                step++;
            if (j + step > K || step > 1e9)
                step = K-j;
            e2 = Q/(1-P) - (Q/(1-P)-e)*pow(P, step);
            // printf("From %d to %d (%d steps): %lf -> %lf\n", j, j+step, step, e, e2);
            e = e2;
            j += step;
        }
    }
    printf("%lf", e+1);
    return 0;
}