对于正整数n，定义f(n)为n所含质因子的最大幂指数。例如f(1960)=f(2^3 * 5^1 * 7^2)=3, f(10007)=1, f(1)=0。
给定正整数a,b，求:
$$\sum_{i=1}^{i<=a}\sum_{j=1}^{j<=b}f(gcd(i,j))$$
bzojP3309
http://www.lydsy.com/JudgeOnline/problem.php?id=3309
化式子：
$$\sum_{i=1}^{i<=a}\sum_{j=1}^{j<=b}f(gcd(i,j))$$
$$\sum_{x=1}^{min(a,b)}f(x)·\sum_{x|d}\mu({d \over x})({a \over d})({b \over d})$$
$$\sum_{x=1}^{min(a,b)}f(x)·\sum_{d=1}^{min(a,b) \over x}\mu(d)({a \over dx})({b \over dx})$$
$$\sum_{dx=1}^{min(a,b)}({a \over dx})({b \over dx})\sum_{d=1}^{d<=dx}f(x)*\mu(d)$$
$$\sum_{x=1}^{min(a,b)}({a \over x})({b \over x})(f*\mu)(x)$$
这个化式子的过程旨在尽可能地把无关a,b的，可预处理的部分提出来；
现在只要预处理出f与$\mu$的卷积即可通过枚举除法做到单次$O(\sqrt{n})的效率$
f和$\mu$可以通过线性筛求出，
如何在较好的时间内求出他们的卷积呢；
不会，
不过打表可以发现：
$$当\mu(i)=1时，(f*\mu)(i)=-1$$
$$当\mu(i)=-1时，(f*\mu)(i)=1$$
$$当\mu(i)=0时，若i的所有质因子齐次，次数为k，则(f*\mu)(i)=(f*\mu)(^{k}\sqrt{i})，否则(f*\mu)(i)=0$$
这样处理好$f*\mu$，然后枚举除法即可；
由于笔者太弱，于是直接上了单次处理$O(log_2log_2n)$的线性筛
2018.01.25 upd:我原来以为暴力计算卷积函数是$O(n\sqrt{n})$的，后来才发现，这其实是$O(n*ln(n))$的，所以暴力计算f*mu应该也能过
2018.01.25 upd:一开始不觉得数列求和分析复杂度可以积分（好像在项数少的时候不太拟合），现在看来好像可以（项数多的时候比较好）
2018.07.12 upd:这里顺便记录一下枚举除法套枚举除法是$O(N^{3\over 4})$的，积分证得
代码：
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long
using namespace std;
const int MAXN=1e7;
LL f_mu[MAXN+];
int pri[MAXN],f[MAXN+],mu[MAXN+],p[MAXN+],cnt;
bool vis[MAXN];
LL sum[MAXN+];
void prime();
LL Sqr(LL ,int );
LL work(int ,int );
int main()
{
int i,j,k,n,a,b;
prime();
sum[]=f_mu[]=sum[]=;
for(i=;i<=MAXN;i++){
if(mu[i]==)
f_mu[i]=-;
if(mu[i]==-)
f_mu[i]=;
if(mu[i]==){
if(i==Sqr(p[i],f[i]))
f_mu[i]=f_mu[p[i]];
else
f_mu[i]=;
}
sum[i]=sum[i-]+f_mu[i];
}
scanf("%d",&n);
for(i=;i<=n;i++){
scanf("%d%d",&a,&b);
printf("%lld\n",work(a,b));
}
return ;
}
void prime(){
int i,j;
vis[]=true,mu[]=;
for(i=;i<=MAXN;i++){
if(!vis[i]){
f[i]=,p[i]=i,mu[i]=-;
pri[++cnt]=i;
}
for(j=;j<=cnt&&pri[j]*i<=MAXN;j++){
vis[i*pri[j]]=true;
if(i%pri[j]){
mu[i*pri[j]]=-mu[i];
f[i*pri[j]]=f[i],p[i*pri[j]]=p[i]*pri[j];
}
else{
mu[i*pri[j]]=;
f[i*pri[j]]=f[i]+(i%Sqr(pri[j],f[i])==);
p[i*pri[j]]=p[i];
break;
}
}
}
}
LL Sqr(LL x,int n){
LL ret=;
while(n){
if(n&)
ret*=x;
n>>=,x*=x;
}
return ret;
}
LL work(int a,int b){
int MIN=min(a,b),i,las;
LL ret=;
for(las=,i=;i<=MIN;las=i,i=min(a/(a/(las+)),b/(b/(las+)))){
ret+=(sum[i]-sum[las])*(1ll*a/i)*(1ll*b/i);
if(i==MIN)break;
}
return ret;
}