模板 - 杜教筛

取模版本:

```cpp #include using namespace std; typedef long long ll;

const int mod=1e9+7;
const int inv2=(mod+1)>>1;

const int MAXN=5e6;

unordered_map<ll,ll> Smu;
unordered_map<ll,ll> Sphi;

ll sphi[MAXN+1];
ll smu[MAXN+1];
int pri[MAXN+1];
int &pritop=pri[0];

void sieve(int n=MAXN) {
sphi[1]=1;
smu[1]=1;
for(int i=2; i<=n; i++) {
if(!pri[i]) {
pri[++pritop]=i;
sphi[i]=i-1;
smu[i]=-1;
}
for(int j=1; j<=pritop; j++) {
int &p=pri[j];
int t=ip;
if(t>n)
break;
pri[t]=1;
if(i%p) {
sphi[t]=sphi[i]
sphi[p];
if(sphi[t]>=mod)
sphi[t]%=mod;
smu[t]=-smu[i];
} else {
sphi[t]=sphi[i]*p;
if(sphi[t]>=mod)
sphi[t]%=mod;
smu[t]=0;
break;
}
}
}
for(int i=2; i<=n; i++) {
sphi[i]+=sphi[i-1];
if(sphi[i]>=mod)
sphi[i]%=mod;
smu[i]+=smu[i-1];
if(smu[i]<0)
smu[i]+=mod;
}
}

inline ll s1(ll n) {
if(n>=mod)
n%=mod;
ll tmp=n(n+1);
if(tmp>=mod)
tmp%=mod;
tmp
=inv2;
if(tmp>=mod)
tmp%=mod;
return tmp;
}

inline ll Mu(ll n) {
if(n<=MAXN)
return smu[n];
if(Smu.count(n))
return Smu[n];
ll ret=1;
for(ll l=2,r; l<=n; l=r+1) {
ll t=n/l;
r=n/t;
ll tmp=r-(l-1);
if(tmp<0)
tmp+=mod;
tmp*=Mu(t);
if(tmp>=mod)
tmp%=mod;
ret-=tmp;
if(ret<0)
ret+=mod;
}
return Smu[n]=ret;
}

inline ll Phi(ll n) {
if(n<=MAXN)
return sphi[n];
if(Sphi.count(n))
return Sphi[n];
ll ret=s1(n);
for(ll l=2,r; l<=n; l=r+1) {
ll t=n/l;
r=n/t;
ll tmp=r-(l-1);
tmp*=Phi(t);
if(tmp>=mod)
tmp%=mod;
ret-=tmp;
if(ret<0)
ret+=mod;
}
return Sphi[n]=ret;
}

int main() {
sieve();
ll n;
scanf("%lld",&n);
printf("%lld\n",Phi(n));
}

</details>

---

2019/6/2更新,杜教筛可以省去哈希表的复杂度和空间。当n不改变的时候,可以发现我们关心的都是n/i,大于sqrt(n)的是不超过sqrt(n)个。
可以换一种记忆方式省去哈希表的复杂度。(虽然实际上好像并没多快,反而还慢了,还巨难写?)但是当n改变时,需要清除多余的记忆。
<details>
```cpp
#include<bits/stdc++.h>
#define ll long long
using namespace std;

const int INF=0x3f3f3f3f;

//N为n^(2/3)最快
int n;

const int N23=pow(1.0*INT_MAX,2.0/3.0);
const int N12=pow(1.0*INT_MAX,1.0/2.0);
const int MAXN=N23;

ll sump[N23+N12+5];
int sumu[N23+N12+5];
int pri[MAXN+5],pritop;
bool notpri[MAXN+5];
//pritop从1开始计数

void sieve(int n=MAXN) {
    //对于欧拉函数,-1是特殊值,对于莫比乌斯函数,INF是特殊值
    notpri[0]=1,sump[0]=0,sumu[0]=0;
    notpri[1]=1,sump[1]=1,sumu[1]=1;
    for(int i=2; i<=n; i++) {
        if(!notpri[i])
            pri[++pritop]=i,sump[i]=i-1,sumu[i]=-1;
        for(int j=1; j<=pritop&&i*pri[j]<=n; j++) {
            notpri[i*pri[j]]=1;
            //略有不同
            if(i%pri[j])
                sump[i*pri[j]]=sump[i]*sump[pri[j]],sumu[i*pri[j]]=-sumu[i];
            else {
                sump[i*pri[j]]=sump[i]*pri[j],sumu[i*pri[j]]=0;;
                break;
            }
        }
    }

    for(int i=2; i<=n; i++) {
        sump[i]+=sump[i-1];
        sumu[i]+=sumu[i-1];
    }

    for(int i=MAXN+1;i<=MAXN+N12;i++){
        sump[i]=-1;
        sumu[i]=INF;
    }
}

inline int idx(int x) {
    return (x<=MAXN)?x:(n/x+MAXN);
}

//杜教筛欧拉函数
inline ll GetSphi(int n) {
    int t=idx(n);
    if(sump[t]!=-1)
        return sump[t];
    //小心溢出
    ll ret=((1ll+n)*n)/2;
    //f*g=id的前缀和
    for(int l=2,r; r<2147483647&&l<=n; l=r+1) {
        r=n/(n/l);
        ret-=(r-l+1)*GetSphi(n/l);
        //同上,因为两个的g都是I
    }
    return sump[t]=ret; // 记忆化
}

//杜教筛莫比乌斯函数
inline int GetSumu(int n) {
    int t=idx(n);
    if(sumu[t]!=INF)
        return sumu[t];
    int ret=1; // 单位元的前缀和就是 1
    for(int l=2,r;r<2147483647&&l<=n;l=r+1) {
        r=n/(n/l);
        ret-=(r-l+1)*GetSumu(n/l);
        //(r-l+1)就是I在[l,r]的和
    }
    return sumu[t]=ret;// 记忆化
}

int main() {
#ifdef Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    sieve();

    int t;
    scanf("%d",&t);
    while(t--) {
        for(int i=MAXN+1;i<=MAXN+N12;i++){
            sump[i]=-1;
            sumu[i]=INF;
        }
        scanf("%d",&n);
        printf("%lld %d\n",GetSphi(n),GetSumu(n));
    }
}

---

因为不会杜教筛(还有质数间隔),被一万六学院调教了。
补充一波杜教筛。

#include<bits/stdc++.h>
#define ll long long
using namespace std;

//N为n^(2/3)最快
int n;
const int MAXN=5e6;

unordered_map<int,int> Smu;
unordered_map<int,ll> Sphi;

ll sump[MAXN+5];
int sumu[MAXN+5];
int pri[MAXN+5],pritop;
bool notpri[MAXN+5];
//pritop从1开始计数

void sieve(int n) {
    notpri[1]=sump[1]=sumu[1]=1;
    for(int i=2; i<=n; i++) {
        if(!notpri[i])
            pri[++pritop]=i,sump[i]=i-1,sumu[i]=-1;
        for(int j=1; j<=pritop&&i*pri[j]<=n; j++) {
            notpri[i*pri[j]]=1;
            //略有不同
            if(i%pri[j])
                sump[i*pri[j]]=sump[i]*sump[pri[j]],sumu[i*pri[j]]=-sumu[i];
            else {
                sump[i*pri[j]]=sump[i]*pri[j],sumu[i*pri[j]]=0;;
                break;
            }
        }
    }
    for(int i=2; i<=n; i++) {
        sump[i]+=sump[i-1];
        sumu[i]+=sumu[i-1];
    }
}

//杜教筛莫比乌斯函数
inline ll GetSumu(int n) {
    if(n <= MAXN)
        return sumu[n]; // sumu是提前筛好的前缀和
    if(Smu.count(n))
        return Smu[n]; // 记忆化
    ll ret = 1ll; // 单位元的前缀和就是 1
    for(int l = 2, r; r<2147483647&&l <= n; l = r + 1) {
        r = n / (n / l);
        ret -= (r - l + 1) * GetSumu(n / l);
        // (r - l + 1) 就是 I 在 [l, r] 的和
    }
    return Smu[n] = ret; // 记忆化
}

//杜教筛欧拉函数
inline ll GetSphi(int n) {
    if(n <= MAXN)
        return sump[n]; // 提前筛好的
    if(Sphi.count(n)) return Sphi[n]; // 记忆化
    ll ret = n;
    if(ret%2==0)
        ret=(ret/2)* (1ll+n);
    else
        ret*=((1ll+n)/2);
    // f * g = id 的前缀和
    for(int l = 2, r; r<2147483647&&l <= n; l = r + 1) {
        r = n / (n / l);
        ret -= (r - l + 1) * GetSphi(n / l);
        // 同上,因为两个的 g 都是 I
    }
    return Sphi[n] = ret; // 记忆化
}

int main() {
    sieve(MAXN);

    int t;
    scanf("%d",&t);
    while(t--) {
        scanf("%d",&n);
        printf("%lld %lld\n",GetSphi(n),GetSumu(n));
    }
}

用来快速求出积性函数的前缀和(因为一般分块也是用前缀和,所以很多莫比乌斯反演可以用杜教筛优化到 \(O(n^{\frac{2}{3}})\) 预处理,然后整除分块 \(O(n^{\frac{1}{2}})\) 查询)。

前置知识

积性函数

积性函数:对于任意互质的整数 \(a,b\)\(f(ab)=f(a)f(b)\) 则称 \(f(x)\) 的数论函数。
完全积性函数:对于任意整数 \(a,b\)\(f(ab)=f(a)f(b)\) 的数论函数。
积性函数的积未必是积性函数,狄利克雷卷积才是。

常见的积性函数:\(\varphi,\mu,\sigma,d\)
常见的完全积性函数:\(\epsilon,I,id\)
这里特殊解释一下 \(\epsilon,I,id\) 分别是什么意思: \(\epsilon(n) = [n=1], I(n) = 1, id(n) = n\)

狄利克雷卷积

\(f, g\) 是两个数论函数,它们的狄利克雷卷积卷积是:\((f*g)(n) = \sum \limits _{d | n} f(d) g(\frac{n}{d})\)

性质:满足交换律,结合律

单位元: \(\epsilon\) (即 \(f*\epsilon=f\)

结合狄利克雷卷积得到的几个性质:

\(\mu * I = \epsilon\)
\(\varphi * I = id\)
\(\mu * id = \varphi\)

杜教筛

设现在要求积性函数 \(f\) 的前缀和, 设 \(\sum \limits_{i=1}^{n} f(i) = S(n)\)

再找一个积性函数 \(g\) ,则考虑它们的狄利克雷卷积的前缀和

\(\sum \limits_{i=1}^{n}(f*g)(i)\)
\(= \sum \limits_{i=1}^{n} \sum \limits_{d|i} f(d)g(\frac{i}{d})\)
\(= \sum \limits_{d=1}^{n} g(d)\sum\limits_{i=1}^{\lfloor \frac{n}{d}\rfloor } f(i)\)
\(= \sum \limits_{d=1}^{n} g(d) S(\lfloor \frac{n}{d} \rfloor)\)

则:
\(g(1)S(n)=\sum\limits_{i=1}^{n}(f*g)(i) - \sum \limits _{i=2}^{n} g(i) S(\lfloor \frac{n}{i} \rfloor)\)

//杜教筛莫比乌斯函数
inline ll GetSumu(int n) {
  if(n <= N) return sumu[n]; // sumu是提前筛好的前缀和
  if(Smu[n]) return Smu[n]; // 记忆化
  ll ret = 1ll; // 单位元的前缀和就是 1
  for(int l = 2, r; l <= n; l = r + 1) {
    r = n / (n / l); ret -= (r - l + 1) * GetSumu(n / l);
    // (r - l + 1) 就是 I 在 [l, r] 的和
  } return Smu[n] = ret; // 记忆化
}

//杜教筛欧拉函数
inline ll GetSphi(int n) {
  if(n <= N) return sump[n]; // 提前筛好的
  if(Sphi[n]) return Sphi[n]; // 记忆化
  ll ret = 1ll * n * (n + 1) / 2; // f * g = id 的前缀和
  for(int l = 2, r; l <= n; l = r + 1) {
    r = n / (n / l); ret -= (r - l + 1) * GetSphi(n / l);
    // 同上,因为两个的 g 都是 I
  } return Sphi[n] = ret; // 记忆化
}

注意一些常用的变换,例如 \(\sum\limits_{d|n}\mu(n)=[n==1]\) ,\(\sum\limits_{d|n}\varphi(n)=n\)


不用哈希表版本还带取模的杜教筛经过了欧拉函数前缀和的验证,跑得飞快。

#include<bits/stdc++.h>
#define ll long long
using namespace std;

const int INF=0x3f3f3f3f;
const int mod=1e9+7;
const int inv2=(mod+1)>>1;

//N为n^(2/3)最快
ll n;

const int MAXN23=pow(1e10,2.0/3.0);
const int MAXN12=pow(1e10,1.0/2.0);
const int MAXN=MAXN23;

ll sump[MAXN23+MAXN12+10];
int pri[MAXN+10],pritop;
bool notpri[MAXN+10];
//pritop从1开始计数

void sieve(int n=MAXN) {
    notpri[0]=1,sump[0]=0;
    notpri[1]=1,sump[1]=1;
    for(int i=2; i<=n; i++) {
        if(!notpri[i])
            pri[++pritop]=i,sump[i]=(i-1);
        for(int j=1,t; j<=pritop&&(t=i*pri[j])<=n; j++) {
            notpri[t]=1;
            //略有不同
            if(i%pri[j]){
                sump[t]=sump[i]*sump[pri[j]];
                if(sump[t]>=mod)
                    sump[t]%=mod;
            }
            else {
                sump[t]=sump[i]*pri[j];
                if(sump[t]>=mod)
                    sump[t]%=mod;
                break;
            }
        }
    }

    //cout<<"done!"<<endl;
    for(int i=2; i<=n; i++) {
        sump[i]+=sump[i-1];
        if(sump[i]>=mod)
            sump[i]-=mod;
    }

    //不用考虑边界,因为比MAXN还大的才会进入这个记忆化
    for(int i=MAXN+1;i<=MAXN+MAXN12;i++){
        sump[i]=-1;
    }
}

inline int idx(ll x) {
    return (x<=MAXN)?x:(n/x+MAXN);
}

inline ll s1(ll n){
    if(n>=mod)
        n%=mod;
    ll t=(n*(n+1)/2);
    if(t>=mod)
        t%=mod;
    return t;
}

//杜教筛欧拉函数
inline ll GetSphi(ll n) {
    int t=idx(n);
    if(sump[t]!=-1)
        return sump[t];
    //小心溢出
    ll ret=s1(n);
    //f*g=id的前缀和
    for(ll l=2,r;l<=n; l=r+1) {
        r=n/(n/l);
        ll tmp=(r-l+1)%mod*GetSphi(n/l)%mod;
        ret-=tmp;
        if(ret<0)
            ret+=mod;
        //同上,因为两个的g都是I
    }
    return sump[t]=ret; // 记忆化
}

int main() {
#ifdef Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    sieve();
    scanf("%lld",&n);
    printf("%lld",GetSphi(n));
}

posted @ 2019-04-22 00:50  韵意  阅读(210)  评论(0编辑  收藏  举报