闲话 23.1.15

闲话

今日推歌：流星一過 / カラスヤサボウ feat. 镜音铃&镜音连

计算几何怎么学？
也不考啊……

没东西好写了，就给大家拜个早年吧！

数学 \(2\)

同样首先感谢大自然的馈赠。

猜拳游戏

修改游戏规则，把平局情况改成一直玩到赢，则己方现在的胜利概率为 \(\frac{win}{win + lose}\)。
然后求解这个东西的最大值可以做 01 分数规划。二分一个 \(\text{mid}\)，然后 dp 出 \(win - \text{mid} * lose\) 的最大值即可。这个意义就是赢了得 \(1\) 分，输了扣 \(k\) 分，看己方是否得分 \(\ge 0\)。
可以倒序 dp，dp 状态可以设 \(f(i, j)\) 为玩到第 \(i\) 轮，之前差距为 \(j\) 场时的最大值。可以第二维 \(+n\) 再做。

得到最大值后做高斯消元即可。方程就是 \(f_i = win\times f_{i + 1} + lose \times f_{i - 1}\)。

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 = 1e3 + 10;
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;

int n, m1, m2, r[N], p[N], s[N];
db ans, f[N][N << 1];
bool check(db mid) {
	rep(i,0,n-1) f[n + 1][i] = -mid, f[n + 1][n] = 0;
	rep(i,n+1,n<<1) f[n + 1][i] = 1;
	pre(i,n,1) rep(j,0,n<<1) {
		db t1 = 0, t2 = 0, t3 = 0;

		if (j < (n << 1)) t1 += 1. * s[i] / 100 * f[i + 1][j + 1];
		t1 += 1. * r[i] / 100 * f[i + 1][j];
		if (j) t1 += 1. * p[i] / 100 * f[i + 1][j - 1];


		if (j < (n << 1)) t2 += 1. * p[i] / 100 * f[i + 1][j + 1];
		t2 += 1. * s[i] / 100 * f[i + 1][j];
		if (j) t2 += 1. * r[i] / 100 * f[i + 1][j - 1];

		if (j < (n << 1)) t3 += 1. * r[i] / 100 * f[i + 1][j + 1];
		t3 += 1. * p[i] / 100 * f[i + 1][j];
		if (j) t3 += 1. * s[i] / 100 * f[i + 1][j - 1];

		f[i][j] = max( { t1, t2, t3 } );
	} return f[1][n] >= 0;
}

db mat[N][N], x[N];
db gauss(int n) {
	for (int r = 1, c = 1; r <= n and c <= n; ++ r, ++ c) {
		int mx = r;
		rep(i,r+1, n) if (fabs(mat[i][c]) > fabs(mat[mx][c])) mx = i;
		if (mx != r) rep(i,c,n+1) swap(mat[mx][i], mat[r][i]);
		if (!mat[r][c]) continue;
		rep(i,r+1,n) {
			db tmp = mat[i][c] / mat[r][c];
			rep(j,c,n+1) mat[i][j] -= tmp * mat[r][j];
		}
	} rep(i,1,n) x[i] = 0;
	pre(i,n,1) {
		db tmp = mat[i][n + 1];
		rep(j,i+1,n) tmp -= mat[i][j] * x[j];
		x[i] = tmp / mat[i][i];
	} return x[m2];
}

signed main() {
	while (1) {
		cin >> n >> m1 >> m2;
		if (n == 0) break;
		rep(i,1,n) cin >> r[i] >> p[i] >> s[i];
		db l = 0, r = 1e6, mid;
		rep(i,1,45) {
			mid = (l + r) / 2;
			if (check(mid)) l = mid;
			else r = mid;
		} db q = 1. / (1. / l + 1.);
		// memset(mat, 0, sizeof mat);
		rep(i,2,m1+m2-2) mat[i][i] = 1, mat[i][i - 1] = q - 1, mat[i][i + 1] = -q;
		mat[m1 + m2 - 1][m1 + m2 - 1] = 1, mat[m1 + m2 - 1][m1 + m2] = q;
		if (m1 + m2 - 1 > 1) mat[1][2] = -q, mat[1][1] = 1, mat[m1 + m2 - 1][m1 + m2 - 2] = q - 1;
		cout << setprecision(5) << fixed << gauss(m1 + m2 - 1) << '\n';
	}
}

---

B君的回忆

奇妙题。

容易发现是迭代 \(k\) 次，每次得到的答案会被作为下一次的下标来计算。考虑如何将下标的范围缩小到一个合理的范围内。
发现这个东西在 \(\bmod\ p\) 下是存在循环的，具体可以通过鸽巢原理以及和斐波那契数列相似的转移发现。
这个东西可以通过 bsgs 来求解，复杂度是根号级别的。

复杂度主要出现在一个大数的根号上，思考如何将大数减小。
发现可以首先将合数质因子分解，然后其循环节长度就是其质因子的循环节的 \(\text{lcm}\)，证明显然。
具体关于循环节长度的上界可以看 yspm 曾阅读过的写 斐波那契模意义下循环节的文章。
可以发现其循环节上界是 \(6\bmod\) 级别的。如果设 \(g(0) = 0, g(1) = 1\) 可以发现 \(g(n)\) 是斐波那契数列的第 \(n\) 项，大概可以得到更紧的上界。

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 N = 1e6 + 10;
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;
int a, b, n, k, p;
unordered_map<ll,int> hsh, cyc;

int phi(int k) {
	int ret = k;
	for (int i = 2; i * i <= k; ++ i) if (k % i == 0) {
		ret = ret / i * (i - 1); while (k % i == 0) k /= i;
	} if (k > 1) ret = ret / k * (k - 1);
	return ret;
}

struct mat {
	int a[2][2];
	mat() { memset(a, 0, sizeof a); }
	int* operator[] (const int& p) { return a[p]; }
	const int* operator[] (const int& p) const { return a[p]; }
	void clear() { memset(a, 0, sizeof a); }
	void I() { a[0][0] = a[1][1] = 1, a[0][1] = a[1][0] = 0; }  
	bool operator == (const mat& b) {
		return a[0][0] == b[0][0] and a[0][1] == b[0][1] and a[1][0] == b[1][0] and a[1][1] == b[1][1];
	}
	ll gethsh() {
		ll ret = 0;
		rep(i,0,1) rep(j,0,1) ret = ret * 131 + a[i][j];
		return ret;
	}
}unit;

void mul(mat& a, mat& b, int mod) {
	mat ret;
	rep(i,0,1) rep(j,0,1) rep(k,0,1) ret[i][j] = (ret[i][j] + 1ll * a[i][k] * b[k][j]);
	rep(i,0,1) rep(j,0,1) a[i][j] = ret[i][j] % mod;
}

int get(int n, int pnow) {
	mat bse; 
	bse[0][0] = 0, bse[0][1] = pnow - 1, bse[1][0] = 1, bse[1][1] = 3;
	int v[2] = { a % pnow, b % pnow };
	ll tmp[2];
	while (n) {
		if (n & 1) {
			tmp[0] = tmp[1] = 0;
			rep(i,0,1) rep(j,0,1) tmp[j] += 1ll * v[i] * bse[i][j];
			rep(i,0,1) v[i] = tmp[i] % pnow;
		} mul(bse, bse, pnow);
		n >>= 1;
	} return v[0];
} 

int find(int x) {
	if (x == 1) return 1;
	if (cyc.count(x)) return cyc[x];
	int up = x * 4, sq = sqrt(up);
	hsh.clear();
	mat bse, now; now.I();
	bse[0][0] = 0, bse[0][1] = x - 1, bse[1][0] = 1, bse[1][1] = 3;
	rep(i,0,sq-1) {
		hsh[now.gethsh()] = i;
		mul(now, bse, x);
		if (now == unit) return cyc[x] = i + 1;
	}
	bse = now, mul(now, now, x);
	rep(i,2,up/sq) {
		if (hsh.count(now.gethsh())) return cyc[x] = i * sq - hsh[now.gethsh()];
		mul(now, bse, x);
	} assert(0);
}

int div(int p) {
	if (p == 1) return 1;
	int ret = -1;
	for (int i = 2; i * i <= p; ++ i) if (p % i == 0) {
		int cnt = 1, tmp;
		while (p % i == 0) cnt *= i, p /= i;
		tmp = find(i) * cnt / i;
		if (~ret) ret = ret / __gcd(ret, tmp) * tmp;
		else ret = tmp;
	} if (p > 1) {
		int tmp = find(p);
		if (~ret) ret = ret / __gcd(ret, tmp) * tmp;
		else ret = tmp; 
	} return ret;
}

int f(int n, int k, int pnow) {
	if (k == 1) return get(n, pnow);
	return get(f(n, k - 1, div(pnow)), pnow);
}

signed main() {
	unit.I();
	multi {
		cin >> a >> b >> n >> k >> p;
		cout << f(n, k, p) << '\n';
	}
} 

---

loj2504

《插值》

简单容斥得到答案

\[\sum_{i=0}^{n - 1} \left(\binom{n - i}{m}^2 - \binom{n - i - 1}{m}^2\right) F(i) \]

然后我们发现这个求和是个差分形式，只需要计算这个玩意的前缀和形式即可
也就是

\[\sum_{i=0}^{n - 1} \binom{n - i}{m}^2 F(i) \]

你发现这玩意上界是模数级别的，不能硬莽。考虑一手 \(\binom{n - i}{m}^2\) 是 \(2m\) 次多项式，\(F\) 是 \(m\) 次多项式，因此乘起来就是 \(3m\) 次，作前缀和就是 \(3m + 1\) 次多项式。
考虑插出来这个东西，因此我们就需要得到前 \(3m + 1\) 次项的系数，随后 \(O(m)\) 解决。
求这个东西的前 \(3m + 1\) 次项的话，首先构造 \(F\) 的前 \(3m + 1\) 次项就好办了。考虑拉格朗日插值，也就是

\[F(n) = \sum_{i = 0}^m F(i) \prod_{j \neq i}^m \frac{n - j}{i - j} \]

套路地设

\[G(i) = F(i) \prod_{j \neq i}^m \frac{1}{i - j} \]

于是能得到

\[\begin{aligned} F(n) &= \sum_{i = 0}^m G(i) \prod_{j \neq i}^m (n - j) \\ &= \sum_{i = 0}^m \frac{1}{i - j} \prod_{j = 1}^m G(i) (n - j) \\ &= \prod_{j = 1}^m (n - j)\sum_{i = 0}^m \frac{G(i) }{i - j} \\ &= \left(\prod_{j = 1}^m (n - j)\right) \left(\sum_{i = 0}^m \frac{G(i) }{i - j} \right) \end{aligned}\]

\(G\) 可以通过逆元前缀积以及讨论 \(-1\) 奇偶性得到，随后只需卷上 \(\dfrac{1}{i}\) 即可。最后可以轻易处理 \(F\)。
随后只需要动态维护 \(\dbinom{i}{m}\) 得到辅助多项式，再卷上 \(F\) 即可得到所需要的前缀和。

最后做两次拉格朗日插值即可。总时间复杂度为 \(O(n\log n)\)。
当然跑得很慢。

Submission.