化學磁牛
https://tioj.ck.tp.edu.tw/problems/1829
Description
現在數線上 $K$ 的地方有隻磁牛，磁牛每秒有 $p = \frac{A}{B}$ 往左移動，$1-p$ 的距離往右移動。當移動到 $0$ 或移動到 $N$ 就會停止移動，問最後停在 $N$ 的機率是多少？
假設答案的最簡分數是 $\frac{P}{Q}$，則你應該輸出 $P$ 除以 $M$ 的餘數和 $Q$ 除以 $M$ 的餘數共兩個數字，其中 $M = 10 ^ 8 + 7$。
有 $T$ 筆測資，$T\leq 10^ 4; K<N\leq 10^ {19}; A, B\leq 3\times 10^ 9;\gcd(A,B)=1$
Solution
本題可以拆成兩個部份。第一部份是推出機率的式子，第二部份是按照式子求出最簡分數的分子分母。
雖然我原本是用 martingale 推的，不過這裡就放一個比較直覺的版本。令 $f(i)$ 表示在位置 $i$ 開始走的時候，最終停在 $N$ 的機率。那麼有
$$
\begin{cases}
f(0) &= 0 \\
f(N) &= 1 \\
f(i) &= pf(i - 1) + (1 - p) f(i + 1), \forall {0 < i < N}
\end{cases}
$$
將最後一行變形可得
$$
p \cdot \left( f(i) - f(i - 1) \right) = (1 - p)\cdot \left( f(i + 1) - f(i) \right)
$$
不妨令 $r = p/(1 - p)$（$p=0,1$ 的情況都特判），$f(1) - f(0) = c$，那麼可以得到
$$
f(i) = c + rc + \dots + r^{i-1}c = \frac{r^{i}-1}{r-1} \cdot c
$$
再由 $f(N)=1$ 可以得到 $f(K) = \frac{r^K-1}{r^N-1}$。為了使用等比級數公式我們需要 $r\neq 1$，因此 $r=1$ 的情形就另外再推一下（此時 $p=1/2$，$f(K) = K / N$）。
接著是第二部份求出最簡分數。利用 $r = (A/B) / (1 - A/B) = A / (B - A)$ 把 $\frac{r^K-1}{r^N-1}$ 化簡一下。令 $C = B - A$。
$$
\frac{r^K-1}{r^N-1}
= \frac{(A / C) ^ K - 1}{(A / C) ^ N - 1}
= \frac{A^K C^{N-K} - C^N}{A^N - C^N}
= \frac{C^{N-K}(A^K - C^K)}{A^N - C^N}
$$
此時上下都是整數了，但有可能上下不互質而不是最簡分數。先假設 $A > C$，這樣上下至少都是正的。因為 $A \perp B$，我們可以推出 $A \perp C$ 以及 $C^{N-K} \perp (A^N-C^N)$。
所以分子跟分母的 gcd 就是 $\gcd(A^K-C^K, A^N-C^N)$，但 $A, C$ 又互質，我們可以知道 $\gcd(A^K-C^K, A^N-C^N) = A^{\gcd(K,N)} - C^{\gcd(K,N)}$。因此精確的最簡分數就是
$$
\frac{C^{N-K}(A^K - C^K) / (A^g - C^g)}{(A^N - C^N) / (A^g - C^g)}
$$
其中 $g = \gcd(K,N)$。這個式子在 $A < C$ 的時候因為正負號會抵銷也剛好是對的。
edit: 原本我以為只要輸出分子是 $C^{N-K}(A^K - C^K) (A^g - C^g) ^ {-1} \pmod M$、分母是 $(A^N - C^N)(A^g - C^g) ^ {-1} \pmod M$ 就好了，但在寫完題解之後發現 $A^g - C^g$ 可能模 $M$ 之後是 $0$，這樣除以 $0$ 會壞掉（也就是說本題測資目前沒有卡這個 case）。不過，令 $K = qg$ 的話可以如下化簡
$$
\frac{A^K - C^K}{A^g - C^g} = \frac{(A^g)^{q} - (A^g)^{q}}{A^g - C^g}
= (A^g)^0 (B^g)^{q-1} + (A^g)^1 (B^g)^{q-2} + \cdots + (A^g)^{q-1}(B^g)^0
$$
然後 $u^0 v^{q-1} + u^1 v^{q-2} + \cdots + u^{q-1} v^0$ 可以用分治求得。
有關 $\gcd(a^m - b^m, a^n - b^m) = a^g - b^g$ 的證明
其中 $\gcd(a, b) = 1, \gcd(m, n) = g$。
從大家學過的 $(a - b) | (a^k - b^k)$ 可以知道 gcd 至少會是 $a^g - b^g$ 的倍數，但恰好是 $a^g - b^g$ 這件事在 AC 之前我是亂矇的。
後來在 stackexchange 查到證明，我理解如下：
假設 $n = qm + r, 0 \leq r < m$。
$$
\begin{align*}
a^n - b^n
&= (a^{qm} - b^{qm}) \cdot a^r + b^{qm} (a^r - b^r) \\
&= (a^m - b^m) \cdot ((a^m) ^ {0} (b^m)^{q-1} + (a^m) ^ {1} (b^m)^{q-2} + \cdots + (a^m) ^ {q-1} (b^m)^{0}) \cdot a^r +
b^{qm}(a^{r} - b^{r})
\end{align*}
$$
又 $b^k$ 跟 $a^n-b^n$ 一定是互質的，因此 $\gcd(a^n - b^n, a^m - b^m) = \gcd(a^m - b^m, a^{r} - b^{r})$，按照輾轉相除法這樣下去就最終會抵達 $a^g - b^g$。
AC code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
| #include <bits/stdc++.h>
using namespace std;
#define all(x) begin(x), end(x)
#ifdef local
#define safe cerr << __LINE__ << " line " << __LINE__ << " safe\n"
#define debug(a...) debug_(#a, a)
#define orange(a...) orange_(#a, a)
template <typename ...T>
void debug_(const char *s, T ...a) {
cerr << "\e[1;32m(" << s << ") = (";
int cnt = sizeof...(T);
(..., (cerr << a << (--cnt ? ", " : ")\e[0m\n")));
}
template <typename I>
void orange_(const char *s, I L, I R) {
cerr << "\e[1;32m[ " << s << " ] = [ ";
for (int f = 0; L != R; ++L)
cerr << (f++ ? ", " : "") << *L;
cerr << " ]\e[0m\n";
}
#else
#define safe ((void)0)
#define debug(...) safe
#define orange(...) safe
#endif
using lld = int64_t;
template <typename T, T MOD> class Modular {
public:
constexpr Modular() : v() {}
template <typename U> Modular(const U &u) { v = static_cast<T>(0 <= u && u < MOD ? u : (u%MOD+MOD)%MOD); }
template <typename U> explicit operator U() const { return U(v); }
T operator()() const { return v; }
#define REFOP(type, expr...) Modular &operator type (const Modular &rhs) { return expr, *this; }
REFOP(+=, v += rhs.v - MOD, v += MOD & (v >> width)) ; REFOP(-=, v -= rhs.v, v += MOD & (v >> width))
// fits for MOD^2 <= 9e18
REFOP(*=, v = static_cast<T>(1LL * v * rhs.v % MOD)) ; REFOP(/=, *this *= inverse(rhs.v))
#define VALOP(op) friend Modular operator op (Modular a, const Modular &b) { return a op##= b; }
VALOP(+) ; VALOP(-) ; VALOP(*) ; VALOP(/)
Modular operator-() const { return 0 - *this; }
friend bool operator == (const Modular &lhs, const Modular &rhs) { return lhs.v == rhs.v; }
friend bool operator != (const Modular &lhs, const Modular &rhs) { return lhs.v != rhs.v; }
friend std::istream & operator>>(std::istream &I, Modular &m) { T x; I >> x, m = x; return I; }
friend std::ostream & operator<<(std::ostream &O, const Modular &m) { return O << m.v; }
Modular inv() const { return inverse(v); }
Modular qpow(uint64_t p) const {
Modular r = 1, e = *this;
while (p) {
if (p & 1) r *= e;
e *= e;
p >>= 1;
}
return r;
}
private:
constexpr static int width = sizeof(T) * 8 - 1;
T v;
static T inverse(T a) {
// copy from tourist's template
T u = 0, v = 1, m = MOD;
while (a != 0) {
T t = m / a;
m -= t * a; std::swap(a, m);
u -= t * v; std::swap(u, v);
}
assert(m == 1);
return u;
}
};
constexpr int mod = 100000007;
using Mint = Modular<int, mod>;
void solve() {
uint64_t N, K;
int64_t rawA, rawB;
cin >> N >> K >> rawA >> rawB;
debug(N, K, rawA, rawB);
if (rawA == 0) {
// p == 0
cout << 0 << ' ' << 1 << '\n';
return;
}
if (rawA == rawB) {
// p == 1
cout << 1 << ' ' << 1 << '\n';
return;
}
if (rawA * 2 == rawB) {
auto g = gcd(K, N);
debug(K / g, N / g);
cout << Mint(K / g) << ' ' << Mint(N / g) << '\n';
return;
}
// start -- (N - K) -- middle -- (K) -- goal
// R_t is martingale -> E[R_{t'} | information <= t <= t'] = R_t
// E[R_{t+1}] = p * r^{x+1} + (1 - p) * r^{x-1} = r^x
// (r^K - 1) * P(\bar{X}_t = K) + (r^{K-N} - 1) * P(\bar{X}_t = K-N) = 0
// p1 + p2 = 1
// Mint p = A / Mint(B);
// Mint r = (1 - p) / p;
// Mint ca = (r.qpow(K) - 1);
// Mint cb = (1 - r.inv().qpow(N - K));
// debug(ca, cb, r);
// auto ans = ca / (ca + cb);
// cout << ans * 63 << '\n';
// p = A / B
// r = (1 - A/B) / (A/B) = (B - A) / A
// ca / (ca + cb) = (r^K-1) / (r^K-1 + 1-r^{K-N})
// = (r^K-1) / (r^K-r^{K-N})
// = (r^{K+N}-r^N) / (r^{K+N} - r^K)
// = (C/A^{K+N}-C/A^N) / (C/A^{K+N} - C/A^K)
// = (C^{K+N} - C^N A^K) / (C^{K+N} - C^K A^N)
// C^N (C^K - A^K)
// C^K (C^N - A^N)
Mint C = rawB - rawA;
Mint A = rawA;
auto gnk = gcd(K, N);
auto G = C.qpow(gnk) - A.qpow(gnk);
auto p = C.qpow(N - K) * (C.qpow(K) - A.qpow(K)) / G;
auto q = (C.qpow(N) - A.qpow(N)) / G;
// auto p = C.qpow(N - K) * (C.qpow(K) - A.qpow(K));
// auto q = (C.qpow(N) - A.qpow(N));
cout << p << ' ' << q << '\n';
}
signed main() {
cin.tie(nullptr)->sync_with_stdio(false);
int T;
cin >> T;
while (T--)
solve();
}
|