[日常摸鱼]积性函数求和——杜教筛

因为博主比较菜所以可能一些地方写的有问题或者不清楚,以及我的废话好像有点多…

本文(大概)会不定期更新…一些东西的证明这里可能暂时没有

以及在这里先感谢下小伙伴 @MoebiusMeow 的帮助~ww

参考资料:

[1]浅谈一类积性函数的前缀和(skywalkert)

[2]杜教筛——省选前的学习1(_rqy)

[3]杜教筛(zzq)

[4]具体数学(Concrete Mathematics)


前置技能(一些定义)

约定$[p]$表示满足条件$p$时这个值为1,否则为0

欧拉函数:用$\phi(n)$表示小于$n$且和$n$互质的整数的个数,若$n=p_1^{q_1}*p_2^{q_2}*...*p_k^{q_k}$那么有一种求法是$\phi(n)=n*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})$

莫比乌斯函数:一种$\mu(n)$的定义是$\mu(1)=1$,对于无平方因子的$n=\prod_{i=1}^k p_i$,$\mu(n)=(-1)^k$,其他情况$\mu(n)=0$,另一种等价的定义是$\sum_{d|n}\mu(d)=[n=1]$

数论函数:若$f:Z^{+} \rightarrow C$,则称$f$为数论函数

积性函数:若一个数论函数$f(n)$对于所有$m_1 \bot m_2$有$f(m_1*m_2)=f(m_1)*f(m_2)$且$f(1)=1$那么称$f(n)$是积性函数,如果对于任意$m_1,m_2$都满足上面那个条件那么就称$f(n)$是完全积性函数,常见的比如莫比乌斯函数、欧拉函数、幂函数($id_k(n)=n^{k})$、除数函数($\sigma_k(n)=\sum_{d|n}d^k$))([1])

如果$f$是一个积性函数,把$n$写成$p_1*p_2*...*p_k$的形式,那么可以根据积性就得到$f(n)=\prod_{i=1}^k f(p_i)$

狄利克雷卷积:对两个函数$f,g$记他们的狄利克雷卷积为:$$ (f*g)(n)=\sum_{d|n}f(d)*g(\frac{n}{d}) $$ 很明显$(f*g)=(g*f)$

狄利克雷卷积的单位元函数为$e(n)=[n=1]$,即$f*e=e*f=f$,如果$f,g$都是积性函数那么$f*g$也是积性函数,函数集和狄利克雷卷积能够构成群

莫比乌斯反演:$g(n)=\sum_{d|n}f(d) \Rightarrow f(n)=\sum_{d|n}\mu(d)g(\frac{n}{d})$,反过来也成立,写成卷积的形式就是$g=1*f \Rightarrow f=\mu*g$,证明就是给$g=1*f$两边卷上$\mu$,变成$g*\mu=1*f*\mu=1*\mu*f=e*f=f$

杜教筛

问题:快速计算一个$f$的前缀和:$S(n)=\sum_{i=1}^n f(i)$(比如欧拉函数,$n$可以到1e10的级别)

假设找到了一个合适的函数$g$,他们卷积的前缀和(QwQ我这里写的会比较繁琐很多…):

$$ \begin{aligned} \sum_{i=1}^n(f*g)(n) & = \sum_{i=1}^n \sum_{d \mid i} g(d)f( \frac{i}{d})  = \sum_{i=1}^n \sum_{d=1}^n [d|i] g(d) f(\frac{i}{d}) \\ & = \sum_{d=1}^n g(d) \sum_{i=1}^n [d|i] f(i)   = \sum_{d=1}^n g(d) \sum_{i=1}^n \sum_{t=1}^n [i=dt] f(i) \\ & = \sum_{d=1}^n g(d) \sum_{t=1}^n \sum_{i=1}^n [i=dt] f(i)  = \sum_{d=1}^n g(d) \sum_{t=1}^{\lfloor\frac{n}{d} \rfloor } \sum_{i=1}^n [i=dt]f(i) \\ & = \sum_{d=1}^n g(d) \sum_{t=1}^{\lfloor \frac{n}{d} \rfloor } f(i) = \sum_{d=1}^n g(d)S(\lfloor \frac{n}{d} \rfloor)\end{aligned}$$

(最后三行是因为当$t>\lfloor \frac{n}{d} \rfloor$时后面的和式为0,所以这时候可以直接把上界缩小,以及对于每一对$d,t$页只有一个$i$与之对应,这时候后面的和式为1)

然后把$d=1$的提出来就得到了:

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

这里的$\lfloor \frac{n}{i} \rfloor$一共只有$O(\sqrt{n})$种取值,如果能够快速处理掉$g(1)$以及算出卷积的前缀和跟每一段连续的$S(\lfloor \frac{n}{i} \rfloor)$的和那么根据[1]里面说的就可以在$O(n^{\frac{3}{4}})$的时间里算出$S(n)$,如果$f$是一个积性函数,那么利用积性可以筛出前$n^{\frac{2}{3}}$项让复杂度降到$O(n^{\frac{2}{3}})$级别。

实现上有个trick就是可以用两个数组把答案存下来,对$<\sqrt{N}$的$x$放到$p1[x]$里,对$>=\sqrt{N}$的放到$p2[n/x]$里[3]

类似这样

inline lint & get(lint x)
{
    if(x<G)return p1[x];
    return p2[n/x];
}

如果有预处理的话这里也可以省去这个p1

举栗子

1.莫比乌斯函数的前缀和(51nod1244)

根据$\mu$的定义卷上一个$1$就得到了$[n=1]$,代到上面的式子直接得到$S(n)=1-\sum_{i=2}^nS(\lfloor \frac{n}{i} \rfloor)$

2.欧拉函数的前缀和(51nod1239)

注意到$\sum_{d|n} \varphi(d)=n$,这就相当于$\varphi * 1=n$,同样直接代进去就可以得到$S(n)=\frac{n(n+1)}{2}-\sum_{i=2}^n S(\lfloor \frac{n}{i} \rfloor)$

(话说bzoj3944那个就是要同时求$\mu$跟$\phi$的前缀和,然后我一开始还想着先全部读进来然后让最大的那个等于$n$然后再去配合前面那个trick算它可以快一点…然后发现这样根本就是错的…根本就不是一个n…)

(bzoj那题的代码)

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
  
typedef long long lint;
  
typedef pair<lint,lint> pa;
  
const lint N=2000005;
const lint G=100005;
  
lint T,tot,n;
lint pri[N],phi[N],mu[N];
bool p[N];
pa p1[G],p2[G];
  
inline pa & get(lint x)
{
    if(x<G)return p1[x];
    return p2[n/x];
}
  
inline void init()
{
    phi[1]=mu[1]=1;p[1]=1;
    for(register int i=2;i<N;i++)
    {
        if(!p[i])
        {
            pri[++tot]=i;phi[i]=i-1;mu[i]=-1;
        }
        for(register int j=1;j<=tot&&i*pri[j]<N;j++)
        {
            p[i*pri[j]]=1;
            if(i%pri[j]==0)
            {
                phi[i*pri[j]]=phi[i]*pri[j];mu[i*pri[j]]=0;break;
            }phi[i*pri[j]]=phi[i]*(pri[j]-1);mu[i*pri[j]]=-mu[i];
        }
    }
    for(register int i=1;i<N;i++)phi[i]+=phi[i-1],mu[i]+=mu[i-1];
}
  
inline pa calc(lint x)
{
    if(x<N)return (pa){phi[x],mu[x]};
    pa & cur=get(x);
    if((cur.first)!=-1)return cur;
    lint res1,res2,pos;res1=x*(x+1)/2;res2=1;
    for(register lint i=2;i<=x;i=pos+1)
    {
        pos=x/(x/i);pa temp=calc(x/i);
        res1-=(pos-i+1)*temp.first;
        res2-=(pos-i+1)*temp.second;
    }return cur=(pa){res1,res2};
}
  
int main()
{
    scanf("%lld",&T);init();memset(p1,-1,sizeof(p1));
    while(T--)
    {
        memset(p2,-1,sizeof(p2));
        scanf("%lld",&n);
        if(n<N)printf("%lld %lld\n",phi[n],mu[n]);
        else
        {
            pa temp=calc(n);
            printf("%lld %lld\n",temp.first,temp.second);
        }
    }
    return 0;
}
View Code

3.$f(n)=\varphi(n)*n$的前缀和$S(n)$

这时候可以构造出$g(n)=n$然后卷上$f$就得到:$(f*g)(n)=\sum_{d|n} \varphi(d)*d*\frac{n}{d}=n \sum_{d|n} \varphi(d)=n^2$,然后就得到了$S(n)=\frac{n(n+1)(2n+1)}{6}-\sum_{i=2}^n i*S(\lfloor \frac{n}{i} \rfloor)$,每一段相同的$S(\lfloor \frac{n}{i} \rfloor)$前面的系数也可以直接算出来

4.已知函数$f(n):\sum_{d|n}f(d)=n^2-3n+2$,求$\sum{i=1}^n f(i)$,多组数据,$n<=1e9$(HDU5608)

注意到已知的条件可以写成$1*f=n^2-3n+2$右边是一个很好求和的东西,那根据上面的结论就得到了$f$的前缀和$S(n)=\sum_{i=1}^n (i^2-3i+2)-\sum_{i=2}^nS(\lfloor \frac{n}{i}\rfloor)$,经过一番化简就得到了$S(n)=\frac{n(n-1)(n-2)}{3}-\sum_{i=2}^nS(\lfloor \frac{n}{i} \rfloor)$,对于前面小的数据用莫比乌斯反演预处理,注意数据范围~

(存个代码)

#include<cstdio>
#include<cstring>

typedef long long lint;

const lint MOD=(1e9+7);
const lint N=1000005;
const lint G=40005;

lint T,n,tot,inv3;
lint p2[G],mu[N],pri[N],sum[N];
bool p[N];

inline lint pow_mod(lint a,lint b,lint p)
{
    lint res=1;
    for(;b;b>>=1,a=(a*a)%p)if(b&1)res=(res*a)%p;
    return res;
}

inline void init()
{
    mu[1]=1;p[1]=1;
    for(register int i=2;i<N;i++)
    {
        if(!p[i])
        {
            pri[++tot]=i;mu[i]=-1;
        }
        for(register int j=1;j<=tot&&i*pri[j]<N;j++)
        {
            lint temp=i*pri[j];p[temp]=1;
            if(i%pri[j]==0){mu[temp]=0;break;}
            mu[temp]=-mu[i];
        }
    }
    for(register lint i=1;i<N;i++)
        for(register lint j=i;j<N;j+=i)sum[j]=(sum[j]+(i*i%MOD-3*i%MOD+2)*mu[j/i])%MOD;
    for(register int i=1;i<N;i++)sum[i]+=sum[i-1],sum[i]=(sum[i]%MOD+MOD)%MOD;
}

inline lint calc(lint x)
{
    if(x<N)return sum[x];
    lint &cur=p2[n/x];
    if(cur!=-1)return cur;
    lint res,pos;
    res=x*(x-1)%MOD*(x-2)%MOD*inv3%MOD;
    for(register lint i=2;i<=x;i=pos+1)
    {
        pos=x/(x/i);
        res-=(pos-i+1)*calc(x/i);
        while(res<0)res+=MOD;
    }res=(res%MOD+MOD)%MOD;
    return cur=res;
}

int main()
{
    inv3=pow_mod(3,MOD-2,MOD);
    scanf("%lld",&T);init();
    while(T--)
    {
        scanf("%lld",&n); 
        if(n<N)printf("%lld\n",sum[n]);
        else
        {
            memset(p2,-1,sizeof(p2));
            printf("%lld\n",calc(n));
        }
    }
    return 0;
}
View Code

待更...

原文地址:https://www.cnblogs.com/yoshinow2001/p/8017364.html