bzoj4176 Lucas的数论 （杜教筛 +莫比乌斯反演）

bzoj4176 Lucas的数论

原题地址：http://www.lydsy.com/JudgeOnline/problem.php?id=4176

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

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

其中 表示ij的约数个数。
他发现答案有点大，只需要输出模1000000007的值。

数据范围
n <= 10^9

题解：
由 d(nm)=∑i∣n∑j∣m[gcd(i,j)==1] 证明见文末
原式
=∑ni=1∑mj=1∑x∣i∑y∣j[gcd(x,y)==1]
=∑ni=1∑mj=1∑x∣i∑y∣j∑t∣x且t∣yμ(t)
=∑nt=1μ(t)∑t∣x∑x∣i∑t∣y∑y∣j1 (i,x,y,j<=n)
观察到后面两个式子的形式几乎完全一样，
=∑nt=1μ(t)(∑t∣x∑x∣i1)2 (i,x<=n)
=∑nt=1μ(t)(∑ntxt=1⌊nx⌋)2
用i替换x/t，得：
=∑nt=1μ(t)(∑nti=1⌊⌊nt⌋i⌋)2

定义，f(n)=∑ni=1⌊ni⌋
那么，
=∑nt=1μ(t) f(⌊nt⌋)2
前面对于μ的前缀和，有：

定义梅滕斯函数M(n)=∑ni=1μ(i)，使用[n=1]=∑d|nμ(d)

1=∑i=1n[i=1]=∑i=1n∑d|iμ(d)=∑i=1n∑d=1⌊ni⌋μ(d)=∑i=1nM(⌊ni⌋)

因此M(n)=1−∑ni=2M(⌊ni⌋)，问题可在O(n23)时间复杂度下解决。

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 ——摘自浅谈一类积性函数的前缀和
预处理 i<=n34的μ前缀，更大的用上述用与e()的卷积推出来的东西，分块+缩小范围。
因为我们最后for的时候是按(n/i)分块一段一段地for，所以访问到的mu前缀和只会是∑(μ(ni))
所以要预先处理n/i>S的 ∑(μ(ni))
因为 n/i>S，所以i<=S，所以数组下标是按i保存。
对于后面f(⌊nt⌋)2的部分，分块计算f(n)=∑ni=1⌊ni⌋
据PoPoQQQ大佬的博客，时间复杂度是O(n1−−√)+O(n2−−√)+...+O(nn√−−−√)=O(n34)
于是总复杂度O(n34)

代码：

#include<cstdio>#include<iostream>#include<cmath>#include<cstring>#include<algorithm>#define LL long longusing namespace std;const int mod=1000000007;const int N=10000007;int n,S,mu[N],sum[N],p[N],ptot=0;bool is[N];void shai(){    memset(is,0,sizeof(is)); is[1]=1; mu[1]=1;     for(int i=2;i<=S;i++)    {        if(!is[i]) {p[++ptot]=i;mu[i]=-1;}        for(int j=1;j<=ptot;j++)        {            int pp=p[j];            if(1LL*pp*i>=1LL*N) break;             int x=pp*i; is[x]=1; mu[x]=mu[i]*mu[pp];            if(i%pp==0) {mu[x]=0; break;}        }    }    mu[0]=0;    for(int i=1;i<=S;i++) mu[i]=(mu[i-1]+mu[i]+mod)%mod;}int cal(int x){    int ans=0;    for(int i=1,ed;i<=x;i=ed+1)    {        ed=(x/(x/i));        ans=(ans+(1LL*(ed-i+1)*(x/i))%mod)%mod;    }    return ans;}int getmu(int x){    if(x<=S) return mu[x];    else return sum[n/x];}void init(){    shai();    int top=1; while(n/top>S) top++;    for(int j=top;j>=1;j--)    {        int now=n/j; sum[j]=1;        for(int i=2,ed;i<=now;i=ed+1)        {            ed=now/(now/i);             sum[j]=(sum[j]-(1LL*(ed-i+1)*getmu(now/i))%mod+mod)%mod;        }    }   }int main(){    scanf("%d",&n); S=ceil(pow((double)n,0.75));    init();     int ans=0;    for(int i=1,ed;i<=n;i=ed+1)    {        ed=n/(n/i); int tmp=cal((n/i));        int ret=(getmu(ed)-getmu(i-1)+mod)%mod;        ret=(1LL*ret*tmp)%mod; ret=(1LL*ret*tmp)%mod;         ans=(ans+ret)%mod;    }    printf("%d\n",ans);    return 0;}

---

补充：

1、d()为约数个数函数
d(nm)=∑i∣n∑j∣m[gcd(i,j)==1]
证明如下：
n,m的任意约数都可表示为i∗j  (i∣n且j∣m)的形式
即 i∗mj  (i∣n且j∣m)
上式即是枚举i,j

为什么必须[gcd(i,j)==1] ?

若i,j都有因子p，
i∗mj等价于ip∗mjp
会算重。

2、⌊⌊ab⌋c⌋=⌊abc⌋