zoukankan      html  css  js  c++  java
  • [日常摸鱼]积性函数求和——杜教筛

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

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

    以及在这里先感谢下小伙伴 @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

    待更...

  • 相关阅读:
    HDU 5059 Help him
    HDU 5058 So easy
    HDU 5056 Boring count
    HDU 5055 Bob and math problem
    HDU 5054 Alice and Bob
    HDU 5019 Revenge of GCD
    HDU 5018 Revenge of Fibonacci
    HDU 1556 Color the ball
    CodeForces 702D Road to Post Office
    CodeForces 702C Cellular Network
  • 原文地址:https://www.cnblogs.com/yoshinow2001/p/8017364.html
Copyright © 2011-2022 走看看