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

    待更...

  • 相关阅读:
    Insert into a Binary Search Tree
    Search in a Binary Search Tree
    Binary Search Tree Iterator
    Validate Binary Search Tree
    Serialize and Deserialize Binary Tree
    图的搜索
    codeforce vk cup2017
    hdu1160dp
    完全背包hdu1114
    最长递增子序列hdu1087
  • 原文地址:https://www.cnblogs.com/yoshinow2001/p/8017364.html
Copyright © 2011-2022 走看看