【BZOJ 2671】 2671: Calc (数论,莫比乌斯反演)
2671: Calc
Time Limit: 10 Sec Memory Limit: 128 MB
Submit: 303 Solved: 157Description
给出N,统计满足下面条件的数对(a,b)的个数:
1.1<=a<b<=N
2.a+b整除a*bInput
一行一个数N
Output
一行一个数表示答案
Sample Input
15
Sample Output
4
HINT
数据规模和约定
Test N Test N
1 <=10 11 <=5*10^7
2 <=50 12 <=10^8
3 <=10^3 13 <=2*10^8
4 <=5*10^3 14 <=3*10^8
5 <=2*10^4 15 <=5*10^8
6 <=2*10^5 16 <=10^9
7 <=2*10^6 17 <=10^9
8 <=10^7 18 <=2^31-1
9 <=2*10^7 19 <=2^31-1
10 <=3*10^7 20 <=2^31-1Source
【分析】
这题的复杂度还挺迷人的。
然后$\sqrt n$也没发现,以为筛$\mu$都要$O(n)$,什么杜教筛的幸好不会。。
首先分析$(a+b)|(a*b) → (a/g+b/g)|(a/g*b/g*g) →(a/g+b/g)|g$
那就是互质的$a',b'$ 找他们的公倍数$g$就行了。
写正常一点就是$$\sum_{j=1}^{N}\sum_{i=1}^{j-1}\dfrac{n}{j*(i+j)} [gcd(i,j)==1]$$
到了这里,我就傻眼了,其实嘛。。。j并不会到$n$,只是到$\sqrt{n}$
$$\sum_{j=1}^{\sqrt{n}}\sum_{i=1}^{j-1}\dfrac{n}{j*(i+j)} [gcd(i,j)==1]$$
然后我又傻眼了,复杂度迷人的东西啊会把我脑子弄得很乱的。
直接枚举j,然后i那里分块,然后就是求[l,r]里面和j互质的数的个数。
差分,先求[1,r]里面的,就是$\sum_{x=1}^{r}1[gcd(x,j)==1]$
即$\sum_{d=1}^{r}\mu(d)*(r/d)$
最后就是$\sum_{d=1}^{r}\mu(d)*(r/d-(l-1)/d)$
枚举约数在前面$\sqrt{j}$枚举去了。。
真的是暴力出奇迹了。。。
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cstring> 4 #include<iostream> 5 #include<algorithm> 6 #include<cmath> 7 using namespace std; 8 #define Maxn 50010 9 #define LL long long 10 11 bool vis[Maxn]; 12 int pri[Maxn],pl,mu[Maxn]; 13 14 int mymin(int x,int y) {return x<y?x:y;} 15 16 void init() 17 { 18 memset(vis,0,sizeof(vis)); 19 pl=0;mu[1]=1; 20 for(int i=2;i<=Maxn;i++) 21 { 22 if(!vis[i]) pri[++pl]=i,vis[i]=1,mu[i]=-1; 23 for(int j=1;j<=pl;j++) 24 { 25 if(i*pri[j]>Maxn) break; 26 vis[i*pri[j]]=1; 27 if(i%pri[j]==0) {mu[i*pri[j]]=0;break;} 28 mu[i*pri[j]]=-mu[i]; 29 } 30 } 31 } 32 33 int sta[Maxn],sl; 34 35 void div(int x) 36 { 37 sl=0; 38 int i; 39 for(i=1;i*i<x;i++) 40 { 41 if(x%i==0) sta[++sl]=i,sta[++sl]=x/i; 42 } 43 if(i*i==x) sta[++sl]=i; 44 } 45 46 int gcd(int a,int b) 47 { 48 if(b==0) return a; 49 return gcd(b,a%b); 50 } 51 52 int main() 53 { 54 int n; 55 scanf("%d",&n); 56 init(); 57 LL ans=0; 58 int sq=(int)(sqrt((double)n)); 59 for(int i=1;i<=sq;i++) 60 { 61 div(i); 62 for(int j=1;j<i;) 63 { 64 int x=n/i/(i+j),r;if(x==0) break; 65 r=mymin(i-1,n/x/i-i); 66 for(int k=1;k<=sl;k++) 67 { 68 ans+=mu[sta[k]]*(r/sta[k]-(j-1)/sta[k])*x; 69 } 70 j=r+1; 71 } 72 } 73 printf("%lld\n",ans); 74 return 0; 75 }
2017-04-06 15:50:26