轨道 反演DP

发表于 更新于 分类于 阅读次数：

貌似是jzoj的题… 题截了个图因为会在主页显示很违和就放链接里了点这里可以看 很有意思的一个题，之前看过一次DP状态设计，自己yy了一个反演… 首先是个DP。$(v,k)=1$，说明$(\prod a_i,k)=k$。令$f_{i,j}$为前$i$个，和$k$的最大公约数为$j$。发现这个状态太大了，把$j$改成第$j$个因子。 $$f_{i,j}=\sum f_{i-1,k_1}f_{1,k_2}$$ 其中，$k_1$是枚举的，$k_2$是算出来的，$k_1k_2$两个因子的积为因子$j$。 问题就在于如何求出$f_{1,x}$ $$f_{1,x} = \sum_{i=1}^m [gcd(i,k) = x]$$ 题解里说是容斥，蒟蒻不会容斥，，，搞了个反演 令$f(x)$为$\sum_{i=1}^m [gcd(i,k) = x]$，$F(x)$为$\sum_{i=1}^m [x \mid gcd(i,k)]$ $$F(x) = \sum_{i=1}^m [x \mid gcd(i,k)] = \lfloor \frac{m}{x} \rfloor [x \mid k]$$ $$f(x) = \sum_{x \mid d} \mu(\frac{d}{x}) F(d) = \sum_{T=1}^{\lfloor \frac{m}{x} \rfloor} \mu(T) F(xT)$$ 化简一下，得到这个式子： $$f(x) = \sum_{T=1}^{\lfloor \frac{m}{x} \rfloor} \mu(T) [xT \mid k] \lfloor \frac{m}{xT} \rfloor$$ 然而，还不能整除分块。两个问题，首先本题$m$的范围是$10^9$，其次，$[xT \mid k]$如何处理。 不过实际上，$k$的范围是$10^7$。也就是说，我们只需要处理到$10^7$就可以了。这可以用一个naive的线性筛随便搞定。至于$[xT \mid k]$，把它乘进$\mu$里即可。 因为用vector被卡常了（悲

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
#pragma GCC optimize(3)
#pragma GCC target("sse","sse2","sse3","sse4","avx","avx2","popcnt") 

#include <bits/stdc++.h>

const int N = 7000, M = 1e7 + 233, P = 10007;

int n, m, k;
int part[N], tot;
int pri[M], vis[M], tt;
int mu[M];
int qaq[M];

void init1() {
    part[++tot] = mu[1] = 1; qaq[1] = 1;
    for (int i = 2; i <= k; ++i) {
        if (k % i == 0)
            part[++tot] = i, qaq[i] = tot;
        if (!vis[i]) {
            mu[i] = -1;
            pri[++tt] = i;
        }
        for (int j = 1; j <= tot && i * pri[j] <= k; ++j) {
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0) {
                mu[i * pri[j]] = 0;
                break;
            }
            mu[i * pri[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= k; ++i)
        mu[i] = mu[i - 1] + (k % i == 0) * mu[i];
}

int f[N / 2][N];
std::vector<std::pair<int, int> > mul[N];

long long calc(int x) {
    long long ret = 0;
    for (int l = 1, r; l <= x; l = r + 1) {
        r = x / (x / l);
        ret += (mu[std::min(k, r)] - mu[std::min(k, l - 1)]) * (x / l);
    }
    return ret % P;
}

void init2() {
    for (int i = 1; i <= tot; ++i) {
        if (m / part[i] == 0)
            break;
        f[1][i] = calc(m / part[i]);
    }
    for (int i = 1; i <= tot; ++i)
        for (int j = 1; j <= tot; ++j) {
            if ((long long)part[i] * part[j] > k)
                break;
            int qwq = part[i] * part[j];
            if (qaq[qwq])
                mul[qaq[qwq]].push_back({ i, j });
        }
}

void solve() {
    for (int i = 2; i <= n; ++i)
        for (int j = 1; j <= tot; ++j) {
            long long sum = 0;
            for (auto y : mul[j]) {
                sum += f[i - 1][y.first] * f[1][y.second];
            }
            f[i][j] = (f[i][j] + sum) % P;
        }

    std::cout << f[n][tot] << std::endl;
}

int main() {
    std::cin >> n >> m >> k;
    init1();
    init2();
    solve();
    return 0;
}