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)=\sum_{i\mid n}\sum_{j\mid m}[gcd(i,j)==1]$ $_{证明见文末}$
原式
$=\sum_{i=1}^n\sum_{j=1}^m\sum_{x\mid i}\sum_{y\mid j}[gcd(x,y)==1]$
$=\sum_{i=1}^n\sum_{j=1}^m\sum_{x\mid i}\sum_{y\mid j}\sum_{t\mid x且 t\mid y}\mu(t)$
$=\sum_{t=1}^n\mu(t)\sum_{t\mid x}\sum_{x\mid i}\sum_{t\mid y}\sum_{y\mid j}1$ (i,x,y,j<=n)
观察到后面两个式子的形式几乎完全一样，
$=\sum_{t=1}^n\mu(t)(\sum_{t\mid x}\sum_{x\mid i}1)^2$ (i,x<=n)
$=\sum_{t=1}^n\mu(t)(\sum_{\frac{x}{t}=1}^{\frac{n}{t}}⌊\frac{n}{x}⌋)^2$
用i替换x/t，得：
$=\sum_{t=1}^n\mu(t)(\sum_{i=1}^{\frac{n}{t}}⌊\frac{⌊\frac{n}{t}⌋}{i}⌋)^2$

定义，$f(n)=\sum_{i=1}^n⌊\frac{n}{i}⌋$
那么，
$=\sum_{t=1}^n\mu(t) \ f(⌊\frac{n}{t}⌋)^2$
前面对于$\mu$的前缀和，有：

定义梅滕斯函数$M(n)=\sum_{i=1}^{n}{\mu(i)}$，使用$[n=1]=\sum_{d|n}{\mu(d)}$

$$1=\sum_{i=1}^{n}{[i=1]}=\sum_{i=1}^{n}\sum_{d|i}{\mu(d)}=\sum_{i=1}^{n}\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}{\mu(d)}=\sum_{i=1}^{n}{M(\lfloor\frac{n}{i}\rfloor)}$$
因此 $M(n)=1-\sum_{i=2}^{n}{M(\lfloor\frac{n}{i}\rfloor)}$，问题可在 $O(n^\frac{2}{3})$时间复杂度下解决。

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

代码：

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#define LL long long
using 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)=\sum_{i\mid n}\sum_{j\mid m}[gcd(i,j)==1]$
证明如下：
n,m的任意约数都可表示为$i*j$ $\ (i \mid n且j \mid m)$的形式
即 $i*\frac mj$ $\ (i \mid n且j \mid m)$
上式即是枚举$i,j$

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

若$i,j$都有因子$p$，
$i*\frac mj等价于\frac ip*\frac m{\frac jp}$
会算重。

2、$⌊\frac{⌊\frac ab⌋}c⌋=⌊\frac{a}{bc}⌋$