「集训队作业2020」Old Problem

2020 年 12 月 17 日

给一个长度为 $n$ 的序列 $a_i$ ，和 $q$ 组询问 $(l,r,x)$ ，表示求 $\displaystyle\prod_{i=l}^r\left(1-\frac{a_i}{x}\right)$ 的值。实数输出，精度要求 $10^{-6}$ 。

$n,q\le6\times10^5,\ 1\leq a_i < x\leq 10^9$ 。

题解

按照 EI 的话说，这是一个误差分析题。

首先需要注意到这个式子可以泰勒展开：

\begin{aligned} &\exp\left(\ln\left(\prod_{i=l}^r\left(1-\frac{a_i}{x}\right)\right)\right) \\ =&\exp\left(\sum_{i=l}^r\ln\left(1-\frac{a_i}{x}\right)\right) \\ =&\exp\left(-\sum_{i=l}^r\sum_{k=1}^{L}\frac{a_i^k}{k\cdot x^k}\right) \\ \end{aligned}

然而，如果直接泰勒展开，需要的 $L$ 是数十万级别的，无益于我们解决问题。

考虑到导致泰勒展开精度损失的主要原因，是因为当 $a_i/x$ 较大时，我们对 $\ln(1-a_i/x)$ 的精度要求很高。然而， $a_i/x$ 较大时，很容易导致答案小于我们要求的精度范围。

故我们不妨设定一个阈值 $R=0.5$ ，当 $a_i/x\le R$ 时，考虑线段树维护泰勒展开；否则，当 $a_i/x>R$ 时，优先找出这些位置并暴力计算。当答案小于精度要求时就退出。暴力计算的次数显然不会超过 $\log_2 10^6$ 次。

这样就得到了一个 $O(n\log^2L)$ 的做法，实践得 $L$ 取 $20$ 左右即可。

可以通过本题，但时间较大。实际上，线段树的部分可以直接换为前缀和。为什么精度还在接受范围内呢？注意到 $x>a_i$ 对于任意 $x$ 和任意 $i$ 都成立。如果前缀和的部分因为 $a_i$ 太小被省略，他本身对泰勒展开的影响也是被省略的级别。换句话说，对于泰勒展开的值，前缀和能保证的精度范围，恰为 double 本身的精度范围。

所以，直接将上述做法中的线段树替换为前缀和就能在 $O(n\log nL)$ 的时间复杂度内解决本题。~~三个 log 年轻人不讲武德。~~

代码

#include <bits/stdc++.h>
using namespace std;
const int N = 6e5 + 10, K = 20;
const double eps = 1e-10;
int n, m, x, a[N], lg[N], st[20][N];
double sum, ans;
mt19937 rng(20040725);
inline int better(int i, int j) { return a[i] < a[j] ? j : i; }
inline int query(int l, int r) {
  if (l == r) return l;
  int k = lg[r - l];
  return better(st[k][l], st[k][r - (1 << k)]);
}
struct segment {
  int l, r, mid;
  double s[K];
} p[N << 2];
void build(int u, int l, int r) {
  p[u].l = l, p[u].r = r, p[u].mid = (l + r) >> 1;
  if (l == r) {
    p[u].s[0] = a[l];
    for (int i = 1; i < K; i++) p[u].s[i] = p[u].s[i - 1] * a[l];
    return;
  }
  build(u << 1, l, p[u].mid);
  build(u << 1 | 1, p[u].mid + 1, r);
  for (int i = 0; i < K; i++) {
    p[u].s[i] = p[u << 1].s[i] + p[u << 1 | 1].s[i];
  }
}
double query(int u, int l, int r) {
  if (p[u].l == l && p[u].r == r) {
    double sum = 0;
    for (int i = K - 1; i >= 0; i--) {
      sum = (sum + p[u].s[i] / (i + 1)) / x;
    }
    return sum;
  }
  if (r <= p[u].mid) return query(u << 1, l, r);
  if (l > p[u].mid) return query(u << 1 | 1, l, r);
  return query(u << 1, l, p[u].mid) + query(u << 1 | 1, p[u].mid + 1, r);
}
void solve(int l, int r) {
  if (l > r || ans < eps) {
    return;
  }
  int pos = query(l, r);
  if (a[pos] < (x >> 1)) {
    ans *= exp(-query(1, l, r));
    return;
  }
  ans *= 1 - a[pos] / (double)x;
  if (rng() & 1) {
    solve(l, pos - 1);
    solve(pos + 1, r);
  } else {
    solve(pos + 1, r);
    solve(l, pos - 1);
  }
}
int main() {
#ifdef memset0
  freopen("1.in", "r", stdin);
#endif
  cin.tie(0)->sync_with_stdio(0);
  lg[0] = -1;
  for (int i = 1; i < N; i++) lg[i] = lg[i >> 1] + 1;
  cin >> n >> m;
  for (int i = 1; i <= n; i++) cin >> a[i];
  for (int i = 1; i < n; i++) {
    st[0][i] = better(i, i + 1);
  }
  for (int i = 1; i < 20; i++)
    for (int j = 1; j + (1 << i) <= n; j++) {
      st[i][j] = better(st[i - 1][j], st[i - 1][j + (1 << (i - 1))]);
    }
  build(1, 1, n);
  cout << fixed << setprecision(12);
  for (int l, r, i = 1; i <= m; i++) {
    cin >> l >> r >> x;
    ans = 1;
    solve(l, r);
    cout << 1 - ans << endl;
  }
}

TABLE OF CONTENTS

题解

代码

题解

代码

评论