cf 1067D Computer Game

题意：有n个任务，每个任务有一次升级机会，升级之前完成任务可以得到$a_i$收益，升级之后完成可以得到$b_i$收益，每个任务可以尝试完成多次，每次完成了以后就可以得到相应的收益和升级机会。可以进行$t$次操作，求最优策略下期望收益。
数据范围：$1\le n\le 10^5$，$1\le t\le 10^{10}$，$a_i<b_i\le 10^8$,$0<p_i<1$
https://codeforces.com/contest/1067/problem/D

解法：这个题比较明显的性质就是后面这个升级，只要得到了，就直接升级$b\times p$最大的那个，一直做一定期望最大。然后我们考虑前面，开始我想的是一定会只用一个，那我就直接枚举，等比数列算一下就好了。交了一发，发现“WA on test 105”。105是最后一个点了，点开发现我的算法是错的……
为什么呢，因为考虑这样一个问题，如果你剩余轮数很多，那最优一定是做那个$p$大的；如果你剩余的轮数不多，但是你还是在不断的做$p$大的那个，就不如$a$大的优了。（居然有人用前面那个算法+特判水过了）
推翻重来，我们考虑一个朴素的做法——dp，方程长这样，设$dp[t]$表示还剩$t$个操作最优策略得到的期望收益，$M$表示升级后每轮期望收益，则
dpt+1=max{pi(tM+ai)+(1−pi)dpt}dp_{t+1} = max \lbrace p_i(tM +a_i) + (1-p_i)dp_t\rbracedpt+1​=max{pi​(tM+ai​)+(1−pi​)dpt​}
我们变一下式子
dpt+1=max{pitM+piai+dpt−pidpt}dp_{t+1} = max \lbrace p_itM + p_i a_i + dp_t -p_i dp_t \rbracedpt+1​=max{pi​tM+pi​ai​+dpt​−pi​dpt​}
dpt+1=max{pi(tM−dpt)+piai+dpt}dp_{t+1} = max \lbrace p_i(tM - dp_t) + p_ia_i + dp_t \rbracedpt+1​=max{pi​(tM−dpt​)+pi​ai​+dpt​}
这个时候$d{p}_t$与这些$p_i,a_i$无关，提出来
dpt+1=max{pi(tM−dpt)+piai}+dptdp_{t+1} = max \lbrace p_i(tM - dp_t) + p_ia_i \rbrace + dp_tdpt+1​=max{pi​(tM−dpt​)+pi​ai​}+dpt​
好，这个时候我们就变成每次求所有形如$y=p_{i}x+p_i{a}_i$直线，覆盖在$tM-d{p}_t$上最大的值是多少,明显这是一个凸包的问题（下凸）。但是轮数还是很多，无法得到明显优化。这个时候我们可以想一下前面的想法，那么思路肯定是先选$p_i$大的，再逐渐选$a_i$大的，$p$是斜率，我们就猜想它有单调性（单调向右），就可以证明一下，那我们要证明
$tM-d{p}_t\le (t+1)M-d{p}_{t+1}$
$tM-d{p}_t\le tM+M-d{p}_{t+1}$
$-d{p}_t\le M-d{p}_{t+1}$
$d{p}_t\ge d{p}_{t+1}-M$
$d{p}_{t+1}\le d{p}_t+M$
好，这个时候我们发现只有$d{p}_{t+1}$和$d{p}_t$了，从意义上入手，多一轮最多就多$M$的收益，那么不等式显然成立。
那我们就可以计算了，凸包的分界是确定的，然后每次考虑一段，利用矩阵乘法或二分等比数列求和，可以求出答案。
时间复杂度$O(n(log(n)+log(t))$
放上代码：

#include<iostream>
#include<cstring>
#include<cmath>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<cstdio>
#include<time.h>
#include<vector>
#include<algorithm>
using namespace std;
#define REP(i,x,y) for(int i=x;i<=y;i++)
#define rep(i,n) REP(i,1,n)
#define rep0(i,n) REP(i,0,n-1)
#define repG(i,x) for(int i=pos[x];~i;i=e[i].next)
#define ll long long
#define db double
const int N=1e5+7;
const ll INF=1e18+7;
const db eps=1e-10;
db p[N],vm[N];
ll a[N],b[N];
db maxn=0;
ll n,m=0,sz=1,t;
struct fac{db k,b;}f[N],dd[N],st[N];
bool cmp(fac x,fac y){return x.k<y.k;}
db meet(fac u,fac v){return (u.b-v.b)/(v.k-u.k);}
db ans[4][4],now[44][4][4],rp[4][4],cr3[4][4];
void cf(db cr0[4][4],db cr1[4][4],db cr2[4][4]){
	rep(i,3){
		rep(j,3){
			cr3[i][j]=0;
			rep(k,3)cr3[i][j]+=cr1[i][k]*cr2[k][j];
		}
	}
	rep(i,3)rep(j,3)cr0[i][j]=cr3[i][j];
}
db check(db q,ll t,db p[4][4]){return q*p[1][1]+t*p[1][2]+p[1][3];}
void work(fac v,ll &nw,db &dp,db R){
	if(nw==t)return;
	if(nw*maxn-dp>=R)return;
	now[0][1][1]=1-v.k;
	now[0][1][2]=v.k*maxn;
	now[0][1][3]=v.b;
	now[0][2][2]=now[0][2][3]=now[0][3][3]=1;
	now[0][2][1]=now[0][3][1]=now[0][3][2]=0;
	memset(ans,0,sizeof(ans));
	ll i=0,cs=0;
	for(;maxn*((1ll<<i)+nw)-check(dp,nw,now[i])<=R&&nw+(1ll<<i)<t;i++)cf(now[i+1],now[i],now[i]);
	rep(i,3)ans[i][i]=1;
	for(;~i;i--){
		rep(j,3)rep(k,3)rp[j][k]=ans[j][k];
		cf(ans,ans,now[i]);
		if((nw+cs+(1ll<<i))*maxn-check(dp,nw,ans)<=R&&nw+cs+(1ll<<i)<t)cs+=(1ll<<i);
		else rep(j,3)rep(k,3)ans[j][k]=rp[j][k];
	}
	cf(ans,ans,now[0]);
	db an=dp*ans[1][1]+nw*ans[1][2]+ans[1][3];
	dp=an;
	nw+=cs+1;
}
int main(){
	scanf("%I64d%I64d",&n,&t);
	rep(i,n)scanf("%I64d%I64d%lf",&a[i],&b[i],&p[i]);
	rep(i,n)f[i]=(fac){p[i],p[i]*a[i]};
	rep(i,n)maxn=max(maxn,b[i]*p[i]);
	sort(f+1,f+n+1,cmp);
	dd[++m]=f[1];
	REP(i,2,n){
		if(f[i].k>f[i-1].k+eps||f[i].k<f[i-1].k-eps)dd[++m]=f[i];
		else dd[m].b=max(dd[m].b,f[i].b);
	}
	rep(i,m)f[i]=dd[i];
	st[1]=f[1];
	rep(i,m){	
		while(sz>1){
			db r1=meet(st[sz-1],f[i]),r2=meet(st[sz],f[i]);
			if(r1<r2)break;
			sz--;
		}
		st[++sz]=f[i];
	}
	rep(i,sz-1)vm[i]=meet(st[i],st[i+1]);
	vm[sz]=INF;
	ll nw=0;
	db dp=0;
	rep(i,sz)work(st[i],nw,dp,vm[i]);
	printf("%.9lf\n",dp);
	return 0;
}