【杜教筛学习笔记】

杜教筛

关于一些前置芝士可以看我上一个博客【数学笔记】
我们构造\(f,g\)为积性函数，令\(S(n) = \sum\limits_{i = 1}^n f(i)\)
\(\sum\limits_{i = 1}^{n}(f * g)(i) = \sum\limits_{i = 1}^n\sum\limits_{d|n}(f(d) * g(n/d)) = \sum\limits_{i = 1} ^ ng(i) * \sum\limits_{j = 1} ^ {\lfloor n / i \rfloor}f(j) = \sum\limits_{i = 1} ^ n g(i) * S(n / i)\)
移项整理
\(g(1) * S(n) = \sum\limits_{i = 1}^n(f * g)(i) - \sum\limits_{i = 2}^ng(i) * S(n/i)\)
然后杜教筛的形式就出来了
\(S(n) = \frac{\sum\limits_{i = 1}^n(f * g)(i) - \sum\limits_{i = 2}^ng(i) * S(n/i)}{g(1)}\)
利用\(\mu * I = \mu(i)\)和\(\varphi * I = Id(n)\)把\(g = I\)
那么我们可以求出关于\(\mu\ \varphi\)的杜教筛公式
\(S(n) = 1 - \sum\limits_{i = 2} ^ nS(n / i)\)
\(S(n) = \frac{(1 + x) * x}{2} - \sum\limits_{i = 2}^nS(n / i)\)
直接写代码的话\(O(n ^ {\frac{3}{4}})\)
可以考虑线筛筛出前面的一些值，\(当线筛到n ^ {\frac{2}{3}}\)时，杜教筛的复杂度是\(O(n ^ {\frac{2}{3}})\)
代码

#include <cstdio>
#include <map>
#define N 7000001
#define ll long long
using namespace std;
int t;
ll n,p[N],mo[N],phi[N];//phi:欧拉函数，mo:莫比乌斯函数
map<ll,ll>f1,f2;
bool used[N];
ll G(ll x)
{
    return x*(x+1)/2;
}
void euler()
{
    mo[1]=1;
    phi[1]=1;
    for(ll i=2;i<N;i++)
    {
        if(used[i]==0)
        {
            p[++p[0]]=i;
            used[i]=1;
            mo[i]=-1;
            phi[i]=i-1;
        }
        for(ll j=1;j<=p[0];j++)
        {
            if(i*p[j]>N)
                break;
            used[i*p[j]]=1;
            if(i%p[j]==0)
            {
                mo[i*p[j]]=0;
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            mo[i*p[j]]=-mo[i];
            phi[i*p[j]]=phi[i]*(p[j]-1);
        }
    }
    for(int i=1;i<N;i++)
    {
        phi[i]=phi[i-1]+phi[i];//前缀和
        mo[i]=mo[i-1]+mo[i];//前缀和
    }
}
void ask(ll x,ll &ans1,ll &ans2)
{
    if(x<N)
    {
        ans1=phi[x],ans2=mo[x];
        return;
    }
    if(f1[x]!=0)
    {
        ans1=f1[x],ans2=f2[x];
        return;
    }
    ans1=G(x),ans2=1;//初始值
    for(ll i=2,j;i<=x;i=j+1)
    {
        j=x/(x/i);
        ll a,b;
        ask(x/i,a,b);//递推
        ans1-=a*(j-i+1);
        ans2-=b*(j-i+1);
    }
    f1[x]=ans1;
    f2[x]=ans2;
    return ;
}
int main()
{
    scanf("%d",&t);
    euler();
    while(t--)
    {
        scanf("%lld",&n);
        ll a,b;
        ask(n,a,b);
        printf("%lld %lld\n",a,b);
    }

    return 0;
}