思路:不难看出可以看出p最多不超过19.
首先考虑这个问题的简化版,求gcd(a,b)小于等于d的数对的数量。
对于这个问题,由莫比乌斯反演我们可以求出gcd为i的数对数量然后从1到j累加起来,但是d太大的话时间复杂度无法承受。
但是我们注意到,对于每个gcd = i,设A(d):gcd(a, b)=d的有多少种,设B(j): gcd(a, b)是j的倍数的有多少种,
有如下形式:
A(1) = μ(1)*B(1) + μ(2)*B(2) + μ(3)*B(3) + ... + μ(p1*p2...)*B(p1*p2...)
A(2) = μ(1)*B(1*2) + μ(2)*B(2*2) + μ(3)*B(3*2) + ... + μ(p1*p2..)*B(p1*p2..*2)
A(d) = μ(1)*B(1*d) + μ(2)*B(2*d) + μ(3)*B(3*d) + ... + μ(p1*p2..)*B(p1*p2..*d)
于是答案就是ans = A(1)+A(2)+...+A(d) = F(1)*B(1) + F(2)*B(2) + ... + F(p1*p2..)*B(p1*p2..)。
而每项的系数和要求的d无关,我们可以用类似筛法nlogn的时间预处理出来,那么这个问题就可以解决了。
考虑本题,只是多了一个p这个参数,那么我们只需要给系数再加一维,表示素因子个数小于等于p时的系数值,这个可以通过一个累加求出,这样我们就预处理出了所要用到的系数,但是这样还是会t,在计算答案时还要加上一个分段优化,即把(n/i)*(m/i)值相同的合并起来并对系数预处理出一个前缀和,也就是将[i, n/(n/i), m/(m/i)]合并起来。
#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#include<vector>
#include<map>
#include<queue>
#include<stack>
#include<string>
#include<map>
#include<set>
#include<ctime>
#define eps 1e-6
#define LL long long
#define pii pair<int, int>
//#pragma comment(linker, "/STACK:1024000000,1024000000")
using namespace std;
const int MAXN = 500000+50;
bool check[MAXN+10];
int prime[MAXN+10];
int mu[MAXN+10], pri_fac[MAXN+10], F[MAXN][20];
int Sum[MAXN][30];
void Moblus()
{
mu[1] = 1;
int tot = 0;
for(int i = 2; i <= MAXN; i++)
{
if( !check[i] )
{
prime[tot++] = i;
mu[i] = -1;
}
for(int j = 0; j < tot; j++)
{
if(i * prime[j] > MAXN) break;
check[i * prime[j]] = true;
if( i % prime[j] == 0)
{
mu[i * prime[j]] = 0;
break;
}
else
{
mu[i * prime[j]] = -mu[i];
}
}
}
for(int i = 0; i < tot; i++) {
for(int j = prime[i]; j < MAXN; j += prime[i]) {
int t = j;
while(t && t%prime[i]==0) pri_fac[j]++, t /= prime[i];
}
}
}
void init() {
for(int i = 1; i < MAXN; i++) {
for(int j = i; j < MAXN; j += i) {
F[j][pri_fac[i]] += mu[j/i];
}
}
for(int i = 1; i < MAXN; i++) {
for(int j = 1; j < 20; j++) F[i][j] += F[i][j-1];
}
for(int i = 1; i < MAXN; i++) {
for(int j = 0; j < 20; j++) {
Sum[i][j] = Sum[i-1][j] + F[i][j];
}
}
}
int main() {
//freopen("input.txt", "r", stdin);
Moblus();
init();
int q, n, m, p;
cin >> q;
while(q--) {
scanf("%d%d%d", &n, &m, &p);
LL ans = 0;
int t = min(n, m), last;
p = min(19, p);
for(int i = 1; i <= t; i = last+1) {
last = min(n/(n/i), m/(m/i));
ans += (LL)(Sum[last][p]-Sum[i-1][p])*(n/i)*(m/i);
}
cout << ans << endl;
}
return 0;
}