闲话 23.1.17

闲话

数学 \(3\)

怎么说？算了不感谢了

题解.txt：skyh's blog
彳亍

解方程

上来以为是分拆数假了一阵子。后来发现这个东西如果没有小于等于 \(a_i\) 的限制后就是插板法算组合数，因此考虑往那个方向转化。

首先为了凑出插板法需要让所有 \(X_i\) 取值从 \(0\) 开始，因此首先让 \(m\) 减去 \(n\)。\(n_1+1\sim n_2\) 没啥用，直接让 \(m\) 减去 \(a_i\) 后这部分的取值是任意的。注意先前已经减去一个 \(1\) 了。

随后考虑容斥。我们直接做子集反演，设集合 \(S\) 对应的状态 \(f(S)\) 表示钦定 \(|S|\) 个值不满足条件，剩下的值随便的计数。答案就是 \(\sum (-1)^{|T|} f(T)\)。
\(f(T)\) 是易求的，我们只需要将不满足条件的值如上地处理掉，剩下的就是随意划分了。假设这种情况下 \(m\) 剩余了 \(sum\)，答案就是 \(\dbinom{sum + n - 1}{n - 1}\)。

由于模数的质因子都是小质数，因此可以采用 exLucas 解决。

code

#include <bits/stdc++.h>
using namespace std; 
#define int long long
using pii = pair<int, int>; using vi = vector<int>; using vp = vector<pii>; using ll = long long; 
using ull = unsigned long long; using db = double; using ld = long double; using lll = __int128_t;
#define multi int T; cin >> T; while ( T -- )
#define timer cerr << 1. * clock() / CLOCKS_PER_SEC << '\n';
#define iot ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define rep(i, s, t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i, s, t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
#define eb emplace_back
#define pb pop_back
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;
int p, n, n1, n2, m, a[50], ans;

int exgcd(int a, int b, int &x, int &y) {
    if (b == 0)  {
    	x = 1, y = 0;
    	return a;
	} int d = exgcd(b, a % b, x, y);
    int tmp = x;
	x = y; y = tmp - a/b*y; 
	return d;
}
inline int lcm(int a, int b) {
    return a / __gcd(a, b) * b;
}
inline int inv(int a, int p) {
    int x, y;
    exgcd(a, p, x, y);
    return (x + p) % p;
}
inline int qp(int a, int b, int p) {
    int ret = 1; a %= p;
    while (b) {
        if (b & 1) ret = ret * a % p;
        b >>= 1;
		a = a * a % p;
    } return ret;
}

inline int F(int n, int P, int PK) {
    if (n == 0) return 1;
    int rou = 1;
    int rem = 1; 
    rep(i, 1, PK) if (i % P != 0) rou = rou * i % PK;
    rou = qp(rou, n/PK, PK);
    rep(i, PK * (n / PK), n) if (i % P != 0) rem = rem * (i % PK) % PK;
    return F(n / P, P, PK) * rou % PK * rem % PK;
}
inline int G(int n, int P) {
    if (n < P) return 0;
    return G(n/P, P) + (n/P);
}

inline int C_PK(int n, int m, int P, int PK) {
    int fz  = F(n, P, PK), 
		fm1 = inv(F(m, P, PK), PK), 
		fm2 = inv(F(n - m, P, PK), PK);
    int mi = qp(P, G(n, P) - G(m, P) - G(n - m, P), PK);
    return 1ll * fz * fm1 % PK * fm2 % PK * mi % PK;
} 

int A[1001], B[1001], tot;
inline int C(int n, int m, int P) {
    int tmp = P; tot = 0;
    for (int i=2; i * i <= P; i++) {
        if (tmp % i == 0) {
            int PK=1;
            while (tmp % i == 0) {
                PK *= i;
				tmp /= i;
            } A[++tot] = PK;
			B[tot] = C_PK(n, m, i, PK);
        }
    } if (tmp != 1) {
        A[++tot] = tmp;
		B[tot] = C_PK(n, m, tmp, tmp);
    } int ans = 0;
    rep(i, 1, tot) {
        int M = P / A[i], T = inv(M, A[i]);
        ans = (ans + 1ll * B[i] * M % P * T % P) % P;
    } return ans;
}

signed main() {
	int T; cin >> T >> p;
	while (T --) {
		ans = 0;
		cin >> n >> n1 >> n2 >> m; m -= n;
		rep(i, 1, n1+n2) cin >> a[i];
		rep(i, n1+1, n1 + n2) m -= a[i] - 1;
		if (n == 0 or m < 0) { puts("0"); continue; }
		rep(i, 0, (1<<n1)-1) {
			int now = m;
			rep(j, 0, n1-1) if ((i >> j) & 1) now -= a[j + 1];
			if (now < 0) continue;
			if (__builtin_popcount(i) & 1) 
				 ans = (ans - C(n + now - 1, n - 1, p) + p) % p;
			else ans = (ans + C(n + now - 1, n - 1, p)) % p;
		} cout << ans << '\n';
	} timer;
}

宇宙序列

刚学完 FWT 就考欸好耶

容易发现 \(a_i\) 序列是 \(a_1\) 序列自卷了 \(i\) 次。
由于异或卷积有结合律，我们可以直接作倍增，这样的时间复杂度是 \(O(pn2^n)\) 的。

由于 FWT 是线性变换，倍增时可以一直通过累加 \(\text{FWT}(a_1)\) 数组来得到答案。我们首先将 \(a_1\) 序列做 FWT，然后每位不断自乘后累加在 \(b\) 序列的对应位上，最后再将 \(b\) 序列做 IFWT 就得到了答案。
这样就优化到了 \(O(n2^n + q2^n)\)。

观察上一个做法，可以发现复杂度出现在不断自乘上。换句话说，我们需要优化的就是

\[b_k = \sum_{i=0}^{p} a_k^{2^{i}} \]

~~因为我想到这时没考虑 p 不是 2 的次方形式，所以~~考虑分治，我们要求的是 \(\sum_{i=0}^{2^k - 1} a^{2^i}\)。
我们做类似矩阵快速幂的做法，将前一半和后一半分开，也就得到了

\[\begin{aligned} &\sum_{i=0}^{2^k - 1} a^{2^i} \\ = \ & \sum_{i=0}^{2^{k-1} - 1} a^{2^i} + \sum_{i=2^{k-1}}^{2^{k} - 1} a^{2^i} \\ = \ & \sum_{i=0}^{2^{k-1} - 1} a^{2^i} + \sum_{i=0}^{2^{k} - 2^{k-1} 1} a^{2^{i+2^{k-1}}} \\ = \ & \sum_{i=0}^{2^{k-1} - 1} a^{2^i} + \sum_{i=0}^{2^{k-1} 1} a^{2^{i} \times 2^{2^{k-1}}} \\ = \ & \sum_{i=0}^{2^{k-1} - 1} a^{2^i} + \sum_{i=0}^{2^{k-1} 1} \left( a^{2^{2^{k-1}}} \right) ^{2^i} \end{aligned}\]

可以发现这个东西前后是同构的。我们设 \(f(a, k) = \sum_{i=0}^{2^k - 1} a^{2^i}\)，则贡献有很奇妙的转移：

\[f(a, k) = f(a, k - 1) + f(a^{2^{2^{k-1}}}, k - 1) \]

由于模数很小，所以这个东西可以 \(O(\text{mod}\log p)\) 地统计得到。

到这我突然发现 \(p\) 不是 \(2^k\)。~~然后我摆了~~
然后其实对于任意 \(p\) 我们可以通过按位处理来得到值，复杂度为 \(O(\log p)\)。
得到 \(O(2^n)\) 个值后我们就可以 IFWT 回去了。

总时间复杂度 \(O(n2^n + \text{mod}\log p + 2^n \log p)\)。
代码挺短的，但是比较难调（~~谁让你摆了一上午最后半小时调啊~~

预处理部分似乎可以 \(O(\text{mod}^2)\) 地找循环节得到。看不懂，有想了解的可以问 lyin。

code

#include <bits/stdc++.h>
using namespace std; 
using pii = pair<int,int>; using vi = vector<int>; using vp = vector<pii>; using ll = long long; 
using ull = unsigned long long; using db = double; using ld = long double; using lll = __int128_t;
#define multi int T; cin >> T; while ( T -- )
#define timer cerr << 1. * clock() / CLOCKS_PER_SEC << '\n';
#define iot ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
#define eb emplace_back
#define pb pop_back
const int N = 1e6 + 10, mod = 10007, phi_1 = 5002;
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;
using poly = vector<int>;
int n, len, p, q, ans, f[10007][32];
poly A, B;

int qp(int a, int b, int p = mod) {
	int ret = 1;
	while (b > 0) {
		if (b & 1) ret = 1ll * ret * a % p;
		a = 1ll * a * a % p;
		b >>= 1;
	} return ret;
}

void FWT(poly& A, const int& len, int typ = 1) {
    if (typ != 1) typ = (mod + 1) >> 1;
    for (int mid = 1, t; mid < len; mid <<= 1) 
        for (int i = 1, j = 0; j < len; j += (mid << 1)) 
            for (int k = 0; k < mid; ++ k) 
                t = A[j + k + mid],
                A[j + k + mid] = 1ll * (A[j + k] - t + mod) * typ % mod,
                A[j + k] = 1ll * (A[j + k] + t) * typ % mod;
}

signed main() {
	iot; cin >> n >> p >> q; len = 1 << n;
    A.resize(len), B.resize(len);
    rep(i,0,len-1) cin >> A[i]; 
	FWT(A, len, 1);
	rep(i,0,mod-1) f[i][0] = i;
	rep(j,1,30) rep(i,0,mod-1) 
		f[i][j] = (f[i][j - 1] + 
			f[qp(i, qp(2, (1 << j - 1) % phi_1, mod - 1), mod)][j - 1]) % mod;
	rep(i,0,len-1) {
		int tmp = 0;
		rep(j,0,30) if ((p + 1) >> j & 1) {
			B[i] = (B[i] + f[qp(A[i] , qp( 2 , tmp , mod - 1 ) , mod)][j]) % mod;
			tmp = (tmp + (1 << j)) % phi_1;
		}
	} 
	FWT(B, len, 0);
	cout << B[q] << '\n';
	timer;
}

exp

奇妙概率 dp。

我们发现每次寻找总会在一个 \('\texttt{.}'\) 处停止，因此我们可以将一段连续子串划分为 \('\texttt{.}'\) 两侧的两个子串，它们之间的答案是不会互相影响的。然后我们可以按照结尾为 \('\texttt{.}'\) 的区间进行 DP。环不太好 DP，我们首先断环成链复制一倍，最后枚举断点即可。容易发现断点就是最后一个被消掉的 \('\texttt{.}'\)。

设 \(f([l, r])\) 表示 \([l,r)\) 已经消完、最后一个 \('\texttt{.}'\) 在 \(r\) 位置对应的期望。并设 \(g([l,r])\) 表示对应的概率。
设 \(p([l,r],k)\) 表示区间 \([l,k)\) 已经消完，区间 \((k, r)\) 已经消完，\(k\) 是最后一个被消的 \('\texttt{.}'\)，\(r\) 未消的概率。

然后设 \([l, k)\) 中有 \(L\) 个 \('\texttt{.}'\)，\((k, r)\) 中有 \(R\) 个。然后可以写出

\[p([l,r],k) = \binom{L + R}{L}\times \left(\frac{k - l + 1}{r - l + 1}\right)^{L} \times g([l, k]) \times \left(\frac{r - k}{r - l + 1}\right)^{R} \times g([k + 1, r]) \times \frac{k - l + 1}{r - l + 1} \]

也就是选择左右两侧的次数合法的方案数除以总方案数，再乘以到达这个子问题的概率。

前半部分可以首先在左侧选 \(L\) 次，再在右侧选 \(R\) 次，然后转自由选择的合法方案数就是 \((k - l + 1)^{L} \times g([l, k]) \times (r - k)^{R} \times g([k + 1, r]) \times \dbinom{L + R}{L}\)。由于我们需要让最后一次合法，也就是最后一次从左侧开始选，选中 \(k\) 位置的 \('\texttt.'\)，可能性再乘以 \((k - l + 1)\)。
总可能性就是 \((r - l + 1)^{L + R + 1}\)。

得到 \(p\) 的贡献式后我们就能得到 \(f, g\) 的式子了：

\[g([l,r]) = \sum_{k = l}^{r-1} p([l,r],k) \]

\[f([l,r]) = \frac{1}{g([l,r])} \sum_{k = l}^{r - 1} p([l,r], k)\times \left(f([l,k]) + f\left([k+1, r]\right) + \frac{k - l}{2}\right) \]

\(f\) 是对应子区间贡献的加权平均值，最后一次从左边选的期望就是 \(\sum_{i=0}^{k - l + 1}i / (k - l+ 1) = (k - l) / 2\)。

初始化直接找到所有只有最右侧一个 \('\texttt{.}'\) 的区间 \([l,r]\)，设它的 \(f = 0, g = 1\) 即可。

枚举最后一个贡献的 \('\texttt{.}'\) 即可，它的贡献是 \(\dfrac{n - 1}{2}\)。假设我们找到一个为 \(.\) 的 \(i> n\)，它对答案的贡献就是 \(f([i - n + 1, i])\times g([i - n + 1, i])\)。求和即可。
由于 \(\sum g = 1\)，取平均值不需要除。

上面是抄 skyh 博客的。

code

#include <bits/stdc++.h>
using namespace std; 
using pii = pair<int,int>; using vi = vector<int>; using vp = vector<pii>; using ll = long long; 
using ull = unsigned long long; using db = double; using ld = long double; using lll = __int128_t;
#define multi int T; cin >> T; while ( T -- )
#define timer cerr << 1. * clock() / CLOCKS_PER_SEC << '\n';
#define iot ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
#define eb emplace_back
#define pb pop_back
const int N = 200 + 10;
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;
int n, u, sum[N << 1];
db ans, f[N << 1][N << 1], g[N << 1][N << 1], C[205][205];
char ch[N << 1];

db qp(db a, int b) {
    db ret = 1;
    while (b) {
        if (b & 1) ret = ret * a;
        a = a * a;
        b >>= 1;
    } return ret;
}

db p(int l, int r, int k) {
    int L = sum[k - 1] - sum[l - 1], R = sum[r - 1] - sum[k];
    return C[L + R][L] * qp(1. * (k - l + 1) / (r - l + 1), L) * g[l][k] * qp(1. * (r - k) / (r - l + 1), R) * g[k + 1][r] * (1. * (k - l + 1) / (r - l + 1));
}

signed main() {
    iot; cout << setprecision(11) << fixed;
    rep(i,0,200) C[i][0] = C[i][i] = 1;
    rep(i,1,200) rep(j,1,i - 1) C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
    multi {
        cin >> ch + 1; n = strlen(ch + 1), ans = 0;
        rep(i,1,n) sum[i + n] = sum[i] = (ch[i] == '.');
        rep(i,1,n<<1) sum[i] += sum[i - 1];
        rep(i,1,n<<1) rep(j,i,i+n-1) 
            if (sum[j] - sum[i - 1] == 1 and sum[j] - sum[j - 1] == 1) f[i][j] = 0, g[i][j] = 1;
            else f[i][j] = 0, g[i][j] = 0;
        rep(len,1,n) for (int l = 1, r = l + len - 1; r <= (n << 1); l ++, r ++) {
            if (sum[r] - sum[r - 1] == 0 or g[l][r]) continue;
            rep(k,l,r-1) if (sum[k] - sum[k - 1] != 0) {
                db pnow = p(l, r, k);
                g[l][r] += pnow;
                f[l][r] += pnow * (f[l][k] + f[k + 1][r] + 1. * (k - l) / 2);
            } if (g[l][r]) f[l][r] /= g[l][r];
        } rep(i,n,(n<<1)-1) if (sum[i] - sum[i - 1] == 1) 
            ans += (f[i - n + 1][i] + 1. * (n - 1) / 2) * g[i - n + 1][i];
        cout << ans << '\n';
    }
}

posted @ 2023-01-17 13:48 joke3579 阅读(93) 评论(0) 收藏举报

刷新页面返回顶部