CF1067D. Computer Game（斜率优化＋倍增＋矩阵乘法）

题目链接

https://codeforces.com/contest/1067/problem/D

题解

首先，如果我们获得了一次升级机会，我们一定希望升级 (b_i imes p_i) 最大的任务，并且之后只完成该任务，这样才能使得期望收益最大。换句话说，当我们完成成功了一次任务之后，决策就固定了。因此，我们实际需要考虑的是还未完成任何任务时的决策。

为了方便，我们记 (maxlimits_{1 leq i leq n}{b_ip_i}) 为 (m)。

我们设 (f_t) 表示还未完成任何任务且还剩 (t) 秒时使用最优决策得到的期望收益。如果我们选择在这一秒完成任务 (i)，那么期望收益即为：(p_i imes ((t - 1)m + a_i) + (1 - p_i)f_{t - 1})。因此，转移即为：

[egin{aligned}f_t &= maxlimits_{1 leq i leq n}{p_i imes ((t - 1)m + a_i) + (1 - p_i)f_{t - 1}} \ &= maxlimits_{1 leq i leq n}{p_i((t - 1)m - f_{t - 1}) + p_ia_i} + f_{t - 1}end{aligned} ]

直接用该转移式做 dp 的时间复杂度是 (O(tn)) 的。我们再将转移式变形：(-((t - 1)m - f_{t - 1})p_i + f_t - f_{t - 1} = p_ia_i)，该式子显然可以使用斜率优化，将 ((t - 1)m - f_{t - 1}) 看做斜率 (k)，那么所求即为斜率为 (k) 且经过点 ((-p_i, p_ia_i)) 的直线在 (y) 轴上的最大的截距。我们预处理出所有 (n) 个点 ((-p_i, p_ia_i)) 构成的上凸包后，每次在凸包上二分找最优点即可。这样，我们得到了一个时间复杂度为 (O(t log n)) 的较为优秀的做法。不过本题中的 (t leq 10^{10})，因此这样的做法依然无法通过。

不过幸运的是，在本题中，随着 (t) 的增加，斜率 (k)（即 ((t - 1)m - f_{t - 1})）是单调不下降的。

首先证明上述结论。对于任意的 (t > 0)，我们希望证明 (tm - f_t geq (t - 1)m - f_{t - 1})。将式子变形，即为 (f_{t} - f_{t - 1} leq m)。考虑该不等式的实际意义，显然 (f_{t - 1}) 与 (f_t) 的差不会超过 (m)，因为我们无法在 (1) 秒内得到大于 (m) 的收益。因此，结论是成立的。

该结论告诉我们了一个重要信息：凸包上的每个点所影响的 (f) 的下标 (t) 一定是一段连续的区间。由于凸包由不超过 (n) 个点构成，因此我们可以依次求出每一个点所对应的 (f) 的下标的范围。分析 (f) 的转移式子：(f_t = maxlimits_{1 leq i leq n}{p_i((t - 1)m - f_{t - 1}) + p_ia_i} + f_{t - 1})，当 (p_i, a_i) 一定时，该转移是可以使用矩阵乘法优化的，因此对于每一个点，我们可以预处理出在该点上的 (2^k) 次转移对应的矩阵，然后倍增求出该点对应的 (f) 下标的范围边界。这样，本题就在 (O(nomega^3 log t)) 的时间内得到了解决，其中 (omega) 是矩阵的大小。

（在下面的代码实现中 (omega = 4)（矩阵构造详见代码），但存在转移矩阵的大小为 (3 imes 3) 的做法）

代码

#include<bits/stdc++.h>

using namespace std;

#define rg register

typedef long long ll;

template<typename T> inline bool checkMax(T& a, const T& b) {
  return a < b ? a = b, true : false;
}

const int N = 1e5 + 10;
const double eps = 1e-12;

int n, a[N], b[N], top;
ll times;
double p[N], m;

struct Point {
  double x, y;
  Point () {}
  Point (double x, double y): x(x), y(y) {}
  bool operator < (const Point& a) const {
    return x == a.x ? y < a.y : x < a.x;
  }
} points[N], _stack[N];

inline double slope(Point a, Point b) {
  return fabs(a.x - b.x) < eps ? 1e100 : (b.y - a.y) / (b.x - a.x);
}

struct Matrix {
  double a[4][4];
  Matrix () {
    memset(a, 0, sizeof a);
  }
  Matrix operator * (const Matrix& b) const {
    Matrix res;
    for (rg int i = 0; i < 4; ++i) {
      for (rg int j = 0; j < 4; ++j) {
        for (rg int k = 0; k < 4; ++k) {
          res.a[i][j] += a[i][k] * b.a[k][j];
        }
      }
    }
    return res;
  }
} mat, trans[40];

int main() {
  scanf("%d%I64d", &n, &times);
  for (rg int i = 1; i <= n; ++i) {
    scanf("%d%d%lf", &a[i], &b[i], &p[i]);
    points[i] = Point (-p[i], a[i] * p[i]);
    checkMax(m, p[i] * b[i]);
  }
  sort(points + 1, points + 1 + n);
  _stack[top = 1] = points[1];
  for (rg int i = 2; i <= n; ++i) {
    for (; top > 1 && slope(_stack[top], points[i]) >= slope(_stack[top - 1], _stack[top]); --top);
    _stack[++top] = points[i];
  }
  ll now = 0;
  trans[0].a[1][0] = 1;
  trans[0].a[1][1] = 1;
  trans[0].a[2][0] = 1;
  trans[0].a[2][2] = 1;
  trans[0].a[3][2] = 1;
  trans[0].a[3][3] = 1;
  for (rg int i = top; i && now ^ times; --i) {
    double r = i == 1 ? 1e100 : slope(_stack[i - 1], _stack[i]);
    double x = -_stack[i].x, y = _stack[i].y;
    mat.a[0][1] = y;
    mat.a[0][2] = now * x * m;
    mat.a[0][3] = x * m;
    trans[0].a[0][0] = 1 - x;
    if (now * m - mat.a[0][0] > r) {
      continue;
    }
    for (rg int j = 1; j < 34; ++j) {
      trans[j] = trans[j - 1] * trans[j - 1];
    }
    for (rg int j = 33; ~j; --j) {
      if (now + (1ll << j) >= times) {
        continue;
      }
      Matrix old = mat;
      mat = mat * trans[j];
      double k = (now + (1ll << j)) * m - mat.a[0][0];
      if (k >= r) {
        mat = old;
        continue;
      }
      now += 1ll << j;
    }
    ++now;
    mat = mat * trans[0];
  }
  printf("%.15lf
", mat.a[0][0]);
  return 0;
}