对于正整数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(nsqrt{n})$的，后来才发现，这其实是$O(n*ln(n))$的，所以暴力计算f*mu应该也能过
2018.01.25 upd:一开始不觉得数列求和分析复杂度可以积分（好像在项数少的时候不太拟合），现在看来好像可以（项数多的时候比较好）
2018.07.12 upd:这里顺便记录一下枚举除法套枚举除法是$O(N^{3over 4})$的，积分证得
代码：
1 #include<cmath> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #define LL long long 6 using namespace std; 7 const int MAXN=1e7; 8 LL f_mu[MAXN+5]; 9 int pri[MAXN],f[MAXN+5],mu[MAXN+5],p[MAXN+5],cnt; 10 bool vis[MAXN]; 11 LL sum[MAXN+5]; 12 void prime(); 13 LL Sqr(LL ,int ); 14 LL work(int ,int ); 15 int main() 16 { 17 int i,j,k,n,a,b; 18 prime(); 19 sum[0]=f_mu[1]=sum[1]=0; 20 for(i=2;i<=MAXN;i++){ 21 if(mu[i]==1) 22 f_mu[i]=-1; 23 if(mu[i]==-1) 24 f_mu[i]=1; 25 if(mu[i]==0){ 26 if(i==Sqr(p[i],f[i])) 27 f_mu[i]=f_mu[p[i]]; 28 else 29 f_mu[i]=0; 30 } 31 sum[i]=sum[i-1]+f_mu[i]; 32 } 33 scanf("%d",&n); 34 for(i=1;i<=n;i++){ 35 scanf("%d%d",&a,&b); 36 printf("%lld ",work(a,b)); 37 } 38 return 0; 39 } 40 void prime(){ 41 int i,j; 42 vis[1]=true,mu[1]=1; 43 for(i=2;i<=MAXN;i++){ 44 if(!vis[i]){ 45 f[i]=1,p[i]=i,mu[i]=-1; 46 pri[++cnt]=i; 47 } 48 for(j=1;j<=cnt&&pri[j]*i<=MAXN;j++){ 49 vis[i*pri[j]]=true; 50 if(i%pri[j]){ 51 mu[i*pri[j]]=-mu[i]; 52 f[i*pri[j]]=f[i],p[i*pri[j]]=p[i]*pri[j]; 53 } 54 else{ 55 mu[i*pri[j]]=0; 56 f[i*pri[j]]=f[i]+(i%Sqr(pri[j],f[i])==0); 57 p[i*pri[j]]=p[i]; 58 break; 59 } 60 } 61 } 62 } 63 LL Sqr(LL x,int n){ 64 LL ret=1; 65 while(n){ 66 if(n&1) 67 ret*=x; 68 n>>=1,x*=x; 69 } 70 return ret; 71 } 72 LL work(int a,int b){ 73 int MIN=min(a,b),i,las; 74 LL ret=0; 75 for(las=0,i=1;i<=MIN;las=i,i=min(a/(a/(las+1)),b/(b/(las+1)))){ 76 ret+=(sum[i]-sum[las])*(1ll*a/i)*(1ll*b/i); 77 if(i==MIN)break; 78 } 79 return ret; 80 }