杜教筛

杜教筛

杜教筛可以在非线性时间内求积性函数前缀和
设现在要求积性函数 \(f\) 的前缀和, 设 \(\sum_{i=1}^{n}f(i)=S(n)\)
找一个积性函数\(g\),考虑它们的狄利克雷卷积的前缀和

\[\sum_{i=1}^{n} (f*g)(i)=\sum_{i=1}^{n}\sum_{d|i}f(d)g(\frac{i}{d})
\]

先枚举 \(d\) 提出 \(g\)

\[\sum_{i=1}^{n} (f*g)(i)= \sum_{d=1}^{n}g(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)
\]

\(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\)替换为\(S(\lfloor\frac{n}{d}\rfloor)\)

\[\sum_{i=1}^{n} (f*g)(i)= \sum_{d=1}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)
\]

可以发现从\(1\)开始的前缀和减去从\(2\)开始的前缀和就是第一项,即

\[g(1)S(n)=\sum_{i=1}^{n}g(i)S(\lfloor\frac{n}{i}\rfloor)-\sum_{i=2}^{n}g(i)S(\lfloor\frac{n}{i}\rfloor)
\]

所以得到杜教筛的核心式子:

\[g(1)S(n)=\sum_{i=1}^{n}(f*g)(i)-\sum_{i=2}^{n}g(i)S(\lfloor\frac{n}{i}\rfloor)
\]

则如果可以找到一个合适的积性函数\(g\),使得可以快速算出\(\sum_{i=1}^{n}(f*g)(i)\)\(g\)的前缀和,便可以用数论分块递归地求解

代码框架如下

il int GetSum_f(cs int n){ // 算 f 前缀和的函数
    ri int as=sum_f_g(n); // sum_f_g(n)是算 f * g 的前缀和
    // 数论分块
    for(ri int l=2,r=0;l<=n;l=r+1){ // 注意从 2 开始
        r=(n/(n/l)); // g_sum 是 g 的前缀和,递归求解
        as-=(sum_g(r)-sum_g(l-1))*GetSum_f(n/l);
    }
    return as;
}

复杂度\(O(n^{\frac{3}{4}})\)

还可以进一步优化,即先线性筛出前 \(m\) 个答案,之后再用杜教筛,复杂度是\(O(\frac{n}{\sqrt{n}})\)
\(m=n^{\frac{2}{3}}\)时,复杂度为\(O(n^\frac{2}{3})\)

可以使用哈希表来存下已经求过的答案,第一篇题解说可以不用。

考虑到上面的求和过程中出现的都是 \(\lfloor\frac{n}{i}\rfloor\)。开一个大小为两倍 \(\sqrt{n}\) 的数组 \(f\) 记录答案。如果现在需要求出 GetSum_f(x) ,若 \(x\le\sqrt{n}\),返回 \(f[x]\) ,否则返回 \(f[sqrt(n)+n/i]\) 即可。

code

\(\sum_{i=1}^n\varphi(i)\)\(\sum_{i=1}^n\mu(i)\)

AC·code
#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
#define ll long long
using namespace std;
namespace Q{
    il int rd(){
        ri int x=0,f=1;ri char c=getchar();
        while(c<'0'||c>'9') f=(c=='-')?-1:1,c=getchar();
        while(c>='0'&&c<='9') x=x*10+(c^48),c=getchar();
        return x*f;
    }
    il void wt(int x){
        if(x<0) x=-x,putchar('-');
        if(x>=10) wt(x/10);
        return putchar(x%10+48),void();
    }
    il void wl(ll x){
        if(x<0) x=-x,putchar('-');
        if(x>=10) wl(x/10);
        return putchar(x%10+48),void();
    }
} 
namespace djs{
    cs int N=5e6;
    int ss[N],mu[N+5],cnt;
    ll fi[N+5];
    il void init(){
        mu[1]=fi[1]=1;
        for(ri int i=2;i<=N;++i){
            if(!fi[i]) ss[++cnt]=i,fi[i]=i-1,mu[i]=-1;
            for(ri int j=1;j<=cnt&&1ll*ss[j]*i<=N;++j){
                if(i%ss[j]){
                    fi[i*ss[j]]=fi[i]*fi[ss[j]];
                    mu[i*ss[j]]=-mu[i];
                }
                else{
                    fi[i*ss[j]]=fi[i]*ss[j];
                    break;
                }
            }
        }
        for(ri int i=2;i<=N;++i){
            fi[i]+=fi[i-1],mu[i]+=mu[i-1];
        }
        return;
    }
    unordered_map<int,int> Mu;
    unordered_map<int,ll> Fi;
    il ll varphi(int n){
        if(n<=N) return fi[n];
        else if(Fi.count(n)) return Fi[n];
        ri ll phi=(1ll+n)*n/2;
        for(ri unsigned int l=2,r=0;l<=n;l=r+1){
            r=n/(n/l),phi-=1ll*(r-l+1)*varphi(n/l);
        }
        return Fi[n]=phi;
    }
    il int get_mu(int n){
        if(n<=N) return mu[n];
        else if(Mu.count(n)) return Mu[n];
        ri int smu=1;
        for(ri unsigned int l=2,r=0;l<=n;l=r+1){
            r=n/(n/l),smu-=(r-l+1)*get_mu(n/l);
        }
        return Mu[n]=smu;
    }

} using namespace djs;

signed main(){
    init();
    for(ri int i=Q::rd(),n;i;--i){
        n=Q::rd();
        Q::wl(varphi(n)),putchar(' ');
        Q::wt(get_mu(n)),putchar('\n');
    }
    return 0;
}

原文链接: https://www.cnblogs.com/windymoon/p/17071163.html

欢迎关注

微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍

    杜教筛

原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/314458

非原创文章文中已经注明原地址,如有侵权,联系删除

关注公众号【高性能架构探索】,第一时间获取最新文章

转载文章受原作者版权保护。转载请注明原作者出处!

(0)
上一篇 2023年2月16日 下午1:22
下一篇 2023年2月16日 下午1:22

相关推荐