CF643C 题解

看见题解里面全是斜率优化，没有李超线段树，所以写一篇纪念。

不想用李超线段树的也别略过啊 QAQ 还是有题目分析的。

题目大意：

长度为 $n$ 的序列 $a$ 代表 $n$ 个关卡，划分成连续的 $m$ 段，记 $i$ 所属段左端点为 $l_i$ 。对于第 $i$ 个关卡用刚好一小时过关概率是 $\frac{a_i}{\sum_{j=l_i}^i a_j}$ 。如果没有过关要继续打这一关到过。求通关时间期望。 $n\leqslant 2\times 10^5,m\leqslant 50$ 。

原题的变量名设的和这里不一样的，在本题解内都以本题解的描述为主。

题目分析

Part1 期望

考虑一小时通过一关的概率是 $\dfrac{a_i}{\sum_{j=l_i}^i a_j}$ 。过一关需要的时间的期望。首先为了书写方便，我们用换元法，记 $x$ 表示分子， $y$ 表示分母。

如果你知道通过一关的时间期望怎么算，请直接跳转 Part2，如果只知道结论不知道推导，建议看一下 ~~因为里面也打了我博客的广告~~。

通过一关的时间的期望：

\begin{aligned} E &= \frac{x}{y}\cdot 1+\frac{y-x}{y}\cdot\frac{x}{y}\cdot 2 +(\frac{y-x}{y})^2\cdot \frac{x}{y}\cdot 3 + \cdots \\ &= \frac{x}{y}[1+ 2\cdot\frac{y-x}{y} + 3\cdot(\frac{y-x}{y})^2 + 4\cdot(\frac{y-x}{y})^3 + \cdots] \\ \end{aligned}

中括号里面的东西就是带系数的等比数列求和（ $\sum_{k=2}^\infty k\cdot p^{k-1}$ ）的形式。我的这篇博客给出了形如 $\sum_{k=0}^\infty p^k$ 和 $\sum_{k=0}^\infty k\cdot p^k$ 的求值，虽然不一样但是变通一下即可。可以得到：

\begin{aligned} \sum_{k=2}^\infty k\cdot p^{k-1} &= \sum_{k=1}^\infty p^k+\sum_{l=1}^\infty\sum_{k=l}^\infty p^k \\ &= \frac{p}{1-p} + \sum_{l=1}^\infty \frac{p^l}{1-p} \\ &= \frac{p}{1-p} + \frac{p}{(1-p)^2} \\ &= \frac{p(2-p)}{(1-p)^2} \end{aligned}

$p=\frac{y-x}{y}=1-\frac{x}{y}$ ，代入：

\begin{aligned} E &= \frac{x}{y}[1 + 2\cdot\frac{y-x}{y} + 3\cdot(\frac{y-x}{y})^2 + 4\cdot(\frac{y-x}{y})^3 + \cdots] \\ &= \frac{x}{y}[1 + \frac{(1-\frac{x}{y})(1+\frac{x}{y})}{(\frac{x}{y})^2}] \\ &= \frac{x}{y} + \frac{y}{x}[1-(\frac{x}{y})^2] \\ &= \frac{x}{y} + \frac{y}{x} - \frac{x}{y} \\ &= \frac{y}{x} \end{aligned}

到这里我们证得，通过一关的时间期望为 $\dfrac{\sum_{j=l_i}^i a_j}{a_i}$ 。

Part2 前缀和与暴力DP

会暴力 $\operatorname{O}(n^2m)$ DP 的可以直接跳转到 Part3。

显然可以前缀和，设 $s_i=\sum_{j=1}^i a_i$ ，通过一关的期望就为：（注：原题用的分段指闭区间，但是为了方便前缀和后文全部用左开右闭的区间表示） $\dfrac{s_i-s_{l_i}}{a_i}$ 。通过一整段关的时间期望是：

\begin{aligned} \sum_{j=l_i}^{r_i} \frac{s_j-s_{l_i}}{a_j} &= \sum_{j=l_i}^{r_i} (\frac{s_j}{a_j}-s_{l_i}\cdot\frac{1}{a_j}) \\ &= \sum_{j=l_i}^{r_i} \frac{s_j}{a_j} - s_{l_i}\sum_{j=l_i}^{r_i} \frac{1}{a_j} \end{aligned}

显然这两个都可以用前缀和，我们记录 $sa_i=\sum_{j=1}^i \frac{s_j}{a_j}$ ， $sb_i=\sum_{j=1}^i \frac{1}{a_j}$ 。

这样发现如果确定了一个分段的左右区间，整段的期望时间都可以 $\operatorname{O}(1)$ 地求出来。到这里不妨可以尝试概率 DP，设 $f_{i,j}$ 表示考虑前 $i$ 个关，且当前已经分出了 $j$ 段（显然按顺序数的第 $j$ 段右端点就是 $i$ ）的最优期望时间。转移时考虑当只分了 $j-1$ 段时，最后一段如何拆成两段即可。转移方程为：

f_{i,j}=\min_{k=1}^{i-1} \{ f_{k,j-1}+ sa_i - sa_k + sb_k\cdot s_k - sb_i\cdot s_k \}

把只和 $i$ 有关的提出来，为：

f_{i,j}=sa_i+\min_{k=1}^{i-1} \{ f_{k,j-1} - sa_k + sb_k\cdot s_k - sb_i\cdot s_k \}

暴力求解是 $\operatorname{O}(n^2m)$ 的。

Part3 李超线段树

我们发现这东西很好用和斜率相关的东西维护。~~我知道你们都用斜率优化吊打我了但是我还是要写李超~~。而最小值函数里面的一陀东西，只有一项和 $i$ 有关，这更加确定了这题可以不用任何骚操作就能李超优化。

我们考虑把 $-s_k$ 当作斜率， $f_{k,j-1} - sa_k + sb_k\cdot s_k$ 当作截距。原转移的最小值函数里面就是一个 $kx+b$ 的形式了。查询这个东西的最小值，不就是查询一条 $x=sb_i$ 的直线，和这些直线交点的最小值吗？

直接李超线段树维护上述的直线即可。但是如果你不幸地 WA 在了第 11 个数据点，明显是因为精度误差爆掉了。（CF 数据上前 10 个点里面还有一个 $n\approx 10^3$ 的，没有挂的话转移基本没问题了）换句话说我们不得不面对动态开点的李超线段树存在较大的精度误差的问题。

其实动态开点的李超线段树不过就是做到了空间复杂度和直线数量（一般就是 $n$ ）同级。考虑到如果改成普通的线段树，就会出现 $nm\cdot 4$ 的尴尬空间复杂度，稳稳地 MLE（为了方便值域的 $10^6$ 我就也写成 $n$ 了）。

一个显而易见的方法可以解决这类问题，考虑到 DP 的状态转移中，对于任意的 $f_{i,j}$ ，它只由 $f_{k,j-1}$ 的状态转移过来，可以进行滚动数组优化。李超线段树也跟着滚动一下，空间复杂度优化至 $n\cdot 8$ ，精度是不容易卡的，这就十分稳定能过了。

接下来考虑第二个问题：不动态开点的李超线段树处理精度问题，一般会需要提前把要查询的浮点数离散化，然后查询，比如这样：

f[i][now]=s1[i]+seg[now].ask(1,1,n,lower_bound(s2+1,s2+1+n,c[i])-s2);
// 这就是我交的 TLE on #52 的代码的转移

虽然找到离散化后对应的数是 $\log n$ ，李超线段树也是 $\log n$ ，但是考虑到他还是 STL，我写的李超线段树也人傻常数大，需要简单改一下：

for(register int i=1;i<=n;++i)
	pos[i]=lower_bound(s2+1,s2+1+n,c[i])-s2;
...... // 中间省略部分语句
f[i][now]=s1[i]+seg[now].ask(1,1,n,pos[i]);

这样就卡过去了，最慢的点差不多 2500ms，理论的时间复杂度是 $nm\log n$ 级别，相对来说还可以。

Part4 代码

#include<cstdio>
#include<algorithm>
#define int long long // 要开 long long！
using namespace std;
const int max_n=1000005,inf=1000000000000000LL;
inline int read(){
   int x=0;bool w=0;char c=getchar();
   while(c<'0' || c>'9') w|=c=='-',c=getchar();
   while(c>='0' && c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
   return w?-x:x;
}
inline int min(int a,int b){return a<b?a:b;}

int n,m,a[max_n],s[max_n];

double s1[max_n],s2[max_n];

struct LiChao{ // 不动态开点李超
   int p[max_n*4];
   #define ls(x) (x<<1)
   #define rs(x) (x<<1|1)

   double b[max_n],k[max_n];

   inline double slp(int j,int i){return b[i]+s2[j]*k[i];}
   // 计算 x=s2[j] 时直线的 y 值。不是算斜率

   inline void add(int x,int l,int r,int t){
   	if(l==r){
   		if(slp(l,t)<slp(l,p[x]))
   			p[x]=t;
   		return;
   	}
   	int mid=(l+r)>>1;
   	if(slp(mid,t)<slp(mid,p[x]))
   		t^=p[x]^=t^=p[x];
   
   	if(slp(l,t)<slp(l,p[x])) add(ls(x),l,mid,t);
   	else if(slp(r,t)<slp(r,p[x])) add(rs(x),mid+1,r,t);
   }
   inline double ask(int x,int l,int r,int t){
   	if(!x) return 0;
   	double res=slp(t,p[x]);
   	if(l==r) return res;
   	int mid=(l+r)>>1;
   	if(t<=mid) return min(res,ask(ls(x),l,mid,t));
   	return min(res,ask(rs(x),mid+1,r,t));
   }
   inline void clear(){ // 清空线段树，状压用
   	for(register int i=0;i<=n;++i)
   		k[i]=b[i]=p[i]=0;
   }
   #undef ls
   #undef rs
}seg[2]; // 滚动数组压在李超上

double f[max_n][2]; // 滚动 j 
double c[max_n];
int pos[max_n];

signed main(){
   n=read(),m=read();
   for(register int i=1;i<=n;++i){
   	a[i]=read();
   	s[i]=s[i-1]+a[i]; // s 即题解的 s
   	s1[i]=1.0*s[i]/a[i], // s1 即题解的 sa
   	s2[i]=1.0/a[i]; // s2 即题解的 sb
   }
   for(register int i=1;i<=n;++i){
   	s1[i]+=s1[i-1],
   	s2[i]+=s2[i-1];
   	c[i]=s2[i]; // c 数组用于离散化
   }
   sort(s2+1,s2+1+n);
   for(register int i=1;i<=n;++i)
   	f[i][0]=inf; // 初始化，分成 0 段是不合法的
   for(register int i=1;i<=n;++i)
   	pos[i]=lower_bound(s2+1,s2+1+n,c[i])-s2;
   for(register int j=1,now=1;j<=m;++j,now^=1){
   	seg[now].clear(); // 滚动掉过去的状态
   	for(register int i=1;i<=n;++i){
   		f[i][now]=s1[i]+seg[now].ask(1,1,n,pos[i]);

   		seg[now].k[i]=-1.0*s[i];
   		seg[now].b[i]=s2[i]*s[i]-s1[i]+f[i][now^1];
   		seg[now].add(1,1,n,i);
   	}
   }
   printf("%.10lf",f[n][m&1]);
   return 0;
}

upd

DP 转移式写错了，已经改好；修正错别字，删除代码冗杂注释。