[bzoj3309] DZY Loves Math

题目大意

对于正整数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}^a \sum_{j=1}^b f((i,j))$

T≤10000
1≤a,b≤$10^7$

分析

看到两个sigma和(i,j)，就要条件反射想到莫比乌斯反演。
首先令a≤b
枚举公约数，得到：

$$Ans=\sum_{d=1}^a f(d)\sum_{i=1}^{\lfloor\frac{a}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{b}{d}\rfloor}[(i,j)=1]$$
$$Ans=\sum_{d=1}^a f(d)\sum_{d'=1}^{\lfloor\frac{a}{d}\rfloor} \mu(d') \lfloor \frac{a}{dd'} \rfloor \lfloor \frac{b}{dd'} \rfloor$$
令T=dd’，得到：
$$Ans=\sum_{T=1}^a \lfloor \frac{a}{T} \rfloor \lfloor \frac{b}{T} \rfloor\sum_{d|T} f(d)*mu(\frac{T}{d})$$
设 $g[T]=\sum_{d|T} f(d)*mu(\frac{T}{d})$，那么只要预处理g数组的前缀和，就可以以根号的复杂度询问了。
g数组的预处理看起来是带log的。但是根据莫比乌斯函数的性质，如果 $\frac{T}{d}$存在平方因子，函数值是等于0的，也就是对答案没有贡献。
那么设 $T=p1^{k1}*p2^{k2}*...*pm^{km}$，T质因数的最大幂是k，那么只有ki=k的质因数有用。又可以设 $\frac{T}{d}=p1^{a1}*p2^{a2}*...*pm^{am}$，其中 $ai∈[0,1]$。
可以发现f(d)只能取到k,k-1，现在令其中一个满足ki=k的质因数为d的最大幂，如果f(d)=k，那么ai=0，其它为0或1均可。然而一个a取0，就相当于给莫比乌斯函数乘1，取1就是乘-1。所以最终答案乘的系数是0。特殊情况：如果T是质数，乘的系数是1（因为没有其它质因数了）。
如果f(d)=k-1，那么ai=1，所有其它满足ki=k的质因数也要让对应的a值取1，这时剩下的质因数也和上面一样，最后得到的系数是0。特殊情况：如果每个ki都等于k，那么由于没有剩下可以取0、1的质因数，它的系数也是1。
这样就可以线性预处理了。

#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;

const int N=1e7;

typedef long long LL;

int T,mu[N+5],tot,p[N],f[N+5];

LL s[N+5],ans;

bool bz[N+5];

char c;

int read()
{
    for (c=getchar();c<'0' || c>'9';c=getchar());
    int x=c-48;
    for (c=getchar();c>='0' && c<='9';c=getchar()) x=x*10+c-48;
    return x;
}

int main()
{
    mu[1]=1;
    for (int i=2;i<=N;i++)
    {
        if (!bz[i])
        {
            mu[i]=-1;
            p[tot++]=i;
        }
        for (int j=0;j<tot && i*p[j]<=N;j++)
        {
            int I=i*p[j];
            bz[I]=1;
            if (i%p[j]==0)
            {
                mu[I]=0; break;
            }
            mu[I]=-mu[i];
        }
    }
    for (int i=2;i<=N;i++) if (mu[i]!=0)
    {
        for (LL j=i;j<=N;j*=i) s[j]=-mu[i];
    }
    for (int i=1;i<=N;i++) s[i]+=s[i-1];
    T=read();
    while (T--)
    {
        int a=read(),b=read();
        if (a>b) a^=b^=a^=b;
        ans=0;
        for (int i=1,j;i<=a;i=j+1)
        {
            j=min(a/(a/i),b/(b/i));
            ans+=(s[j]-s[i-1])*(a/i)*(b/i);
        }
        printf("%lld\n",ans);
    }
    return 0;
}