[BZOJ3309]DZY Loves Math

题目大意：
　　定义$f(n)$为$n$所含质因子的最大幂指数，共$q(qleq10000)$组询问，对于给定的$n,m(n,mleq10^7)$，求$sum_{i=1}^nsum_{j=1}^mf(gcd(i,j))$。
思路：
$$
egin{align*}
原式&=sum_{d=1}^{min(n,m)}f(d)sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{m}{d} floor}[gcd(i,j)=1]\
&=sum_{d=1}^{min(n,m)}f(d)sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{m}{d} floor}sum_{p|gcd(i,j)}mu(p)\
&=sum_{d=1}^{min(n,m)}f(d)sum_{p=1}^{min(lfloorfrac{n}{d} floor,lfloorfrac{m}{d} floor)}mu(p)sum_{i=1}^{lfloorfrac{n}{dp} floor}sum_{j=1}^{lfloorfrac{m}{dp} floor}1\
&=sum_{d=1}^{min(n,m)}f(d)sum_{p=1}^{min(lfloorfrac{n}{d} floor,lfloorfrac{m}{d} floor)}mu(p)lfloorfrac{n}{dp} floorlfloorfrac{m}{dp} floor\
&=sum_{d=1}^{min(n,m)}lfloorfrac{n}{d} floorlfloorfrac{m}{d} floorsum_{p|d}f(p)mu(frac{d}{p})
end{align*}
$$
　　令$g(d)=sum_{p|d}f(p)mu(frac{d}{p})$，则
$$
原式=sum_{d=1}^{min(n,m)}lfloorfrac{n}{d} floorlfloorfrac{m}{d} floor g(d)
$$
　　考虑预处理$g(d)$。（具体过程略，参见https://www.cnblogs.com/RogerDTZ/p/8227485.html）

#include<cstdio>
#include<cctype>
#include<algorithm>
typedef long long int64;
inline int getint() {
    register char ch;
    while(!isdigit(ch=getchar()));
    register int x=ch^'0';
    while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0');
    return x;
}
const int N=10000001,M=N;
bool vis[N];
int p[M],a[N],b[N];
int64 g[N],sum[N];
inline void sieve() {
    for(register int i=2;i<N;i++) {
        if(!vis[i]) {
            p[++p[0]]=i;
            a[i]=g[i]=1;
            b[i]=i;
        }
        for(register int j=1;j<=p[0]&&i*p[j]<N;j++) {
            vis[i*p[j]]=true;
            if(i%p[j]==0) {
                a[i*p[j]]=a[i]+1;
                b[i*p[j]]=b[i]*p[j];
                if(i==b[i]) {
                    g[i*p[j]]=1;
                } else {
                    g[i*p[j]]=(a[i*p[j]]==a[i/b[i]])?-g[i/b[i]]:0;
                }
                break;
            } else {
                a[i*p[j]]=1;
                b[i*p[j]]=p[j];
                g[i*p[j]]=(a[i]==1)?-g[i]:0;
            }
        }
    }
    for(register int i=1;i<N;i++) {
        sum[i]=sum[i-1]+g[i];
    }
}
int main() {
    sieve();
    for(register int T=getint();T;T--) {
        const int n=getint(),m=getint(),lim=std::min(n,m);
        int64 ans=0;
        for(register int i=1,j;i<=lim;i=j+1) {
            j=std::min(n/(n/i),m/(m/i));
            ans+=(sum[j]-sum[i-1])*(n/i)*(m/i);
        }
        printf("%lld
",ans);
    }
    return 0;
}