[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 }

 

转载于:https://www.cnblogs.com/skylee03/p/8467489.html