2021-07-15 集训题解
光影交错
Description
Solution
可以看出，假设我们设置一个临界点 \(n\) ，当 \(n\) 足够大的时候，\(n\) 步操作之后对答案的影响就不再精度考虑范围之类了。
我们设 \(f(i)\) 分别表示 \(i\) 次操作时非中性灵气的期望出现次数，\(g(i)\) 表示 \(i\) 次非中性灵气是总体阳间的概率。
那么答案就是 \(\sum_{i=1}^{n} f(i)\times g(i)\) 。
可以得到：
\[g(x)=\sum_{i=\lfloor\frac{x}{2}\rfloor+1}^{x} \binom{x}{i}qL^iqR^{x-i}
\]
我们可以得到递推式：
\[\left\{\begin{array}{l}
g(x)=0,\text{when }x=0\\
g(x)=g(x-1)-pD\binom{x-1}{\lfloor\frac{x}{2}\rfloor} pL^{x/2}pD^{x/2-1},\text{when }x\equiv 0\pmod{2}\\
g(x)=g(x-1)+pL\binom{x-1}{\lfloor\frac{x}{2}\rfloor} pL^{(x-1)/2}pD^{(x-1)/2},\text{when }x\equiv 1\pmod{2}
\end{array}\right.
\]
考虑如何求解 \(f(x)\)，可以考虑对其构造生成函数 \(F(x)\) 。那么我们就有：
\[F(x)=\sum_{i=0}^{\infty} f(i)\times x^i=\sum_{i=1}^{\infty} (1-p)^{i-1}((pL+pD)x+(1-pL-pD))^i=\frac{1-p}{1-(1-p)[(pL+pD)x+(1-pL-pD)]}
\]
然后你就可以得到：
\[f(x)=\frac{f(x-1)(1-p)(pL+pD)}{1-(1-p)(1-pL-pD)}
\]
那么你就可以 \(\Theta(n)\) 解决这个问题，\(n\) 大概是 \(1e7\) 左右。
update
这里提供 神·tly 的一种精确求值 \(\Theta(1)\) 做法。
可以写出答案就是：
\[\sum_{l>r,e} \binom{l+r+e}{l+r}\binom{l+r}{l}pL^lpD^r(1-pL-pD)^e(1-p)^{l+r+e-1}
\]
考虑提出 \(e\)，那么答案就是：
\[\frac{1}{1-p}\sum_{l>r} \binom{l+r}{l}(pL\times (1-p))^l(pD\times (1-p))^r\sum_{e} \binom{l+r+e}{e}(1-pL-pD)^e(1-p)
\]
你发现 \(e\) 的贡献只与 \(l+r\) 有关，考虑它的生成函数，然后将 \(pL\times (1-p)\)写成 \(pL^{'}\)，\(pD\times (1-p)\) 写成 \(pD^{'}\)，发现答案可以写成：
\[\sum_{l>r} \binom{l+r}{l}(pL^{'})^l(pD^{'})^r\frac{1}{(1-pE)^{l+r+1}}
\]
其中 \(pE=(1-pL-pD)\times (1-p)\)，将 \(pL^{'}\times (1-pE)\) 写成 \(pL^{''}\)，同理替换成 \(pD^{''}\)，那么式子就是：
\[\frac{1}{(1-p)(1-pE)}\sum_{l>r} \binom{l+r}{l}(pL^{''})^l(pD^{''})^r
\]
考虑提出 \(pD^{''}\)，那么它的贡献就是：
\[S_D(n)=\sum_{i=0}^{n-1}(pD^{''})^i\binom{n+i}{n}
\]
\(n\) 也就是原式中的 \(l\)。
你发现通过错位相减可以得到：
\[S_D(n)(1-pD^{''})=S_D(n-1)+\binom{2(n-1)}{n-1}(pD^{''})^{n-1}-\binom{2n-1}{n}(pD^{''})^n
\]
为了方便，我们可以将 \(\binom{2(n-1)}{n-1}(pD^{''})^{n-1}-\binom{2n-1}{n}(pD^{''})^n\) 设为 \(a_n\)，那么 \(r\) 的贡献就可以表示为：
\[\sum_{i=1}^{l} \frac{a_i}{(1-pD^{''})^{l-i+1}}
\]
考虑交换求和符号，答案就是（忽略了前面的常数）：
\[\sum_{i=1}a_i\times (1-pD^{''})^{i-1}\sum_{l=i}^{\infty} (\frac{pL^{''}}{1-pD^{''}})^l
\]
将 \(\frac{pL^{''}}{1-pD^{''}}\) 设为 \(pH\)，你发现答案就是：
\[\sum_{i=1} a_i\times (1-pD^{''})^{i-1}\frac{pH^i}{1-pH}
\]
\[=\frac{1}{(1-pH)(1-pD^{''})}\sum_{i=1} a_i\times pL^{'}
\]
考虑将 \(a_i\) 展开，你发现就是（忽略系数）：
\[pL^{''}\sum_{i=0}^{\infty} \binom{2i}{i}(pD^{''})^i(pL^{''})^{i}-\sum_{i=1}^{\infty} \binom{2i-1}{i}(pD^{''})^i(pL^{''})^i
\]
对于前面一个你发现假设 Catalan 数列的生成函数为 \(C(x)\) ，那么前面一个的生成函数（设 \(x=pD^{''}pL^{''}\)）就是 \((C(x)\times x)^{'}\)，而我们知道 \(C(x)=C^2(x)\times x+1\)，所以前面一个的贡献就是 \(qL^{''}\frac{1}{\sqrt{1-4x}}\)，后面一个类似，可以推出就是 \(\frac{1}{2x}(\frac{1}{\sqrt{1-4x}}-1)\)。
至此，你就可以做到 \(\Theta(1)\) 计算了。
Code for O(n)
#include <bits/stdc++.h>
using namespace std;
#define Int register int
#define MAXN 10000005
template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c != ' ' && c != '\n') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
template <typename T> void chkmax (T &a,T b){a = max (a,b);}
template <typename T> void chkmin (T &a,T b){a = min (a,b);}
#define eps 1e-10
double p,pl,pd,f[MAXN],g[MAXN],drw[MAXN],tmp[MAXN];
void work (){
scanf ("%lf%lf%lf",&pl,&pd,&p);
double now = 1,ans = 0;int n = 0;tmp[0] = 1,tmp[1] = f[0] = f[1] = 0;
while (now > eps){
for (Int i = 1;~i;-- i){
tmp[i] *= (1 - pl - pd);
if (i) tmp[i] += tmp[i - 1] * (pl + pd);
f[i] += tmp[i] * now;
}
now *= (1 - p),n ++;
}
double h = 1 - (1 - p) * (1 - pl - pd);
for (Int i = 2;i <= n;++ i) f[i] = f[i - 1] * (1 - p) * (pl + pd) / h;
double sum = pl + pd;pl /= sum,pd /= sum;
for (Int i = 0;i <= n;++ i){
if (i & 1){
g[i] = g[i - 1] + drw[i - 1] * pl;
if (i == 1) drw[i] = pl;
else drw[i] = drw[i - 2] * pl * pd * i * (i - 1) / (i / 2) / (i / 2 + 1);
ans += g[i] * f[i];
}
else{
if (i == 0) drw[i] = 1;
else drw[i] = 4 * drw[i - 2] * (i - 1) / i * pl * pd;
if (i) g[i] = g[i - 1] - drw[i - 1] * pd;
ans += g[i] * f[i];
}
}
printf ("%.10f\n",ans);
}
signed main(){
freopen ("augury.in","r",stdin);
freopen ("augury.out","w",stdout);
int num,T;read (num,T);
while (T --> 0) work ();
return 0;
}
Code for O(1)
#include <bits/stdc++.h>
using namespace std;
#define Int register int
template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c != ' ' && c != '\n') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
template <typename T> void chkmax (T &a,T b){a = max (a,b);}
template <typename T> void chkmin (T &a,T b){a = min (a,b);}
double pL,pD,p;
void Work (){
scanf ("%lf%lf%lf",&pL,&pD,&p);
if (p == 1) return printf ("%.10f\n",pL),void ();
if (pL == 0) return printf ("%.10f\n",0.0),void();
double pE = (1 - pL - pD) * (1 - p),pL2 = pL * (1 - p) / (1 - pE),pD2 = pD * (1 - p) / (1 - pE);
double pH = pL2 / (1 - pD2),r = pD2 * pL2,u = 1.0 / sqrt(1 - 4 * r);
printf ("%.10f\n",pL2 / (1 - p) / (1 - pE) / (1 - pH) / (1 - pD2) * (!pD2 || !pL2 ? 1 : u - pD2 / (2 * r) * (u - 1)));
}
signed main(){
freopen ("augury.in","r",stdin);
freopen ("augury.out","w",stdout);
int num,T;read (num);read (T);
while (T --> 0) Work ();
return 0;
}
七星连珠
Description
Solution
md，处理 \(k\) 进制 \(\text{FWT}\) 根本不是难点好吧。。。
考虑 \(2\) 进制时如何处理，可以发现，如果我们如果对于每个位置 \((i,j)\)，构造一个数组，并将 \(a_{i,j}\) 赋值，再做 \(\text{DWT}\)，对整个矩阵求个行列式（可以看出这是一个 \(\text{DWT}\) 数组），拿行列式再 \(\text{IDWT}\) 回去，有值的位置就很有可能是正确答案。为什么说很有可能呢？因为行列式要乘 \((-1)^{\pi(p)}\)，就有可能把本来有值的消成 \(0\)，所以我们将 \(a_{i,j}\) 位置赋值的时候需要赋上一个随机值。
至于 \(3\) 进制，随便高维 \(\text{FWT}\) 搞搞就行了。因为要求单位根，所以找个模数 \(10^9+9\)，原根就是 \(13\)。
Code
#include <bits/stdc++.h>
using namespace std;
#define Int register int
#define mod 1000000009
#define MAXN 55
template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c == '-') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
template <typename T> void chkmax (T &a,T b){a = max (a,b);}
template <typename T> void chkmin (T &a,T b){a = min (a,b);}
int N,K,w[3][3],val[MAXN][MAXN];
int mul (int a,int b){return 1ll * a * b % mod;}
int dec (int a,int b){return a >= b ? a - b : a + mod - b;}
int add (int a,int b){return a + b >= mod ? a + b - mod : a + b;}
int qkpow (int a,int b){
int res = 1;for (;b;b >>= 1,a = mul (a,a)) if (b & 1) res = mul (res,a);
return res;
}
int inv (int x){return qkpow (x,mod - 2);}
void Add (int &a,int b){a = add (a,b);}
void Sub (int &a,int b){a = dec (a,b);}
#define MAXM 2187
int ivk,lim;
struct node{//储存一个fwt数组
int a[MAXM];
node(){memset (a,0,sizeof (a));}
int & operator [](const int key){return a[key];}
void fwt (int type){
for (Int i = 1;i < lim;i *= K)
for (Int j = 0;j < lim;j += i * K)
for (Int k = 0;k < i;++ k){
int tx[3] = {};
for (Int d = 0;d < K;++ d) for (Int e = 0;e < K;++ e) Add (tx[d],mul (type == 1 ? w[d][e] : inv(w[d][e]),a[j + k + i * e]));
for (Int d = 0;d < K;++ d) a[j + k + i * d] = mul (tx[d],type == 1 ? 1 : ivk);
}
}
}t,fuc[MAXN][MAXN];
int mat[MAXN][MAXN];
int getdet (){
int ans = 1;
for (Int i = 1;i <= N;++ i){
int pos = i;
for (Int j = i + 1;j <= N;++ j) if (mat[j][i]){pos = j;break;}
if (pos != i) swap (mat[i],mat[pos]),ans = mod - ans;
for (Int j = i + 1,iv = qkpow (mat[i][i],mod - 2);j <= N;++ j){
int del = mul (mat[j][i],iv);
for (Int k = i;k <= N;++ k) Sub (mat[j][k],mul (mat[i][k],del));
}
ans = mul (ans,mat[i][i]);
}
return ans;
}
signed main(){
freopen ("astrology.in","r",stdin);
freopen ("astrology.out","w",stdout);
srand (time(NULL));
int num;read (num,N,K),ivk = qkpow (K,mod - 2),lim = qkpow (K,7);
for (Int i = 1;i <= N;++ i)
for (Int j = 1;j <= N;++ j)
read (val[i][j]);
if (K == 2) w[0][0] = w[0][1] = w[1][0] = 1,w[1][1] = mod - 1;
else{
int gi = qkpow (13,(mod - 1) / 3);
for (Int i = 0;i < 3;++ i) w[i][0] = 1,w[i][1] = qkpow (gi,i),w[i][2] = mul (w[i][1],w[i][1]);
}
for (Int i = 1;i <= N;++ i)
for (Int j = 1;j <= N;++ j)
fuc[i][j][val[i][j]] = rand(),fuc[i][j].fwt(1);
for (Int x = 0;x < lim;++ x){
for (Int i = 1;i <= N;++ i)
for (Int j = 1;j <= N;++ j)
mat[i][j] = fuc[i][j][x];
t[x] = getdet();
}
t.fwt(-1);
for (Int i = 0;i < lim;++ i) if (t[i]) write (i),putchar (' ');
putchar ('\n');
return 0;
}
魔法咒语
Description
Solution
可以发现，答案就是：
\[\sum_{h=0}^{N-1} f(h)\sum_{n=h+1}^{N} (-1)^{n-h-1}\binom{n}{h}
\]
其中，\(f(h)=\sum_{j=0}^{k} h^j\)。
考虑后面那个东西如何递推。
因为存在 \(\binom{h}{i}=\binom{h-1}{i-1}+\binom{h-1}{i}\)，通过移项我们可以得到 \(\binom{h-1}{i-1}=\binom{h}{i}-\binom{h-1}{i}\Rightarrow \binom{h}{i}=\binom{h+1}{i+1}-\binom{h}{i+1}\) 。
将其代入原式中，我们可以得到：
\[g(i)=\sum_{h=i}^{N} (-1)^{h-i}\binom{h}{i}
\]
\[=\sum_{h=i}^{n} (-1)^{h-i}\binom{h+1}{i+1}+\sum_{h=i}^{n} (-1)^{h-i-1}\binom{h}{i+1}
\]
\[=2\times g_{i+1}+(-1)^{N-i}\binom{N+1}{i+1}
\]
然后你发现 \(i^{k+1}\) 用线性筛筛出来就好了，复杂度 \(\Theta(n)\)。
Code
#include <bits/stdc++.h>
using namespace std;
#define Int register int
#define int long long
#define MAXN 10000005
template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c != ' ' && c != '\n') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
template <typename T> void chkmax (T &a,T b){a = max (a,b);}
template <typename T> void chkmin (T &a,T b){a = min (a,b);}
bool vis[MAXN];
int N,k,tot,mod,g[MAXN],pw[MAXN],fac[MAXN],ifac[MAXN],prime[MAXN];
int mul (int a,int b){return 1ll * a * b % mod;}
int dec (int a,int b){return a >= b ? a - b : a + mod - b;}
int add (int a,int b){return a + b >= mod ? a + b - mod : a + b;}
int qkpow (int a,int b){
b %= (mod - 1),a %= mod;int res = 1;for (;b;b >>= 1,a = mul (a,a)) if (b & 1) res = mul(res,a);
return res;
}
void Add (int &a,int b){a = add (a,b);}
void Sub (int &a,int b){a = dec (a,b);}
void Euler (int up){
pw[1] = 1;
for (Int i = 2;i <= up;++ i){
if (!vis[i]) prime[++ tot] = i,pw[i] = qkpow (i,k + 1);
for (Int j = 1;j <= tot && i * prime[j] <= up;++ j){
vis[i * prime[j]] = 1,pw[i * prime[j]] = mul (pw[i],pw[prime[j]]);
if (i % prime[j] == 0) break;
}
}
}
int inv (int x){return mul (ifac[x],fac[x - 1]);}
int binom (int a,int b){return a >= b ? mul (fac[a],mul (ifac[b],ifac[a - b])) : 0;}
int f (int x){return x == 0 ? 1 : (x == 1 ? (k + 1) % mod : mul (pw[x] - 1,inv (x - 1)));}
signed main(){
freopen ("abracadabra.in","r",stdin);
freopen ("abracadabra.out","w",stdout);
int num;read (num,N,k,mod),Euler (N);
fac[0] = 1;for (Int i = 1;i <= N + 1;++ i) fac[i] = mul (fac[i - 1],i);
ifac[N + 1] = qkpow (fac[N + 1],mod - 2);for (Int i = N + 1;i;-- i) ifac[i - 1] = mul (ifac[i],i);
for (Int i = N;i >= 0;-- i) g[i] = add (add (g[i + 1],g[i + 1]),mul (N - i & 1 ? mod - 1 : 1,binom (N + 1,i + 1)));
int ans = 0;for (Int h = 0;h <= N;++ h) Add (ans,f(h)),Sub (ans,mul (f(h),g[h]));
write (ans),putchar ('\n');
return 0;
}
