Problem
题意概要：给定一个有重边有自环 (n) 点 (m) 边的有向无环图（DAG），每条边有其权值，每当你走到一个点 (x) 时，所有从 (x) 连出去的边上的权值会互相随机打乱，问从 (S) 到 (T) 最短路长度的期望
(n,mleq 10^3)
Solution
首先第一步很明显是按照 DAG 的拓扑序一个个地转移，只需考虑处理每个点怎么转移，设 (f[x]) 表示从 (x) 走到 (T) 的最短路长度期望
先暂不考虑重边和自环
设当前在点 (x)，且共有 (K) 条出边，那么 (f[x]) 应该会有 (K^2) 种取值（(K) 种出边边权和 (K) 个 (f[v](x ightarrow v)) 值自由搭配）
考虑计算若当前选择了 边权 (w) 与 出点(f[v](x ightarrow v)) 搭配，需要计算 (w+f[v]) 作为所有 (K) 种搭配下最小的权值的概率
对于 (x) 的每一个不同于 (v) 的出点 (t)，需要计算出有多少边权与其搭配使得这组搭配花费总和不小于 (w) 与 (v) 的这组搭配，记作 (ret[t])
一个结论是 (ret[t]) 与 (f[t]) 正相关（因为若 (f[t]) 越大，则给边权 (t_w) 的限制越小，(f[t]+t_wgeq w+f[v]Leftrightarrow t_wgeq w+f[v]-f[t])）。而对于 (f[i]leq f[j])，定有 (ret[i] leq ret[j])，且两者集合之间为包含关系（后者包含前者）
则考虑将所有出点按照 (f) 值从小至大排序后，(t) 处在第 (i) 位，则减去其包含的区间，这一位能选的配对有 (ret[t]-(i-1)) 种
所以边权 (w) 与出点 (f[v]) 的配对为 (K) 个配对最小的情况，有 (prod_t (ret[t]-i+1)) 种，这样就能在 (O(K^3)) 的时间内计算出每个点的 (f) 值；事实上，按照配对的权值排序，计算增量出现的概率，每一次只会改变两个值，可以 (O(1)) 解决，算上排序，可以在 (O(K^2log K)) 的时间内计算出每个点的 (f) 值
考虑上重边，发现 (x) 若有 (k) 条边指向 (v)，则只要假定 (x) 有 (k) 个 (f) 值为 (f[v]) 的出点即可
考虑上自环，可以使用二分 (f) 的方式，每次假定 (f=mid)，和由此算出来的 (f') 进行对比，若 (f<f') 则代表 (mid) 设得比较小，反之则越大（事实上二分中也可以用一点点小优化，若 (f<f') 则将 (l) 设为 (f') 而非 (mid)，这样实际上算一种小迭代？速度上升一倍）
加上二分后就需要将对配对的排序提到外头来，内层仅处理自环的部分后线性归并一下，否则复杂度会变成两个 (log)
时间复杂度为 (O(nlog p+m^2log m)) （(p) 视精度要求而定）
Code
#include <bits/stdc++.h>
using namespace std;
#define For(x,y) for(int x=1;x<=y;++x)
const double eps = 1e-8;
const int N = 1010;
struct Edge {int v, nxt; double w;} a[N+N];
int head[N], Head[N], _;
int deg[N], q[N];
double f[N];
int n, m, S, T;
inline void ad() {
static int x, y; static double w; scanf("%d%d%lf",&x,&y,&w);
a[++_].v = y, a[_].w = w, a[_].nxt = head[x], head[x] = _;
a[++_].v = x, a[_].nxt = Head[y], Head[y] = _;
if(x != y) ++deg[x];
}
typedef pair<double,int> pr;
pr brr[N*N], arr[N*N], crr[N*N];
int arc, brc, crc;
double ew[N*N]; int _ew;
double ot[N*N]; int _ot;
int ret[N], self_circle;
void init(int x) {
arc = _ew = _ot = self_circle = 0;
for(int i=head[x];i;i=a[i].nxt) {
ew[++_ew] = a[i].w;
if(a[i].v == x) ++self_circle;
else ot[++_ot] = f[a[i].v];
}
sort(ew + 1, ew + _ew + 1);
sort(ot + 1, ot + _ot + 1);
For(i, _ew) For(j, _ot)
arr[++arc] = make_pair(ew[i] + ot[j], j);
sort(arr + 1, arr + arc + 1);
}
double calc(int x, double sw) {
f[x] = sw, brc = 0;
int t = 0; while(t < _ot and ot[t+1] < sw) ++t;
For(i, _ew) For(j, self_circle)
brr[++brc] = make_pair(ew[i] + sw, t + j);
For(i, arc) if(arr[i].second > t) arr[i].second += self_circle;
merge(arr+1, arr+arc+1, brr+1, brr+brc+1, crr+1);
For(i, arc) if(arr[i].second > t) arr[i].second -= self_circle;
crc = arc + brc;
For(i, _ew) ret[i] = _ew;
double coe = 1, ans = 0;
for(int i=1, id; i <= crc and coe > eps; ++i) {
ans += coe * (crr[i].first - crr[i-1].first);
id = crr[i].second;
coe /= ret[id] - (id - 1);
--ret[id];
coe *= ret[id] - (id - 1);
}
return ans;
}
void solve(int x) {
init(x);
if(!_ot) return f[x] = 1e9, void();
double l = calc(x, 0), r = calc(x, 1e7), mid, res;
while(l + eps < r) {
mid = 0.5 * (l + r);
res = calc(x, mid);
if(res > mid) l = res;
else r = res;
}
f[x] = l;
}
int main() {
scanf("%d%d%d%d",&n,&m,&S,&T);
while(m--) ad();
int he = 1, ta = 0;
For(i, n) {
if(!deg[i]) q[++ta] = i;
f[i] = 1e9;
}
while(he <= ta) {
int x = q[he++];
for(int i=Head[x];i;i=a[i].nxt)
if(!(--deg[a[i].v]))
q[++ta] = a[i].v;
if(x == T) f[x] = 0;
else solve(x);
}
if(f[S] < 1e8) printf("%.7lf
",f[S]);
else puts("-1");
return 0;
}