Description
As we know, any positive integer C ( C >= 2 ) can be written as the multiply of some prime numbers:
C = p1×p2× p3× ... × pk
which p1, p2 ... pk are all prime numbers.For example, if C = 24, then:
24 = 2 × 2 × 2 × 3
here, p1 = p2 = p3 = 2, p4 = 3, k = 4
Given two integers P and C. if k<=P( k is the number of C's prime factors),
we call C a lucky number of P.
Now, XXX needs to count the number of pairs (a, b), which 1<=a<=n ,
1<=b<=m, and gcd(a,b) is a lucky number of a given P ( "gcd"
means "greatest common divisor").
Please note that we define 1 as lucky number of any non-negative integers
because 1 has no prime factor.
Input
The first line of input is an integer Q
meaning that there are Q test cases.
Then Q lines follow, each line is a test case and each test case contains three
non-negative numbers: n, m and P (n, m, P <= 5×10 5. Q
<=5000).
Output
For each test case, print the number of pairs (a, b), which 1<=a<=n , 1<=b<=m, and gcd(a,b) is a lucky number of P.
Sample Input
2
10 10 0
10 10 1
Sample Output
63
93
题目大意就是求k <= p，满足gcd(x, y) == k的对数。
首先可以预处理出任意数i其质因子个数luck[i]。
这个的话可以初始化luck[1] = 0，然后类似于bfs的过程，把当前i能到达的状态prime*i标记，即luck[prime*i] = luck[i]+1。然后用vis数组来避免重复标记。
然后就是根据莫比乌斯，之前写过类似的题解，结论就是gcd(x, y) == k的对数f(k)满足：
f(k) = sum(u(j/k)*(m/j)*(n/j)) (k | j)
这样对于满足的k，就是上式的k求和，之前也说过sum(sum(u(j/k)*(m/j)*(n/j)) (k | j))中两个sum的次序可以进行交换的。因为一个k对应好多j，等同于一个j对应好多k。
这样就是求满足条件的k的sum((m/j)*(n/j)*sum(u(j/k))) (k | j)
这样可以预先处理出所有j对应的sum(u(j/k))，
但是还不够，这样对于重复的p，会重复计算，这样需要多开一维sum(p, u(j/k))表示满足p条件下的sum，相当于把所有情况的sum离线下来。这样做主要的前提是p很小，因为5*10^5内质因子的个数最多有log2(5*10^5)个，大约18个。
在求sum[p][j]时，可以先计算出gcd(x, y) = p的情况，然后求前缀和就是gcd(x, y) <= p的情况了。
不过这题到这里为止还没有结束。。。
预处理复杂度能过去，不过后面的数据规模比较大，O(n)时间内过不去。
需要对上述答案进行分组加速。
因为对于n/j来说，必然会存在一种情况，n/j == n/(j+1) = …..
就是连续的几个j的结果一致，因为这里的除法取的是下整。
这样的话，若n/j == p，
那么p*j <= n < (p+1)*j
如果x也满足n/j == n/x == p
那么自然p*x <= n < (p+1)*x
于是x <= n/p。
也就是x最大可以取n/p，由于取的是下整，可以放心取等号。
那么从j到n/(n/j)区间内的n/j的结果都是一致的。然后j到m/(m/j)的m/j的结果又是一致的。
自然从j到min(n/(n/j), m/(m/j))内(n/j)*(m/j)的结果一致。
那么这一段应该直接计算sum[p](j~ min(n/(n/j), m/(m/j)))* (n/j)*(m/j),然后可以事先对sum[p][j]求和，就可以通过作差求出前半段。
最后需要注意的一点是当p >= log2(len)时，所有(x, y)对就都能满足条件了，答案自然是m*n。
代码：
#include <iostream> #include <cstdio> #include <cstdlib> #include <cstring> #include <cmath> #include <queue> #define LL long long using namespace std; const int maxN = 500010; int n, m, p; int luck[maxN], sum[20][maxN]; int prime[maxN], cnt, u[maxN]; bool vis[maxN]; void mobius() { memset(vis, false,sizeof(vis)); u[1] = 1; cnt = 0; for(int i = 2; i < maxN; i++) { if(!vis[i]) { prime[cnt++] = i; u[i] = -1; } for(int j = 0; j < cnt && i*prime[j] < maxN; j++) { vis[i*prime[j]] = true; if(i%prime[j]) u[i*prime[j]] = -u[i]; else { u[i*prime[j]] = 0; break; } } } } void calLuck() { memset(vis, false, sizeof(vis)); luck[1] = 0; queue<int> q; q.push(1); vis[1] = true; int k, v; while (!q.empty()) { k = q.front(); q.pop(); for (int i = 0; i < cnt; ++i) { if ((LL)prime[i]*k >= maxN) break; v = prime[i]*k; if (!vis[v]) { luck[v] = luck[k]+1; q.push(v); vis[v] = true; } } } } void init() { mobius(); calLuck(); memset(sum, false, sizeof(sum)); //计算sum[p][j]的值，gcd个数为p的i的倍数j系数的总和 for (int i = 1; i < maxN; ++i) { for (LL j = i; j < maxN; j += i) sum[luck[i]][j] += u[j/i]; } //计算sum[p][j]的值，gcd个数小于等于p的i的倍数j系数的总和 for (int i = 1; i < maxN; ++i) { for (int k = 1; k < 20; ++k) sum[k][i] += sum[k-1][i]; } //计算sum[p][j]的j前缀和，用于分组加速 for (int k = 0; k < 20; ++k) for (int i = 1; i < maxN; ++i) sum[k][i] += sum[k][i-1]; } void work() { int len = min(m, n); if (p >= log2(len)) { printf("%I64d ", (LL)m*n); return; } LL ans = 0; for (int i = 1, j; i <= len; i = j+1) { j = min(n/(n/i), m/(m/i)); ans += (LL)(sum[p][j]-sum[p][i-1])*(m/i)*(n/i); } printf("%I64d ", ans); } int main() { //freopen("test.in", "r", stdin); init(); int T; scanf("%d", &T); for (int times = 1; times <= T; ++times) { scanf("%d%d%d", &n, &m, &p); work(); } return 0; }