题解 - [SDOI2017] 龙与地下城

发表于 更新于 分类于 ， 本文字数： 1.2k 阅读时长 ≈ 4 分钟

很有 OI 味的一道题

题目链接

原始题面

题目描述

小 Q 同学是一个热爱学习的人，但是他最近沉迷于各种游戏，龙与地下城就是其中之一

在这个游戏中，很多场合需要通过掷骰子来产生随机数，并由此决定角色未来的命运，因此骰子堪称该游戏的标志性道具

骰子也分为许多种类，比如 4 面骰、6 面骰、8 面骰、12 面骰、20 面骰，其中 20 面骰用到的机会非常多。当然，现在科技发达，可以用一个随机数生成器来取代真实的骰子，所以这里认为骰子就是一个随机数生成器

在战斗中，骰子主要用来决定角色的攻击是否命中，以及命中后造成的伤害值。举个例子，假设现在已经确定能够命中敌人，那么 \(YdX\) (也就是掷出 \(Y\) 个 \(X\) 面骰子之后所有骰子显示的数字之和) 就是对敌人的基础伤害。在敌人没有防御的情况下，这个基础伤害就是真实伤害

众所周知，骰子显示每个数的概率应该是相等的，也就是说，对于一个 \(X\) 面骰子，显示 \(0, 1, 2,\dots ,X−1\) 中每一个数字的概率都是 \(\frac {1}{x}\)

更形式地说，这个骰子显示的数 \(W\) 满足离散的均匀分布，其分布列为

\(W\) \(0\) \(1\) \(2\) \(\dots\) \(X-1\)
\(P\) \(\frac{1}{X}\) \(\frac{1}{X}\) \(\frac{1}{X}\) \(\dots\) \(\frac{1}{X}\)

除此之外还有一些性质

  • \(W\) 的一阶原点矩 (期望) 为 \(v_1(W)=E(W)=\sum_{i=0}^{X-1}iP(W=i)=\frac {X-1}{2}\)

  • \(W\) 的二阶中心矩 (方差) 为 \(\mu_2(W)=E((W-E(W))^2)=\sum_{i=0}^{X-1}(i-E(W))^2P(W=i)=\frac {X^2-1}{12}\)

言归正传，现在小 Q 同学面对着一个生命值为 A 的没有防御的敌人，能够发动一次必中的 \(YdX\) 攻击，显然只有造成的伤害不少于敌人的生命值才能打倒敌人。但是另一方面，小 Q 同学作为强迫症患者，不希望出现 overkill, 也就是造成的伤害大于 \(B\) 的情况，因此只有在打倒敌人并且不发生 overkill 的情况下小 Q 同学才会认为取得了属于他的胜利

因为小 Q 同学非常谨慎，他会进行 10 次模拟战，每次给出敌人的生命值 \(A\) 以及 overkill 的标准 \(B\), 他想知道此时取得属于他的胜利的概率是多少，你能帮帮他吗？

输入格式

第一行是一个正整数 \(T\), 表示测试数据的组数，

对于每组测试数据:

第一行是两个整数 \(X\) , \(Y\), 分别表示骰子的面数以及骰子的个数；

接下来 10 行，每行包含两个整数 \(A\) , \(B\), 分别表示敌人的生命值 \(A\) 以及 overkill 的标准 \(B\)

输出格式

对于每组测试数据，输出 10 行，对每个询问输出一个实数，要求绝对误差不超过 \(0.013579\),

也就是说，记输出为 \(a\), 答案为 \(b\), 若满足 \(|a-b|\leq 0.013579\), 则认为输出是正确的

样例 #1

样例输入 #1

1
2
3
4
5
6
7
8
9
10
11
12
1
2 19
0 0
0 1
0 2
0 3
0 4
0 5
0 6
0 7
0 8
0 9

样例输出 #1

1
2
3
4
5
6
7
8
9
10
0.000002
0.000038
0.000364
0.002213
0.009605
0.031784
0.083534
0.179642
0.323803
0.500000

提示

对于 \(100\%\) 的数据，\(T \leq 10\), \(2 \leq X \leq 20\), \(1 \leq Y \leq 200000\), \(0 \leq A \leq B \leq (X-1)Y\), 保证满足 \(Y > 800\) 的数据不超过 \(2\) 组

测试点编号 \(X\) \(Y\)
1 \(\leq20\) \(\leq40\) \(X^Y\leq10000000\)
2, 3, 4 \(\leq20\) \(\leq1600\)
5, 6, 7, 8, 9, 10 \(\leq20\) \(\leq8000\)
11, 12 \(=2\) \(\leq200000\)
13, 14, 15, 16, 17, 18, 19, 20 \(\leq20\) \(\leq200000\)

解题思路

对于单次询问，答案显然是

\[ \sum_{k=A}^B[x^k]\left(\sum_{i=0}^{X-1}\frac{1}{X}x^i\right)^Y \]

这个时间复杂度是 \(O(XY\log(XY))\), 对于大数据来说不能接受

注意到本题对误差要求极为宽泛，这提示我们寻找一些拟合的方法

由中心极限定理 , 设 \(n\) 个独立同分布的随机变量 \(x_1,x_2,...,x_n\) 的均值 \(E(x_i)=\mu\), 方差 \(E((x_i-E(x_i))^2)=\sigma^2\ne 0\), 则随机变量 \(Y=\sum_{i=1}^nx_i\) 的分布在 \(n\to\infty\) 时服从 \(N(n\mu,n\sigma^2)\)

所以在 \(XY\) 较大时我们可以用正态分布去拟合结果，求解时直接上自适应 Simpson 积分即可

代码参考

Show code

Luogu_P3779view raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
/*
 * @Author: Tifa
 * @Description: From <https://github.com/Tiphereth-A/CP-archives>
 * !!! ATTENEION: All the context below is licensed under a 
 *  GNU Affero General Public License, Version 3.
 *  See <https://www.gnu.org/licenses/agpl-3.0.txt>.
 */
#include <bits/stdc++.h>
#define for_(i, l, r, vars...) \
  for (decltype(l + r) i = (l), i##end = (r), ##vars; i <= i##end; ++i)
using namespace std;
class SimpsonSolver {
protected:
  const double eps__;
  const int least_dep__;
  template <typename Func>
  constexpr double simpson(double l, double r, Func &&func) const {
    return (func(l) + 4 * func(r - (r - l) * .5) + func(r)) * (r - l) / 6;
  }
  template <typename Func>
  double
  asr(double l, double r, double area, double eps, int dep, Func &&func) const {
    double mid = r - (r - l) * .5;
    double ls = simpson(l, mid, func), rs = simpson(mid, r, func);
    if (std::abs(ls + rs - area) <= 15 * eps && dep < 0)
      return ls + rs + (ls + rs - area) / 15;
    return asr(l, mid, ls, eps * .5, dep - 1, func) +
           asr(mid, r, rs, eps * .5, dep - 1, func);
  }

public:
  SimpsonSolver(double eps, int least_dep)
    : eps__(eps), least_dep__(least_dep) {}
  template <typename Func>
  double operator()(double l, double r, Func &&func) const {
    return asr(l, r, simpson(l, r, func), eps__, least_dep__, func);
  }
} simpson_solver(1e-6, 8);
namespace FFT {
using CP = std::complex<double>;
const double PI = acos(-1.0);
size_t n = 0;
std::vector<size_t> rev;
void init(size_t m) {
  if (n > m) return;
  n = 1;
  int k = 0;
  while (n <= m) {
    n <<= 1;
    ++k;
  }
  rev.resize(n);
  for (int i = 0; i < n; ++i)
    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
}
void dft(std::vector<CP> &a, int op) {
  for (int i = 0; i < n; ++i)
    if (rev[i] > i) std::swap(a[rev[i]], a[i]);
  for (int i = 1; i < n; i <<= 1) {
    CP now(std::cos(PI / i), op * std::sin(PI / i));
    for (int j = 0, p = (i << 1); j < n; j += p) {
      CP w(1, 0);
      for (int k = 0; k < i; ++k, w = w * now) {
        CP x = a[j + k], y = a[j + k + i] * w;
        a[j + k] = x + y;
        a[j + k + i] = x - y;
      }
    }
  }
  if (op == -1)
    for (int i = 0; i < n; ++i) a[i] = {a[i].real() / n, a[i].imag() / n};
}
std::vector<double> conv(std::vector<double> const &a,
                         std::vector<double> const &b) {
  init(a.size() + b.size() - 1);
  std::vector<CP> a_(n), b_(n), c_(n);
  for (size_t i = 0; i < a.size(); ++i) a_[i] = a[i];
  for (size_t i = 0; i < b.size(); ++i) b_[i] = b[i];
  dft(a_, 1);
  dft(b_, 1);
  for (size_t i = 0; i < n; ++i) c_[i] = a_[i] * b_[i];
  dft(c_, -1);
  std::vector<double> c(a.size() + b.size() - 1);
  for (size_t i = 0; i < a.size() + b.size() - 1; ++i) c[i] = c_[i].real();
  return c;
}
}  // namespace FFT
using FFT::conv;
std::vector<double> qpow(std::vector<double> a, uint64_t b) {
  std::vector<double> res{1};
  for (; b; b >>= 1, a = conv(a, a))
    if (b & 1) res = conv(res, a);
  return res;
};
const uint32_t OFFSET = 5;
const uint32_t N = 3e5 + OFFSET;
auto solve([[maybe_unused]] int t_ = 0) -> void {
  int x, y;
  cin >> x >> y;
  if (x * y < N) {
    vector<double> f = qpow(vector<double>(x, 1. / x), y);
    for_(i, 1, f.size() - 1) f[i] += f[i - 1];
    for_(i, 1, 10, a, b) {
      cin >> a >> b;
      cout << f[b] - (a ? f[a - 1] : 0) << '\n';
    }
  } else {
    const double ymu = (x - 1) * .5 * y,
                 ystdi = 1. / sqrt((x * x - 1) / 12. * y);
    auto norm = [ymu, ystdi](double x) -> double { return (x - ymu) * ystdi; };
    auto f = [](double x) -> double {
      return exp(-x * x * .5) / sqrt(2 * acos(-1.));
    };
    for_(i, 1, 10, a, b) {
      cin >> a >> b;
      cout << simpson_solver(norm(a), norm(b), f) << '\n';
    }
  }
}
int main() {
  std::ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  int i_ = 0;
  std::cout << std::fixed << std::setprecision(12);
  int t_ = 0;
  std::cin >> t_;
  for (i_ = 0; i_ < t_; ++i_) solve(i_);
  return 0;
}