[BZOJ3309]DZY Loves Math
题目大意：
定义$f(n)$为$n$所含质因子的最大幂指数，共$q(q\leq10000)$组询问，对于给定的$n,m(n,m\leq10^7)$，求$\sum_{i=1}^n\sum_{j=1}^mf(\gcd(i,j))$。
思路：
$$
\begin{align*}
原式&=\sum_{d=1}^{\min(n,m)}f(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1]\\
&=\sum_{d=1}^{\min(n,m)}f(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{p|\gcd(i,j)}\mu(p)\\
&=\sum_{d=1}^{\min(n,m)}f(d)\sum_{p=1}^{\min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\mu(p)\sum_{i=1}^{\lfloor\frac{n}{dp}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{dp}\rfloor}1\\
&=\sum_{d=1}^{\min(n,m)}f(d)\sum_{p=1}^{\min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\mu(p)\lfloor\frac{n}{dp}\rfloor\lfloor\frac{m}{dp}\rfloor\\
&=\sum_{d=1}^{\min(n,m)}\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor\sum_{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)}\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor 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\n",ans); 55 } 56 return 0; 57 }
