【BZOJ4176】Lucas的数论-杜教筛

求$$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}f(ij)$$,其中$f(x)$表示$x$的约数个数,$0\leq n\leq 10^9$,答案膜$10^9+7$

题解

首先有个妙不可言(被hjw污染了)的结论:$$f(nm)=\sum\limits_{i|n}\sum\limits_{j|m}[gcd(i,j)=1]$$

证明:咕

那么大力推一波式子:

$$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}f(ij)$$

$$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{a|i}\sum\limits_{b|j}[gcd(a,b)=1]$$

$$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{a|i}\sum\limits_{b|j}\sum\limits_{d|gcd(a,b)}\mu(d)$$

$$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{a|i}\sum\limits_{b|j}\sum\limits_{d|a\& d|b}\mu(d)$$

$$=\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{a=1}^{\lfloor\frac{n}{id}\rfloor}\sum\limits_{b=1}^{\lfloor\frac{n}{jd}\rfloor}\mu(d)$$

$$=\sum\limits_{d=1}^{n}\mu(d)(\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{id}\rfloor)$$

杜教筛+莫比乌斯反演解决

时间复杂度:$O(n^{\frac{2}{3}}logn+n^{\frac{3}{4}})$

代码:

 

 1 #include<iostream>
 2 #include<cstring>
 3 #include<cstdio>
 4 #include<cmath>
 5 #include<map>
 6 #define mod 1000000007
 7 using namespace std;
 8 typedef long long ll;
 9 ll n,pri=0,p[1000001],miu[1000001],pre[1000001],ans=0;
10 bool isp[1000001];
11 map<ll,ll>HASH;
12 void _(){
13     miu[1]=pre[1]=1;
14     for(int i=2;i<=1000000;i++){
15         if(!isp[i]){
16             p[++pri]=i;
17             miu[i]=-1;
18         }
19         for(int j=1;j<=pri&&i*p[j]<=1000000;j++){
20             isp[i*p[j]]=true;
21             if(i%p[j]==0){
22                 miu[i*p[j]]=0;
23                 break;
24             }
25             miu[i*p[j]]=-miu[i];
26         }
27     }
28     for(int i=2;i<=1000000;i++){
29         pre[i]=(pre[i-1]+miu[i]+mod)%mod;
30     }
31 }
32 ll work1(ll x){
33     if(x<=1000000)return pre[x];
34     if(HASH.count(x))return HASH[x];
35     ll ret=1;
36     for(int i=2,j;i<=x;i=j+1){
37         j=x/(x/i);
38         ret=(ret-(j-i+1)*work1(x/i)%mod+mod)%mod;
39     }
40     HASH[x]=ret;
41     return ret;
42 }
43 ll work2(ll x){
44     ll ret=0;
45     for(int i=1,j;i<=x;i=j+1){
46         j=x/(x/i);
47         ret=(ret+(x/i)*(j-i+1))%mod;
48     }
49     return ret*ret%mod;
50 }
51 int main(){
52     _();
53     scanf("%lld",&n);
54     for(int i=1,j;i<=n;i=j+1){
55         j=n/(n/i);
56         ans=(ans+(work1(j)-work1(i-1)+mod)%mod*work2(n/i))%mod;
57     }
58     printf("%lld",ans);
59     return 0;
60 }

 

posted @ 2018-08-05 14:30  DCDCBigBig  阅读(224)  评论(0编辑  收藏  举报