feat(solve): monte_carlo_integral 落地，修掉体积符号真 bug

实现由用户本人写；Codex 只碰工程侧 + 一处纯 bug（决策 31）。

- 真 bug：体积写成 lows - highs，奇数维为负（1D x^2 于 [0,1] 得 -0.3324、3D 得 -0.5005）；
  2D 因两个负号相乘"看着对"（+0.2499），把 bug 藏住。改成 widths = highs - lows 后
  复测 0.3324 / 0.5005 / 1.9971（[0,2]）全部正确
- docstring：补 func / n / random_state 与 Raises；写清 metrics["mean"] 才是积分估计值、
  data 是 n 个样本各自的贡献（决策 51）；修掉跑不通的 Examples（data 是数组，bool(数组) 会抛）
- tests/test_solve.py：摘 1 条 skip + 新增 4 条（3D 奇数维 / [0,2] 宽度 / 2D+params / 非法 bounds）
- 验证：197 收集 / 166 passed / 31 skipped / 0 xfail，ruff 0，pyright 0
- docs：STATUS 决策 51 + 3 条待办 + 下一站 random_walk；CHANGELOG；易错点 22/23，第 3 条复发 +1

src/mcm_toolkit/solve/simulate.py

@@ -213,7 +213,7 @@ def monte_carlo_probability(
 
         warning_message = ""
         if p_value <= 1e-3 or p_value >= 1 - 1e-3:
-            warning_message = "小概率事件 P < 1e-3 不建议使用monte_carlo_probability"
+            warning_message = "小概率事件 P < 1e-3 或者高概率事件 P > 1 - 1e-3,不建议使用monte_carlo_probability"
 
         metrics = {
             "probability": p_value,
@@ -253,8 +253,26 @@ def monte_carlo_integral(
 
     Parameters
     ----------
+    func:
+        `f(X) -> 一批函数值`。`X` 是形状 `(n, d)` 的样本矩阵（一行一个样本），
+        返回长度 `n` 的数组（写成 `X[:, 0] ** 2` 这种逐行写法，不要写 Python 循环）。
     bounds:
-        `[(lo, hi), ...]`，维度即积分重数。
+        `[(lo, hi), ...]`，维度即积分重数，要求每个 `lo < hi`。
+    n:
+        模拟次数，正整数。
+    random_state:
+        随机种子。`None` 时跟随 `core.set_seed` 的全局流，见 `core.random.get_rng`。
+
+    Returns
+    -------
+    `Result`：**`metrics["mean"]` 才是积分估计值** ≈ `∫f dV`；
+    `data` 是 `n` 个样本各自的贡献（`func(X_i) × 体积`），可拿去画随 `n` 的收敛图。
+    `params` 里额外记了 `bounds` 与 `dimension`。
+
+    Raises
+    ------
+    ValueError:
+        `bounds` 里出现 `lo >= hi`。
 
     Notes
     -----
@@ -265,11 +283,42 @@ def monte_carlo_integral(
     --------
     >>> import numpy as np
     >>> r = monte_carlo_integral(lambda X: X[:, 0] ** 2, bounds=[(0, 1)], n=200000, random_state=0)
-    >>> abs(r.data - 1 / 3) < 0.01
+    >>> abs(r.metrics["mean"] - 1 / 3) < 0.01   # ∫₀¹ x² dx = 1/3
     True
     """
-    # TODO(L1)  ★ S 档：配完整函数卡片
-    raise NotImplementedError
+    bounds_arr = np.asarray(bounds)
+
+    lows = bounds_arr[:, 0]
+    highs = bounds_arr[:, 1]
+    d = len(lows)
+    if np.any(lows >= highs):
+        raise ValueError("积分边界非法：下界必须严格小于上界")
+
+    widths = highs - lows
+
+    volume = np.prod(widths)
+
+    def sampler(rng, size=n):
+        return rng.uniform(low=lows, high=highs, size=(size, d))
+
+    def simulator(X):
+        vals = func(X)
+        if len(vals) != n:
+            raise ValueError(f"func 返回的结果数量({len(vals)})必须与模拟次数 n({n})一致")
+        return np.asarray(vals, dtype=float).ravel() * volume
+
+    result = monte_carlo(
+        sampler=sampler,
+        simulator=simulator,
+        n=n,
+        random_state=random_state,
+    )
+    result.name = "monte_carlo_integral"
+    result.params["bounds"] = bounds
+    result.params["dimension"] = d
+    result.message = f"蒙特卡洛多重积分成功：维度 {d}，体积 {volume:.4f}，模拟 {n} 次。"
+
+    return result
 
 
 def random_walk(