【bzoj 4176】 Lucas的数论莫比乌斯反演（杜教筛）

Description

去年的Lucas非常喜欢数论题，但是一年以后的Lucas却不那么喜欢了。

在整理以前的试题时，发现了这样一道题目“求Sigma(f(i)),其中1<=i<=N”，其中表示i的约数个数。他现在长大了，题目也变难了。

求如下表达式的值：

一行一个整数ans，表示答案模1000000007的值。

Sample Input

Sample Output

HINT

对于100%的数据n <= 10^9。

题解：

　　解锁新技能：杜教筛。

　　再复习一下：

　　若$F(n)=sum_{i=1}^{n}f(i)，g(i)=sum_{j|i}f(i),G(n)=sum_{i=1}^{n}g(i)$，

　　则有：$G(n)=sum_{i=1}^{n}F(lfloorfrac{n}{i} floor)$。

　　即：$F(n)=G(n)-sum_{i=2}^{n}F(lfloorfrac{n}{i} floor)$

　　然后就可以上杜教筛了。

　　话归正题：
　　对于本题，其实和bzoj3994式子一样……

$sum_{i=1}^{n}sum_{j=1}^{n}d(ij)$

$=sum_{i=1}^{n}sum_{j=1}^{n}sum_{p=1}^{n^{2}}[p|ij]$

$=sum_{i=1}^{n}sum_{j=1}^{n}sum_{p=1}^{n^{2}}[frac{p}{(p,i)}|j]$

$=sum_{t=1}^{n}sum_{i=1}^{lfloorfrac{n}{t} floor}sum_{j=1}^{n}sum_{p=1}^{lfloorfrac{n^{2}}{t} floor}[p|j]*[(i,p)==1]$

$=sum_{t=1}^{n}sum_{i=1}^{lfloorfrac{n}{t} floor}sum_{p=1}^{lfloorfrac{n^{2}}{t} floor}sum_{d|(i,p)}mu(d)lfloor frac{n}{p} floor$

$=sum_{d=1}^{n}mu(d)sum_{t=1}^{lfloorfrac{n}{d} floor}sum_{i=1}^{lfloorfrac{n}{td} floor}sum_{p=1}^{lfloorfrac{n^{2}}{td} floor}lfloorfrac{n}{pd} floor$

$n>=pd,n^{2}>=tdn$

$ herefore Ans=sum_{d=1}^{n}mu(d)sum_{t=1}^{lfloorfrac{n}{d} floor}sum_{i=1}^{lfloorfrac{n}{td} floor}sum_{p=1}^{lfloorfrac{n}{td} floor}lfloorfrac{n}{pd} floor$

$Ans=sum_{d=1}^{n}mu(d)sum_{t=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{td} floorsum_{p=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{pd} floor$

$设f(i)=sum_{x=1}^{i}lfloorfrac{i}{x} floor$

$ herefore Ans=sum_{d=1}^{n}mu(d)f^2(lfloorfrac{n}{d} floor)$

　　然后分块+杜教筛即可。

代码：

 1 #include<cstdio>
 2 using namespace std;
 3 typedef long long ll;
 4 const ll N=(ll)3e6+1;
 5 const ll mod=1000000007ll;
 6 short miu[N];
 7 int prim[N/10],num;
 8 bool vis[N];
 9 int sum[N];
10 inline void init(){
11     miu[1]=sum[1]=1;
12     for(int i=2;i<N;i++){
13         if(!vis[i]){
14             prim[++num]=i;
15             miu[i]=-1;
16         }for(int j=1;j<=num&&prim[j]*i<N;j++){
17             vis[i*prim[j]]=1;
18             if(i%prim[j]==0){
19                 miu[i*prim[j]]=0;
20                 break;
21             }
22             else
23                 miu[i*prim[j]]=-miu[i];
24         }
25       //  printf("miu[%d]=%d
",i,miu[i]);
26     }
27     for(int i=2;i<N;i++)
28         sum[i]=sum[i-1]+miu[i];
29 }
30 struct edges{
31     ll v;int w;edges *last;
32 }edge[N/10],*head[76545];int cnt;
33 const int limit=76543;
34 inline void push(int u,ll v,int w){
35     edge[++cnt]=(edges){v,w,head[u]};head[u]=edge+cnt;
36 }
37 inline ll Get_sum(ll x){
38     if(x<N) return sum[x];
39     int t=x%limit;
40     for(edges *i=head[t];i;i=i->last)
41         if(i->v==x) return i->w;
42     ll ans=1;
43     for(ll i=2,pos;i<=x;i=pos+1){
44         pos=x/i;
45         pos=x/pos;
46         ans-=(pos-i+1)*Get_sum(x/pos);
47     }
48     push(t,x,ans);
49     return ans;
50 }
51 inline ll Get_F(ll x){
52     ll ans=0;
53     for(ll i=1,pos;i<=x;i=pos+1){
54         pos=x/i;
55         pos=x/pos;
56         ll t=x/i;
57         ans+=t*(pos-i+1);
58         ans%=mod;
59     }return ans;
60 }
61 inline ll solve(ll x){
62     ll ans=0;
63     for(ll i=1,pos;i<=x;i=pos+1){
64         pos=x/i;
65         pos=x/pos;
66         ll t=Get_F(x/i);
67         ans+=(Get_sum(pos)-Get_sum(i-1))%mod*t%mod*t%mod;
68         ans%=mod;
69     }
70     return (ans%mod+mod)%mod;
71 }
72 int main(){
73     init();
74     ll a;
75     scanf("%lld",&a);
76     printf("%lld
",solve(a));
77 }

【bzoj 4176】 Lucas的数论 莫比乌斯反演（杜教筛）

代码：

相关推荐

【bzoj 4176】 Lucas的数论莫比乌斯反演（杜教筛）