NKOJ 3441 Lucas的数论（杜教筛）

P3441【HN Training 2015 Round5】lucas的数论

问题描述

数据范围

对于100%的数据:n<=1000000000

---

直接推式子，用到一个公式，这个公式也比较显然，就是根据定义直接得到，注意到不互质的数对乘积也一定会被乘积相同的一个互质数对算到。

$$Ans=\sum_{i=1}^{N}\sum_{j=1}^{N}\tau (i\times j)，已知\tau (i \times j)=\sum_{x|i}\sum_{y|j}1[gcd(x,y)=1]$$
$$Ans=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{x|i}\sum_{y|j}\sum_{p|gcd(x,y)}\mu(p)=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{p|gcd(i,j)}\tau(\frac{i}{p})\tau(\frac{j}{p})\mu(p)$$
$$Ans=\sum_{p=1}^{N}\mu(p)\sum_{i=1}^{\lfloor{\frac{N}{p}}\rfloor}\tau(i)\sum_{j=1}^{\lfloor{\frac{N}{p}}\rfloor}\tau(j)=\sum_{p=1}^{N}\mu(p)[\sum_{i=1}^{\lfloor{\frac{N}{p}}\rfloor}\tau(i)]^2$$

注意到后面是取整的形式，如果知道$\tau(i)$和$\mu(i)$的前缀和，那么可以分块在$O(\sqrt{n})$内求解
考虑如何快速求出$\tau(i)$的前缀和与$\mu(i)$的前缀和，考虑杜教筛。

对于$\mu(i)$，注意到$\sum_{d|n}\mu(d)=[n==1]$，那么有$\sum_{i=1}^{n}\sum_{d|i}\mu(d)=1$，于是$\sum_{i=1}^{n}\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}\mu(d)=1$

$$\sum_{i=1}^{n}\mu(d)=1-\sum_{i=2}^{n}\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}\mu(d)，令\sum_{i=1}^{n}\mu(i)=M(n)$$
$$那么M(n)=1-\sum_{i=2}^{n}M(\lfloor\frac{n}{i}\rfloor)$$
上式显然可以分块迭代处理，预处理一部分后做到 $O(n^{\frac{2}{3}})$，从狄利克雷卷积的角度就是 $\mu*I=e$

对于$\tau(i)$，注意到$\tau(i)=\sum_{d|i}1$，即$\tau=I*I$，那么$\mu*\tau=I*e=I$，即$\sum_{d|n}\mu(\frac{n}{d})\tau(d)=1$，同样有$\sum_{i=1}^{n}\sum_{d|i}\mu(\frac{n}{d})\tau(d)=\sum_{i=1}^{n}\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}\mu(i)\tau(d)=1$，于是

$$\sum_{i=1}^{n}\tau(i)=1-\sum_{i=2}^{n}\mu(i)\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}\tau(d)，令T(n)=\sum_{i=1}^{n}\tau(i)$$
$$那么T(n)=1-\sum_{i=2}^{n}\mu(i)T(\lfloor\frac{n}{i}\rfloor)$$
上式显然还是分块迭代处理，但是需要用到 $\mu(i)$的前缀和，复杂度不好说，但还是比较快的。

关于预处理一部分$tau(i)$，线性筛的时候需要增加两个数组，一个记录$i$的最小质因子，另一个记录最小质因子的指数，由于线性筛每次筛掉一个数一定是用他的最小质因数，因此可以方便的转移，具体可以看代码。

事实上，求$\tau(i)$的前缀和有更快的方法，因为$\sum_{i=1}^{n}\tau(i)=\sum_{i=1}^{n}\sum_{d|i}1=\sum_{i=1}^{n}\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}1=\sum_{i=1}^{n}\lfloor\frac{n}{i}\rfloor$，直接分块+预处理可以做到$O(\sqrt{n})$，最终复杂度$O(n^{\frac{3}{4}})$，复杂度并不会证，看看就好。

---

代码（杜教筛求$\tau(i)$）：

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#define N 12345678
#define ll long long
using namespace std;
const ll mod=1e9+7;
map<ll,ll>Mu,G;
ll n,K,tot,P[N],g[N],mu[N],pc[N],pd[N];
void EU()
{
    ll i,j;mu[1]=g[1]=1;
    for(i=2;i<=K;i++)
    {
        if(!g[i])P[++tot]=i,mu[i]=-1,g[i]=2,pc[i]=1,pd[i]=tot;
        for(j=1;j<=tot&&i*P[j]<=K;j++)
        {
            if(i%P[j])
            {
                pc[i*P[j]]=1;
                pd[i*P[j]]=j;
                mu[i*P[j]]=-mu[i];
                g[i*P[j]]=g[i]<<1;
            }
            else
            {
                pc[i*P[j]]=j==pd[i]?pc[i]+1:1;
                pd[i*P[j]]=j;
                g[i*P[j]]=j==pd[i]?(g[i]/(pc[i]+1)*(pc[i]+2)):(g[i]<<1);
            }
        }
    }
    for(i=2;i<=K;i++)g[i]+=g[i-1],mu[i]+=mu[i-1],g[i]%=mod;
}
ll Gmu(ll x)
{
    if(x<=K)return mu[x];
    if(Mu[x])return Mu[x];
    ll i,j,ans=1;
    for(i=2;i<=x;i=j+1)
    {
        j=x/(x/i);
        ans-=(j-i+1)*Gmu(x/i);
    }
    return Mu[x]=ans;
}
ll Gg(ll x)
{
    if(x<=K)return g[x];
    if(G[x])return G[x];
    ll i,j,ans=x;
    for(i=2;i<=x;i=j+1)
    {
        j=x/(x/i);
        ans-=(Gmu(j)-Gmu(i-1))*Gg(x/i)%mod;ans%=mod;
    }
    return G[x]=ans;
}
int main()
{
    ll ans=0,i,j;
    scanf("%lld",&n);
    K=pow(n,0.66);EU();
    for(i=1;i<=n;i=j+1)
    {
        j=n/(n/i);
        ans+=(Gmu(j)-Gmu(i-1))*Gg(n/i)%mod*Gg(n/i)%mod;ans%=mod;
    }
    printf("%lld",(ans+mod)%mod);
}

另附$O(\sqrt{n})$求$\tau$做法，实测并快不了多少，但如果是算单个应该快很多。

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#define N 12345678
#define ll long long
using namespace std;
const ll mod=1e9+7;
map<ll,ll>Mu,G;
ll n,K,tot,P[N],g[N],mu[N],pc[N],pd[N];
void EU()
{
    ll i,j;mu[1]=g[1]=1;
    for(i=2;i<=K;i++)
    {
        if(!g[i])P[++tot]=i,mu[i]=-1,g[i]=2,pc[i]=1,pd[i]=tot;
        for(j=1;j<=tot&&i*P[j]<=K;j++)
        {
            if(i%P[j])
            {
                pc[i*P[j]]=1;
                pd[i*P[j]]=j;
                mu[i*P[j]]=-mu[i];
                g[i*P[j]]=g[i]<<1;
            }
            else
            {
                pc[i*P[j]]=j==pd[i]?pc[i]+1:1;
                pd[i*P[j]]=j;
                g[i*P[j]]=j==pd[i]?(g[i]/(pc[i]+1)*(pc[i]+2)):(g[i]<<1);
            }
        }
    }
    for(i=2;i<=K;i++)g[i]+=g[i-1],mu[i]+=mu[i-1],g[i]%=mod;
}
ll Gmu(ll x)
{
    if(x<=K)return mu[x];
    if(Mu[x])return Mu[x];
    ll i,j,ans=1;
    for(i=2;i<=x;i=j+1)
    {
        j=x/(x/i);
        ans-=(j-i+1)*Gmu(x/i);
    }
    return Mu[x]=ans;
}
ll Gg(ll x)
{
    if(x<=K)return g[x];
    if(G[x])return G[x];
    ll ans=0,i,j;
    for(i=1;i<=x;i=j+1)
    {
        j=x/(x/i);
        ans+=(j-i+1)*(x/i)%mod;ans%=mod;
    }
    return G[x]=ans;
}
int main()
{
    ll ans=0,i,j;
    scanf("%lld",&n);
    K=pow(n,0.66);EU();
    for(i=1;i<=n;i=j+1)
    {
        j=n/(n/i);
        ans+=(Gmu(j)-Gmu(i-1))*Gg(n/i)%mod*Gg(n/i)%mod;ans%=mod;
    }
    printf("%lld",(ans+mod)%mod);
}