[bzoj3309] DZY Loves Math
Description
对于正整数n，定义f(n)为n所含质因子的最大幂指数。例如f(1960)=f(2^3 * 5^1 * 7^2)=3, f(10007)=1, f(1)=0。
给定正整数a,b，求sigma(sigma(f(gcd(i,j)))) (i=1..a, j=1..b)。
Input
第一行一个数T，表示询问数。
接下来T行，每行两个数a,b，表示一个询问。
Output
对于每一个询问，输出一行一个非负整数作为回答。
Sample Input
4
7558588 9653114
6514903 4451211
7425644 1189442
6335198 4957
Sample Output
35793453939901
14225956593420
4332838845846
15400094813
HINT
【数据规模】
T<=10000
1<=a,b<=10^7
solution
前置知识：莫比乌斯反演
题目让求的是：
[ans=sum_{i=1}^nsum _{j=1}^mf(gcd(i,j))
]
然后通过莫比乌斯反演可以得到：
[ans=sum_{T=1}^{min(n,m)}lfloorfrac{n}{T}
floorlfloorfrac{m}{T}
floorsum_{d|T}f(d)mu(frac{T}{d})
]
然后令式子后面一块为(g)，即：
[g(T)=sum_{d|T}f(d)mu(frac{T}{d})
]
显然这种题不线筛出这个东西不让过，
然后分解下(d)和(T)：
[T=prod_{i=1}^ka_i^{p_i},d=prod_{i=1}^ka_i^{x_i}
]
考虑(mu(frac{T}{d}))，显然对于任意(i),满足(x_i=p_i)或者(x_i=p_i-1)，否则(mu)为0。
设(res=f(T))，有(num)个(p_i)等于(res)，然后分类讨论：
若(num=k)，则(f(d)=res)的答案为：
[egin{align}
&sum_{f(d)=res}f(d)mu(frac{T}{d})\
=&res*sum_{f(d)=res}mu(frac{T}{d})\
=&res*sum_{i=1}^k(-1)^{k-i}inom{k}{i}\
=&res*(1+(-1))^k-res*(-1)^k\
=&res*(-1)^{k+1}
end{align}
]
(f(d)=res-1)的答案为：
[egin{align}
&sum_{f(d)=res-1}f(d)mu(frac{T}{d})=(res-1)*(-1)^k
end{align}
]
综合起来，答案为(res*(-1)^{k+1}+(res-1)*(-1)^k=(-1)^{k+1})。
否则，即(num eq f(T))，则(f(d)=res)的答案为：
[egin{align}
&sum_{f(d)=res}f(d)mu(frac{T}{d})\
=&res*sum_{i=1}^{num}(-1)^{k-i}inom{num}{i}*sum_{i=0}^{k-num}(-1)^iinom{k-num}{i}\
=&res*sum_{i=1}^{num}(-1)^{k-i}inom{num}{i}*0=0
end{align}
]
同理(f(d)=res-1)的答案也为0，所以综合起来答案为0。
由于线筛的时候每个数是被最小的质因子筛掉的，所以考虑线筛的时候对于每个数最小质因子出现了多少次，以及不算这个质因子剩下的是多少，然后对于(g)讨论下就行了。
剩下的随便数论分块下就做完了，时间复杂度(O(n+qsqrt{n}))。
#include<bits/stdc++.h>
using namespace std;
#define ll long long
void read(int &x) {
x=0;int f=1;char ch=getchar();
for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
void print(int x) {
if(x<0) x=-x,putchar('-');
if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}
const int maxn = 1e7+1;
int f[maxn],pri[maxn/20],tot,num[maxn],lst[maxn];
bool vis[maxn];
void sieve() {
for(int i=2;i<maxn;i++) {
if(!vis[i]) pri[++tot]=i,f[i]=num[i]=lst[i]=1;
for(int j=1;j<=tot&&i*pri[j]<maxn;j++) {
vis[i*pri[j]]=1;
if(!(i%pri[j])) {
lst[i*pri[j]]=lst[i];
num[i*pri[j]]=num[i]+1;
if(lst[i*pri[j]]==1) f[i*pri[j]]=1;
else f[i*pri[j]]=(num[lst[i*pri[j]]]==num[i*pri[j]]?-f[lst[i]]:0);
break;
}
lst[i*pri[j]]=i;
num[i*pri[j]]=1;
f[i*pri[j]]=(num[i]==1?-f[i]:0);
}
}
//for(int i=1;i<=30;i++) printf("%d %d %d %d
",i,f[i],lst[i],num[i]);;
for(int i=1;i<maxn;i++) f[i]=f[i-1]+f[i];
//for(int i=1;i<=20;i++) write(f[i]);
}
signed main() {
int t;read(t);
//int PRE=clock();
sieve();
//cerr << (double)(clock()-PRE)/CLOCKS_PER_SEC << endl;
while(t--) {
int n,m;read(n),read(m);
int T=1;ll ans=0;
while(T<=n&&T<=m) {
int pre=T;T=min(n/(n/T),m/(m/T));
ans+=1ll*(n/T)*(m/T)*(f[T]-f[pre-1]);T++;
}
printf("%lld
",ans);
}
return 0;
}