細菌
https://tioj.ck.tp.edu.tw/problems/2116
Description
科學家發現了一種神奇的細菌，他們的繁殖方式很奇怪。這個世界上目前有 $N$ 隻細菌，第 $i$ 隻細菌的體積是 $A_i$ 立方公分。每過一年，對於一隻體積為 $x$ 立方公分的細菌，假設 $x$ 的所有正因數分別是
$$1=d_1<d_2<…<d_k=x$$
那這隻細菌會生出 $k-1$ 隻細菌，體積分別是 $d_1,d_2,…,d_{k-1}$ 立方公分。
科學家發現這種繁殖方式太快了，而且這種細菌是永生的，很快的地球會被這種細菌佔滿!請告訴科學家們，過$K$年後，細菌們的總體積是多少。由於總體積可能太大了，所以請輸出總體積除以 $10^ {9}+7$ 的餘數。
$N \leq 10^6, A_i \leq 10^6, K \leq 10^9$
Solution
因為細菌之間互不干擾，不妨令 $f_t(x)$ 表示一個體積 $x$ 的細菌經過 $t$ 秒的分裂之後分裂出的所有細菌的體積總和是多少。那麼可以列出遞迴式
$$
\begin{cases}
f_t(x) &= \sum _ {d | x} f_{t-1} (d) & \text{ if } t > 0 \\
f_0(x) &= x
\end{cases}
$$
最後的答案是 $\sum _ i f_K(A_i)$。
注意到 $f_0 = \textrm{id}$，而且每次 $t$ 增加 $1$ 就相當於跟 $\mathbb{1}: n\mapsto 1$ 做一次 dirichlet 捲積，即 $f_t = f_{t-1} * \mathbb{1}$。因此 $f_K = \textrm{id} * \mathbb{1}^K$，且因為 $\textrm{id}, \mathbb{1}$ 都是積性函數，$f_K$ 也是一個積性函數。
因此，我們只需要計算出 $f_K$ 的質數冪次的值，就可以 $\mathcal{O}(C)$ 推出所有 $x \leq C$ 的 $f_K(x)$ 了。
現在考慮計算 $f_K(1), f_K(p), f_K(p^2), f_K(p^3), \dots$。因為 $\mathrm{id}(p^k) = p^k$ 以及 $\mathbb{1}(p^k) = 1$，這相當於計算
$$
(1 + px + p^2x^2 + p^3x^3 + \cdots) (1 + x + x^2 + \cdots) ^ K
$$
的前幾項。可以用快速冪計算，也可以直接利用 $( 1+x+x ^ 2+\cdots ) ^ K = (1-x) ^ {-K} = \sum _ i \binom{i+K-1}{i} x ^ i$ 的性質算一次乘法計算。因為 degree 不大（小於 $\log_2{10^6} \approx 20$）所以基本上是寫 $\mathcal{O}(d^2)$ 的多項式乘法。
積性函數和積性函數的乘積是積性函數的證明
設 $f * g = h$，即 $h(n) = \sum _ {d | n} f(d) g(\frac{n}{d})$。
引理：
對於互質的兩數 $a, b$ 以及一個任意的函數 $s$，$\sum _ {d | ab} s(d) = \sum _ {x | a} \sum _ {y | b} s(xy)$
- 證明：令 $D_n \triangleq \{d : d | n\}$。
考慮 $\phi: D_a \times D_b \to D_{ab},\ \phi(x, y) = xy$。$\phi$ 是雙射的，因為存在反函數
$$
\phi^{-1}: D_{ab} \to D_a\times D_b,\ \phi^{-1}(d) = (\gcd(d,a), \gcd(d,b))
$$
因此自然的 $\sum _ {(x,y)\in D_a\times D_b} s(xy) = \sum _ {(x,y)\in D_a\times D_b} s(\phi(x,y)) = \sum _ {d\in D_{ab}} s(d)$ - 反函數的驗證：
若 $x|a, y|b$，則令 $a=ux, b=vy$，那麼 $\gcd(xy, a) = \gcd(xy, ux) = x \cdot \gcd(y, u)$，而 $\gcd(y, u) | \gcd(vy, ux) = \gcd(a, b) = 1$，故 $\gcd(y,u) = 1$，$\gcd(xy, a)$ 確實等於 $x$。$\gcd(xy, b) = y$ 也類似。
由引理便得
$$
\begin{align*}
h(ab)
&= \sum _ {d | ab} f(d) g(\frac{ab}{d}) \\
&= \sum _ {x | a} \sum _ {y | b} f(xy) g(\frac{ab}{xy}) \\
&= \sum _ {x | a} \sum _ {y | b} f(x)f(y) g(\frac{a}{x})g(\frac{b}{y}) \\
&= \left(\sum _ {x | a} f(x)g(\frac{a}{x}) \right)\left(\sum _ {y | b} f(y)g(\frac{b}{y}) \right) \\
&= h(a) h(b)
\end{align*}
$$
對於所有 $\gcd(a, b) = 1$ 都成立。
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
141
| #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 = 1000000007;
using Mint = Modular<int, mod>;
signed main() {
cin.tie(nullptr)->sync_with_stdio(false);
// f(x, k) = \sum _ {d | x} f(d, k - 1)
// f(X, K) = id * 1^K
//
// f_p = 1 / (1 - px) (1 - x)^K
int N, K;
cin >> N >> K;
const int maxn = 1000025;
vector<int> lpf(maxn), pc(maxn), primes;
vector<Mint> f(maxn);
f[1] = 1;
for (int i = 2; i < maxn; i++) {
if (lpf[i] == 0) primes.emplace_back(lpf[i] = pc[i] = i);
for (int p : primes) {
if (i * p >= maxn) break;
lpf[i * p] = p;
if (lpf[i] == p) {
pc[i * p] = pc[i] * p;
break;
} else {
pc[i * p] = p;
}
}
}
vector<Mint> choose(60, 1);
for (int i = 0; i < 60; i++) {
// binom(i + K - 1, i);
for (int j = 1; j <= i; j++)
choose[i] *= K + j - 1;
for (int j = 1; j <= i; j++)
choose[i] /= j;
}
for (int p : primes) {
vector<Mint> u;
for (int64_t prod = 1, c = 0; prod < maxn; prod *= p, ++c) {
u.emplace_back(prod);
}
vector<Mint> v(u.size());
for (size_t i = 0; i < v.size(); i++) {
for (size_t j = 0; j <= i; j++) {
v[i] += u[j] * choose[i - j];
}
}
for (int64_t prod = 1, c = 0; prod < maxn; prod *= p, ++c) {
f[prod] = v[c];
}
}
for (int i = 2; i < maxn; i++)
if (i != pc[i])
f[i] = f[pc[i]] * f[i / pc[i]];
Mint ans = 0;
for (int i = 0; i < N; i++) {
int x;
cin >> x;
ans += f[x];
}
cout << ans << '\n';
}
|
很懶惰就貼了 modint 的模板。pc 是維護 least prime factor 在 $n$ 裡面的冪次究竟是多少。記得初始化 f[1] = 1。