题目大意：
定义$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）
1 #include<cstdio> 2 #include<cctype> 3 #include<algorithm> 4 typedef long long int64; 5 inline int getint() { 6 register char ch; 7 while(!isdigit(ch=getchar())); 8 register int x=ch^'0'; 9 while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0'); 10 return x; 11 } 12 const int N=10000001,M=N; 13 bool vis[N]; 14 int p[M],a[N],b[N]; 15 int64 g[N],sum[N]; 16 inline void sieve() { 17 for(register int i=2;i<N;i++) { 18 if(!vis[i]) { 19 p[++p[0]]=i; 20 a[i]=g[i]=1; 21 b[i]=i; 22 } 23 for(register int j=1;j<=p[0]&&i*p[j]<N;j++) { 24 vis[i*p[j]]=true; 25 if(i%p[j]==0) { 26 a[i*p[j]]=a[i]+1; 27 b[i*p[j]]=b[i]*p[j]; 28 if(i==b[i]) { 29 g[i*p[j]]=1; 30 } else { 31 g[i*p[j]]=(a[i*p[j]]==a[i/b[i]])?-g[i/b[i]]:0; 32 } 33 break; 34 } else { 35 a[i*p[j]]=1; 36 b[i*p[j]]=p[j]; 37 g[i*p[j]]=(a[i]==1)?-g[i]:0; 38 } 39 } 40 } 41 for(register int i=1;i<N;i++) { 42 sum[i]=sum[i-1]+g[i]; 43 } 44 } 45 int main() { 46 sieve(); 47 for(register int T=getint();T;T--) { 48 const int n=getint(),m=getint(),lim=std::min(n,m); 49 int64 ans=0; 50 for(register int i=1,j;i<=lim;i=j+1) { 51 j=std::min(n/(n/i),m/(m/i)); 52 ans+=(sum[j]-sum[i-1])*(n/i)*(m/i); 53 } 54 printf("%lld ",ans); 55 } 56 return 0; 57 }