CS Academy Round #10 - Expected Tree Degrees
文章目录
【注意】最后更新于 August 14, 2020，文中内容可能已过时，请谨慎使用。
题目描述
有一棵 $n$ 个点的树，由如下方式生成：
对于每个点 $i$ ($i\ge 2$)，在 $1$ 到 $i-1$ 中等概率选择一个点作为点 $i$ 的父亲。
求每个点的度数的平方和的期望，要求绝对误差不超过 $10^{-6}$，$n\le 10^6$。
状态设计与转移
令 $f_{i, j}$ 表示 $i$ 个点的树中度数为 $j$ 的点的个数的期望。
那么：$f_{i, j} = f_{i - 1, j} - \frac{f_{i - 1, j}}{i - 1} + \frac{f_{i - 1, j - 1}}{i - 1} + [j=1]$
其中 $\frac{f_{i - 1, j}}{i - 1}$ 表示点 $i$ 的父亲在只有前 $i-1$ 个点时度数为 $j$ 的概率，$[j=1]$ 表示 j == 1 ? 1 : 0。
然而，直接 DP 复杂度是 $O(n^2)$，无法通过此题。有多种处理方式，见下文。
求答案的近似值
可以猜测 $j$ 比较大时 $f_{i, j}$ 非常小，可以忽略不计。实际上 $f_{10^6, 50}$ 只有大约 $10^{-14}$，所以只用计算 $j\le 50$ 的 $f_{i, j}$。
用生成函数和矩阵求答案的精确值
虽然上面这种求近似值的方法可以通过本题，但本题实际上是可以线性地求答案的精确值的。
首先，每个 $f_{i, 1..i}$ 可以视为一个生成函数，记作 $F_i$，也就是说 $F_i$ 的 $x^j$ 项的系数表示 $f_{i, j}$。
那么，转移方程可以表示为：
$$ \begin{aligned} F_i &= F_{i - 1} - \frac{F_{i - 1}}{i - 1} + x\frac{F_{i - 1}}{i - 1}+x\\ &= \frac{i-2+x}{i-1}F_{i - 1} + x \end{aligned} $$
然而，这样算出 $F_n$ 的每一项的系数依然是 $O(n^2)$ 的。
但是，我们要求的并不是 $F_n$ 的每一项的系数，而是每一项的系数乘上次数的平方的和。
注意到 $\{1^2, 2^2, 3^2, \ldots\}$ 这个数列可以用下面这个矩阵的幂计算：
$$ A= \begin{bmatrix} 1&0&0\\ 1&1&0\\ 1&2&1 \end{bmatrix} $$
$$ A \begin{bmatrix} 1\\ i\\ i^2 \end{bmatrix}= \begin{bmatrix} 1\\ (i+1)\\ (i+1)^2 \end{bmatrix} $$
所以，把 $A$ 作为 $x$ 带入生成函数中计算即可得到答案。
参考代码
// Problem: Expected Tree Degrees
// Contest: CS Academy
// URL: https://csacademy.com/contest/archive/task/expected-tree-degrees/
// Memory Limit: 128 MB
// Time Limit: 1000 ms
// Powered by CP Editor (https://github.com/cpeditor/cpeditor)
#ifndef OUUAN
#pragma GCC optimize("Ofast,no-stack-protector,unroll-loops,fast-math")
#endif
#include<bits/stdc++.h>
#define int ll
// #define FAST_IOSTREAM 1
#define For(i,l,r)for(int i=(l),i##end=(r);i<=i##end;++i)
#define FOR(i,r,l)for(int i=(r),i##end=(l);i>=i##end;--i)
#define SON(i,u)for(int i=head[u];i;i=nxt[i])
#define BE(x)(x).begin(),(x).end()
#define BE1(x)((x).begin()+1),(x).end()
#define fi first
#define se second
#define pb push_back
#define eb emplace_back
#define pq priority_queue
#define min minOfDifferentTypes
#define max maxOfDifferentTypes
#define y1 why_is_there_a_function_called_y1
using namespace std;typedef long long ll;typedef pair<int,int>pii;typedef vector<int>vi;typedef vector<vi>vvi;typedef vvi v2i;typedef vector<vvi>v3i;typedef vector<v3i>v4i;typedef vector<bool>vb;typedef vector<vb>vvb;typedef vvb v2b;typedef vector<vvb>v3b;typedef vector<v3b>v4b;typedef vector<pii>vpii;typedef vector<vpii>vvpii;typedef vvpii v2pii;typedef vector<vvpii>v3pii;typedef vector<v3pii>v4pii;typedef vector<double>vd;typedef vector<vd>vvd;typedef vvd v2d;typedef vector<v2d>v3d;typedef vector<v3d>v4d;const double inf=1e121;const double eps=1e-10;mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());int randint(int l,int r){int out=rng()%(r-l+1)+l;return out>=l?out:out+r-l+1;}template<typename A,typename B>string to_string(pair<A,B>p);template<typename A,typename B,typename C>string to_string(tuple<A,B,C>p);template<typename A,typename B,typename C,typename D>string to_string(tuple<A,B,C,D>p);string to_string(const string&s){return'"'+s+'"';}string to_string(const char*s){return to_string((string)s);}string to_string(bool b){return(b?"true":"false");}string to_string(vector<bool>v){bool first=true;string res="{";for(int i=0;i<static_cast<int>(v.size());i++){if(!first){res+=", ";}first=false;res+=to_string(v[i]);}res+="}";return res;}template<size_t N>string to_string(bitset<N>v){string res="";for(size_t i=0;i<N;i++){res+=static_cast<char>('0'+v[i]);}return res;}template<typename A>string to_string(A v){bool first=true;string res="{";for(const auto&x:v){if(!first){res+=", ";}first=false;res+=to_string(x);}res+="}";return res;}template<typename A,typename B>string to_string(pair<A,B>p){return"("+to_string(p.first)+", "+to_string(p.second)+")";}template<typename A,typename B,typename C>string to_string(tuple<A,B,C>p){return"("+to_string(get<0>(p))+", "+to_string(get<1>(p))+", "+to_string(get<2>(p))+")";}template<typename A,typename B,typename C,typename D>string to_string(tuple<A,B,C,D>p){return"("+to_string(get<0>(p))+", "+to_string(get<1>(p))+", "+to_string(get<2>(p))+", "+to_string(get<3>(p))+")";}template<typename A,typename B,typename C,typename D,typename E>string to_string(tuple<A,B,C,D,E>p){return"("+to_string(get<0>(p))+", "+to_string(get<1>(p))+", "+to_string(get<2>(p))+", "+to_string(get<3>(p))+","+to_string(get<4>(p))+")";}void debug_out(){cerr<<endl;}template<typename Head,typename...Tail>void debug_out(Head H,Tail...T){cerr<<" "<<to_string(H);debug_out(T...);}template<typename T>struct is_pair{static const bool value=false;};template<typename T,typename U>struct is_pair<std::pair<T,U>>{static const bool value=true;};
#ifdef OUUAN
#define debug(...)cerr<<"["<<#__VA_ARGS__<<"]:",debug_out(__VA_ARGS__)
#else
#define debug(...)42
#endif
#ifdef int
const int INF=0x3f3f3f3f3f3f3f3fll;
#else
const int INF=0x3f3f3f3f;
#endif
#ifdef FAST_IOSTREAM
#define br cout<<'\n'
#define sp cout<<' '
#define fl cout.flush()
ll read(){ll x;cin>>x;return x;}template<typename T>void read(T&x){cin>>x;}template<typename T>void write(const T&x){cout<<x;}
#else
#define br putchar('\n')
#define sp putchar(' ')
#define fl fflush(stdout)
template<typename T>typename enable_if<!is_integral<T>::value&&!is_pair<T>::value,void>::type read(T&x){cin>>x;}ll read(){char c;ll out=0,f=1;for(c=getchar();!isdigit(c)&&c!='-';c=getchar()){}if(c=='-'){f=-1;c=getchar();}for(;isdigit(c);c=getchar())out=(out<<3)+(out<<1)+c-'0';return out*f;}template<typename T>typename enable_if<is_integral<T>::value,T>::type read(T&x){char c;T f=1;x=0;for(c=getchar();!isdigit(c)&&c!='-';c=getchar()){}if(c=='-'){f=-1;c=getchar();}for(;isdigit(c);c=getchar())x=(x<<3)+(x<<1)+c-'0';return x*=f;}char read(char&x){for(x=getchar();isspace(x);x=getchar()){}return x;}double read(double&x){scanf("%lf",&x);return x;}template<typename T>typename enable_if<!is_integral<T>::value&&!is_pair<T>::value,void>::type write(const T&x){cout<<x;}template<typename T>typename enable_if<is_integral<T>::value,void>::type write(const T&x){if(x<0){putchar('-');write(-x);return;}if(x>9)write(x/10);putchar(x%10+'0');}void write(const char&x){putchar(x);}void write(const double&x){printf("%.12lf",x);}
#endif
template<typename T>typename enable_if<is_pair<T>::value,void>::type read(T&x){read(x.fi);read(x.se);}template<typename T>typename enable_if<is_pair<T>::value,void>::type write(const T&x){write(x.fi);sp;write(x.se);}template<typename T,typename...Args>void read(T&x,Args&...args){read(x);read(args...);}template<typename OutputIt,typename=typename enable_if<is_same<output_iterator_tag,typename iterator_traits<OutputIt>::iterator_category>::value||(is_base_of<forward_iterator_tag,typename iterator_traits<OutputIt>::iterator_category>::value&&!is_const<OutputIt>::value)>::type>void read(OutputIt __first,OutputIt __last){for(;__first!=__last;++__first)read(*__first);}template<typename InputIt,typename=typename enable_if<is_base_of<input_iterator_tag,typename iterator_traits<InputIt>::iterator_category>::value>::type>void wts(InputIt __first,InputIt __last){bool isFirst=true;for(;__first!=__last;++__first){if(isFirst)isFirst=false;else sp;write(*__first);}br;}template<typename InputIt,typename=typename enable_if<is_base_of<input_iterator_tag,typename iterator_traits<InputIt>::iterator_category>::value>::type>void wtb(InputIt __first,InputIt __last){for(;__first!=__last;++__first){write(*__first);br;}}template<typename T>void wts(const T&x){write(x);sp;}template<typename T>void wtb(const T&x){write(x);br;}template<typename T>void wte(const T&x){write(x);exit(0);}template<typename T,typename...Args>void wts(const T&x,Args...args){wts(x);wts(args...);}template<typename T,typename...Args>void wtb(const T&x,Args...args){wts(x);wtb(args...);}template<typename T,typename...Args>void wte(const T&x,Args...args){wts(x);wte(args...);}template<typename T1,typename T2>inline bool up(T1&x,const T2&y){return x<y?x=y,1:0;}template<typename T1,typename T2>inline bool dn(T1&x,const T2&y){return y<x?x=y,1:0;}template<typename T1,typename T2,typename T3>inline bool inRange(const T1&x,const T2&l,const T3&r){return!(x<l)&&!(r<x);}template<typename T1,typename T2>inline auto minOfDifferentTypes(const T1&x,const T2&y)->decltype(x<y?x:y){return x<y?x:y;}template<typename T1,typename T2>inline auto maxOfDifferentTypes(const T1&x,const T2&y)->decltype(x<y?y:x){return x<y?y:x;}template<typename T1,typename T2,typename T3>inline T1&madd(T1&x,const T2&y,const T3&modulo){return x=(ll)(x+y+modulo)%modulo;}template<typename T1,typename T2,typename T3>inline T1&mmul(T1&x,const T2&y,const T3&modulo){return x=(ll)x*y%modulo;}inline int modadd(int x,int y,int modulo){return(x+y)>=modulo?x+y-modulo:x+y;}inline int isinf(int x){return x<INF?x:-1;}inline void yesno(bool x){wtb(x?"Yes":"No");}
/* ------------------------------------------------------------------------------------------------------------------- */
typedef long double ld;
struct Matrix
{
ld a[3][3];
Matrix() { memset(a, 0, sizeof(a)); }
ld *operator[](int x) { return a[x]; }
const ld *operator[](int x) const { return a[x]; }
Matrix operator*(const Matrix &b) const
{
Matrix res;
For (i, 0, 2)
{
For (j, 0, 2)
{
For (k, 0, 2)
{
res[i][j] += a[i][k] * b[k][j];
}
}
}
return res;
}
Matrix operator*(ld x) const
{
Matrix res;
For (i, 0, 2)
{
For (j, 0, 2)
{
res[i][j] = a[i][j] * x;
}
}
return res;
}
Matrix operator+(const Matrix &b) const
{
Matrix res;
For (i, 0, 2)
{
For (j, 0, 2)
{
res[i][j] = a[i][j] + b[i][j];
}
}
return res;
}
};
signed main()
{
#ifdef FAST_IOSTREAM
ios_base::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
#endif
int n = read();
Matrix sq;
sq[0][0] = sq[1][0] = sq[1][1] = sq[2][0] = sq[2][2] = 1;
sq[2][1] = 2;
Matrix one;
one[0][0] = one[1][1] = one[2][2] = 1;
Matrix res = one;
For (i, 2, n) res = res * (sq + one * (i - 2)) * (1.0l / (i - 1)) + sq;
cout << fixed << setprecision(20) << res[2][0];
return 0;
}
推公式求答案的精确值
在 评论区 看到的：答案为 $6(n-1) - 4(1 + \frac{1}{2} + \frac{1}{3} + \cdots + \frac{1}{n - 1})$。
使用转移方程证明公式
评论区除了这个公式啥都没有..我就自己试着把它推出来了..
考虑增加一个点对答案的贡献，分为两部分：
- 新加的点度数为 $1$，贡献也为 $1$；
- 若父亲的度数为 $x$，则父亲的贡献增加 $2x+1$（因为 $(x+1)^2-x^2=2x+1$）。
所以，在一棵有 $i$ 个点的树上增加一个点的期望贡献为：（$f_{i, j}$ 为 上文中设计的状态）
$$ 1 + \sum\limits_{j=0}^i (2j+1) \frac{f_{i, j}}{i} $$
为了方便，在下文中设 $g(i)=\sum\limits_{j=0}^i (2j+1) f_{i, j}$。
这时再来观察一下转移方程：
- $f_{i-1, j}$，对 $g(i)$ 的贡献为 $g(i-1)$；
- $\frac{f_{i - 1, j - 1}}{i - 1} - \frac{f_{i - 1, j}}{i - 1}$，相当于把 $f_{i, j}$ 的 $\frac{1}{i-1}$ “移动”到了 $f_{i, j+1}$，从而使 $g$ 中的系数从 $2j+1$ 变成了 $2j+3$，即每一项的系数加二，所以对 $g(i)$ 的贡献为 $\frac{2}{i-1}\sum_{j=0}^{i-1}f_{i-1,j}=2$；
- $[j=1]$，对 $g(i)$ 的贡献为 $2j+1=3$。
所以，$g(i)=g(i-1)+5$。又因为由定义计算可得 $g(1)=1$，所以 $g(i)=5i-4$。
所以，在一棵有 $i$ 个点的树上增加一个点的期望贡献为 $1+\frac{g(i)}{i}=6-\frac{4}{i}$，累加一下就可以得到最上面那个公式了。
枚举两条边证明公式
这个做法是 @AprilGrimoire 告诉我的。
枚举两条边，对答案的贡献为这两条边期望相同的顶点个数（除了两条边是同一条边的情况，即两条边有公共顶点的概率）。
那么，答案为：
$$ 2(n-1) + 2\sum\limits_{i=3}^n\sum\limits_{j=2}^{i-1}\frac{2}{i-1} $$
其中，$2(n-1)$ 是枚举的两条边相同的贡献。
后面则是枚举每两条不同的边：$i$ 到其父亲和 $j$ 到其父亲。有公共顶点的概率是编号小的点（$j$）的父亲随便选，而编号大的点（$i$）的父亲需要是 $j$ 或 $j$ 的父亲，概率为 $\dfrac{2}{i-1}$。由于实际上需要枚举有序数对，再乘二。
化简一下就可以得到公式。
坑
不用 long double 可能会被卡精度..原因大概是题目要求的是 绝对 误差不超过 $10^{-6}$… 怎么还有人出题要求绝对精度误差的 /fad
只不过标算是求近似值，要是出题人会求精确值的做法就可以取模了 /fad
评论正在加载中...如果评论较长时间无法加载，你可以 搜索对应的 issue 或者 新建一个 issue 。