闲话
今日推歌：流星一過 / カラスヤサボウ 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';
}
}
《插值》
简单容斥得到答案
\[\sum_{i=0}^{n - 1} \left(\binom{n - i}{m}^2 - \binom{n - i - 1}{m}^2\right) F(i) \]然后我们发现这个求和是个差分形式，只需要计算这个玩意的前缀和形式即可
也就是
你发现这玩意上界是模数级别的，不能硬莽。考虑一手 \(\binom{n - i}{m}^2\) 是 \(2m\) 次多项式，\(F\) 是 \(m\) 次多项式，因此乘起来就是 \(3m\) 次，作前缀和就是 \(3m + 1\) 次多项式。
考虑插出来这个东西，因此我们就需要得到前 \(3m + 1\) 次项的系数，随后 \(O(m)\) 解决。
求这个东西的前 \(3m + 1\) 次项的话，首先构造 \(F\) 的前 \(3m + 1\) 次项就好办了。考虑拉格朗日插值，也就是
套路地设
\[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)\)。
当然跑得很慢。