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