hdu 4746 Mophues 莫比乌斯 分块优化

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

210 10 010 10 1 

 

Sample Output

6393 

然后由莫比乌斯反演，我们知道f(d)=g(d)*u(1)+g(2*d)*u(2)+g(3*d)*u(3)+....

考虑结果ans，ans为所有f(d)且F(d)<=P的和，

我们枚举1<=i<=n，如果i是某个d的倍数且F(d)<=P，那么ans+=g(i)*u(i/d)=[n/i]*[m/i]*u(i/d)。那么这个怎么计算能更快一点？

我们设G(i)为容斥因子：G(i)=sum{u(i/d) | F(d)<=P} 这个值可以nlogn预处理出来，然后我们只需要ans+=G(i)*[n/i]*[m/i]即可

这样的话总的复杂度为O(n*q)还是会T的样子

然后我们注意到[n/i]*[m/i]在一定的范围内是不变的，这个范围是[i,min(n/(n/i),m/(m/i)]，这样我们可以预处理出G(i)的前缀和，然后加快运算（复杂度网上说是sqrt(n)的。。。）

这样总的复杂度是O(q*sqrt(n)+nlog(n))大概这样.

#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<vector>
#include<set>
#include<map>
#include<queue>
#include<bitset>
#define fi first
#define se second
#define pii pair<int,int>
#define ll long long
#define inf 1<<30
#define eps 1e-8
using namespace std;
const int maxn=500010;
int num[maxn];
int cnt[19][maxn];
int n,m,p;
int mu[maxn],prime[maxn],tot;
bool check[maxn];
void getmu()
{
int N=500000;
memset(check,0,sizeof(check));
mu[1]=1,tot=0;
for(int i=2;i<=N;i++){
if(!check[i])
prime[tot++]=i,mu[i]=-1;
for(int j=0;j<tot;j++){
if(i*prime[j]>N)
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=1;i<=N;i++) {
int u=i;
int v=0;
for(int j=0;j<tot;j++) {
if(prime[j]*prime[j]>u)
break;
while(u%prime[j]==0) {
u/=prime[j];
v++;
}
}
if(u>1) v++;
num[i]=v;
}
}
void init()
{
memset(cnt,0,sizeof(num));
int N=500000;
for(int i=1;i<=N;i++) {
for(int j=i;j<=N;j+=i) {
cnt[num[i]][j]+=mu[j/i];
}
}
for(int i=1;i<=N;i++) {
for(int j=1;j<=18;j++) {
cnt[j][i]+=cnt[j-1][i];
}
}
for(int i=0;i<=18;i++) {
for(int j=1;j<=N;j++) {
cnt[i][j]+=cnt[i][j-1];
}
}
}
int main()
{
getmu();
init();
int t;
scanf("%d",&t);
while(t--) {
scanf("%d%d%d",&n,&m,&p);
if(p>18) {
printf("%I64d\n",(ll)n*m);
continue;
}
ll ans=0;
if(n>m) swap(n,m);
for(int i=1;i<=n;) {
int j=min(n/(n/i),m/(m/i));
//cout<<i<<" "<<j<<endl;
ans+=(ll)(cnt[p][j]-cnt[p][i-1])*(n/i)*(m/i);
i=j+1;
}
printf("%I64d\n",ans);
}
return 0;
}