杜教筛

　　字符串好难啊,还是数论比较有趣.

　　从洛谷试炼场找了一个莫反的题做,非常套路,熟练的做完之后发现线性复杂度过不了,所以不得不学一下这种奇妙的筛法了.

　　设要求前缀和的函数$f$,是一个积性函数.

　　凭借灵感或者是经验找出另一个积性函数$g$来,两个函数卷一下.

　　$h=f imes g$

　　$sum_{i=1}^nh(i)=sum_{i=1}^nsum_{d|i}f(frac{i}{d})g(d)$

　　莫反常见套路:改变枚举顺序;

　　$sum_{i=1}^nh(i)=sum_{d=1}^ng(d)sum_{k=1}^{frac{n}{d}}f(k)$

　　$sum_{i=1}^nh(i)=sum_{d=1}^ng(d)S(frac{n}{d})$

　　$sum_{i=1}^nh(i)=g(1)S(n)+sum_{d=2}^ng(d)S(frac{n}{d})$

　　$S(n)g(1)=sum_{i=1}^nh(i)-sum_{d=2}^ng(d)S(frac{n}{d})$

　　那么$sum_{i=1}^nh(i)$怎么求呢?积性函数的卷积还是积性函数,显然可以线性筛!

　　不过要是这样我们为什么不直接筛$f$...

　　这就是考察技巧的时候了,需要选择一个巧妙的函数$g$使得$h$的前缀和非常好算,举几个例子:

　　$$sum_{i=1}^nvarepsilon(n)=[n>=1]$$

　　$$sum_{i=1}^n1=n$$

　　$$sum_{i=1}^ni=frac{n(n+1)}{2}$$

　　$$sum_{i=1}^ni^2=frac{n(n+1)(2n+1)}{6}$$

　　$$sum_{i=1}^n i^3=frac{n^2(n+1)^2}{4}$$

　　自然数的$k$次幂求和是一个关于$n$的$k+1$次多项式.可以用高斯消元暴力求每一项的系数,也可以用拉格朗日插值法做...好像说远了,回到杜教筛吧.

　　再说一点:今天烜神仙问我为什么一定就是$k+1$次多项式?目瞪口呆.jpg,因为之前学的课件上也只是写了"通过观察可以发现",但是经过不懈的yy,终于想到一个靠谱一点的证明:

　　　　首先明确一点:N维空间中的直线解析式是一个N次多项式,没有什么疑问吧.

　　　　对于k次方求和,我们将它表示到k+1维空间中,坐标为(n,n,...,n^k),构成了一个$n+1$维"正方体",那么它的体对角线长度就是

　　　　$$sqrt[n+1]{(n+1)x^{n+1}}=sqrt[n+1]{n+1}x$$

　　　是一个关于$x$的线性函数,而且每个立方体体对角线的角度都是相同的,沿着体对角线的方向这些点就构成了一条直线,它们的和就是在这条直线上一些等距离的点的求和,是一个k+1次多项式.

　　　　这么说可能不是很好理解,举两个简单的例子:

　　　　$i^1$: $i^2$:用电脑不大好画...

　　$$S(n)g(1)=sum_{i=1}^nh(i)-sum_{d=2}^ng(d)S(frac{n}{d})$$

　　前半部分$O(1)$求,后半部分除法分块,递归调用这个式子,一定要记忆化！

　　单纯考察杜教筛的题目不是很多，但是部分莫反的最后20-40分必须要低于线性才能做，每当看到$n<=10^{9}$或者$n<=10^{10}$的时候... ...

　　复杂度分析:不会积分.jpg，所以直接抄结论了,首先线性筛预处理出$n^{frac{2}{3}}$的答案,小于时直接返回,记忆化时不用$map$,而是用$f[x]$表示$f[n/x]$,可以做到$O(n^{frac{2}{3}}$的复杂度.出于这样的分析,我们可以发现需要记忆化的东西并不是很多,即使是多组询问也很快,复杂度几乎不变.实际运用中为了更快,可以在在不影响复杂度的情况下多线筛一些.注意...除法分块应从2开始分，至于为什么...看定义...

　　莫比乌斯函数求和：https://www.51nod.com/Challenge/Problem.html#!#problemId=1244

　　题意概述:求$sum_{i=a}^{b}mu(i),a,b<=10^{10}$

　　首先给莫比乌斯函数卷上一个别的函数,使得前缀和好求.

　　这个函数还是比较好找的,毕竟 $mu$ 的某个定义式就已经给了我们一些启发,

　　$$sum_{d|n}mu(d)=varepsilon(n)$$

　　你可能说这哪里卷积了? 我们可以假装它乘了一个一嘛.

　　$$sum_{d|n}mu(d)1(frac{n}{d})=varepsilon(n)$$

　　然后套模板就好了,不过这道题存在的意义就是为了讲解模板写法.

 1 # include <cstdio>
 2 # include <iostream>
 3 # include <cstring>
 4 # include <string>
 5 # define ll long long
 6 # define R register int
 7 
 8 using namespace std;
 9 
10 const int maxn=5000000;
11 int N=5000000,pri[maxn],h,vis[maxn+4],dvis[100005];
12 ll m[maxn+3],f[100005];
13 ll a,b,n;
14 
15 void init()
16 {
17     m[1]=1;
18     for (R i=2;i<=N;++i)
19     {
20         if(!vis[i]) pri[++h]=i,m[i]=-1;
21         for (R j=1;j<=h&&i*pri[j]<=N;++j)
22         {
23             vis[ i*pri[j] ]=1;
24             if(i%pri[j]==0) break;
25             m[ i*pri[j] ]=-m[i];
26         }
27     }
28     for (R i=2;i<=N;++i)
29         m[i]+=m[i-1];
30 }
31 
32 ll get_mu (ll x) { return (x<=N)?m[x]:f[n/x]; }
33 
34 void solve (ll x)
35 {
36     if(x<=N) return;
37     ll i=2,l,t=n/x;
38     if(dvis[t]) return;
39     dvis[t]=1;
40     f[t]=1;
41     while(i<=x)
42     {
43         l=x/(x/i);
44         solve(x/i);
45         f[t]-=(l-i+1)*get_mu(x/i);
46         i=l+1;
47     }
48 }
49 
50 int main()
51 {
52     scanf("%lld%lld",&a,&b);
53     if(a>b) swap(a,b);
54     a--;
55     if(a<0) a=0;
56     init();
57     ll ans=0;
58     if(b<=N)
59     {
60         ans+=m[b];
61         n=a;
62         if(a<=N) ans-=m[a]; else solve(a),ans-=f[1];
63     }
64     else
65     {
66         n=b;
67         solve(b);
68         ans+=f[1];
69         n=a;
70         if(a<=N) ans-=m[a];
71          else 
72         {
73             memset(dvis,0,sizeof(dvis));
74             solve(a);
75             ans-=f[1];
76         }
77     }
78     printf("%lld",ans);
79     return 0;
80 }

1244

　　欧拉函数求和:https://www.51nod.com/Challenge/Problem.html#!#problemId=1239

　　题意概述:求$sum_{i=1}^{n}varphi(i),n<=10^{10}$

　　考虑之前证明过的一个性质,可以快速的想出合理的卷积函数.

　　$$id=varphi imes 1$$

　　然后就是套板子啦.注意因为欧拉函数的和的范围会比较大,要开$ll$,有的时候要用慢速乘.

 1 # include <cstdio>
 2 # include <iostream>
 3 # include <cstring>
 4 # include <map>
 5 # define R register int
 6 # define ll long long
 7 
 8 using namespace std;
 9 
10 const int maxn=5000006;
11 const ll mod=1000000007;
12 int N=5000000,pri[maxn],h,phi[maxn];
13 ll n,vis[100005],f[100005];
14 map <ll,ll> m;
15 
16 void init()
17 {
18     phi[1]=1;
19     for (R i=2;i<=N;++i)
20     {
21         if(!phi[i]) pri[++h]=i,phi[i]=i-1;
22         for (R j=1;j<=h&&i*pri[j]<=N;++j)
23         {
24             if(i%pri[j]) phi[ i*pri[j] ]=1LL*phi[i]*(pri[j]-1)%mod;
25             else 
26             {
27                 phi[ i*pri[j] ]=1LL*phi[i]*pri[j]%mod;
28                 break;
29             }
30         }
31     }
32     for (R i=1;i<=N;++i)
33         phi[i]=(phi[i]+phi[i-1])%mod;
34 }
35 
36 ll get_phi (ll x) 
37 {
38     return (x<=N)?phi[x]:f[n/x]; 
39 }
40 
41 ll mul (ll a,ll b)
42 {
43     ll s=0;
44     while(b)
45     {
46         if(b&1LL) s=(s+a)%mod;
47         a=(a+a)%mod;
48         b>>=1;
49     }
50     return s;
51 }
52 
53 void solve (ll x)
54 {
55     if(x<=N) return;
56     ll i=2,l,t=n/x;
57     if(vis[t]) return;
58     f[t]=0;
59     vis[t]=1;
60     if(m[x]) return;
61     while(i<=x)
62     {
63         l=x/(x/i);
64         solve(x/i);
65         f[t]=(f[t]+(l-i+1)*get_phi(x/i)%mod)%mod;
66         i=l+1;
67     }
68     ll y=x+1;
69     if(y%2LL) x>>=1; else y>>=1;
70     f[t]=((mul(x,y)-f[t])%mod+mod)%mod;
71 }
72 
73 int main()
74 {
75     scanf("%lld",&n);
76     if(n<=N) N=n;
77     init();
78     if(n<=N) printf("%d",phi[n]);
79     else
80     {
81         solve(n);
82         printf("%lld",f[1]);
83     }
84     return 0;
85 }

1239

　　Sum:https://www.lydsy.com/JudgeOnline/problem.php?id=3944

　　题意概述:上述两个题的二合一.

 1 # include <cstdio>
 2 # include <iostream>
 3 # include <cstring>
 4 # define R register int
 5 # define ll long long
 6 
 7 using namespace std;
 8 
 9 const int maxn=4000000;
10 int t,N=4000000,pri[maxn+2],h,vis[100005];
11 ll n,mu[maxn+2],phi[maxn+2],m[100005],p[100005];
12 
13 void init()
14 {
15     mu[1]=1,phi[1]=1;
16     for (R i=2;i<=N;++i)
17     {
18         if(!phi[i]) phi[i]=i-1,pri[++h]=i,mu[i]=-1;
19         for (R j=1;j<=h&&i*pri[j]<=N;++j)
20         {
21             if(i%pri[j]==0)
22             {
23                 phi[ i*pri[j] ]=phi[i]*pri[j];
24                 break;
25             }
26             phi[ i*pri[j] ]=phi[i]*phi[ pri[j] ];
27             mu[ i*pri[j] ]=-mu[i];
28         }
29     }
30     for (R i=1;i<=N;++i)
31         mu[i]+=mu[i-1],phi[i]+=phi[i-1];
32 }
33 
34 ll smu (ll x) { return x<=N?mu[x]:m[n/x]; }
35 ll sphi (ll x) { return x<=N?phi[x]:p[n/x]; }
36 
37 void solve (ll x)
38 {
39     if(x<=N) return;
40     ll i=2,l,t=n/x;
41     if(vis[t]) return;
42     vis[t]=true;
43     p[t]=(x*(x+1))/2LL,m[t]=1;
44     while(i<=x)
45     {
46         l=x/(x/i);
47         solve(x/i);
48         p[t]-=sphi(x/i)*(l-i+1);
49         m[t]-=smu(x/i)*(l-i+1);
50         i=l+1;
51     }
52 }
53 
54 int main()
55 {
56     scanf("%d",&t);
57     init();
58     while(t--)
59     {
60         scanf("%lld",&n);
61         if(n<=N)
62             printf("%lld %lld
",phi[n],mu[n]);
63         else
64         {
65             memset(vis,0,sizeof(vis));
66             solve(n);
67             printf("%lld %lld
",p[1],m[1]);
68         }
69     }
70     return 0;
71 }

3944

　　许多莫反的题都可以用杜教筛来优化,举个例子:

　　GCDSUM：https://www.51nod.com/Challenge/Problem.html#!#problemId=1237

　　题意概述:求$sum_{i=1}^nsum_{j=1}^n(i,j),n<=10^{10}$.

　　套路题，枚举答案，前面除法分块，后面杜教筛求和。注意，因为这道题其实是变相的多组询问，所以不要用之前的技巧(要多次清空数组)，反而是unordered_map快一些.

　　$$sum_{d=1}^{n}dsum_{i=1}^{frac{n}{d}} varphi(i)$$

　　$$sum_{d=1}^nS(frac{n}{d})d$$

　　这样化出来的式子会少算一半,答案乘二后,如果$i,j$相同又会多算一次,减掉就好了.

 1 # include <cstdio>
 2 # include <iostream>
 3 # include<tr1/unordered_map>
 4 # define ll long long
 5 # define R register int
 6 
 7 using namespace std;
 8 
 9 const ll mod=1000000007;
10 const int maxn=5000006;
11 int N=5000000,h,cnt,pri[maxn],phi[maxn];
12 ll n;
13 tr1::unordered_map <ll,ll> m;
14 
15 void init()
16 {
17     phi[1]=1;
18     for (R i=2;i<=N;++i)
19     {
20         if(!phi[i]) pri[++h]=i,phi[i]=i-1;
21         for (R j=1;j<=h&&i*pri[j]<=N;++j)
22         {
23             if(i%pri[j])
24                 phi[ i*pri[j] ]=1LL*phi[i]*(pri[j]-1)%mod;
25             else
26             {
27                 phi[ i*pri[j] ]=1LL*phi[i]*pri[j]%mod;
28                 break;
29             }
30         }
31     }
32     for (R i=1;i<=N;++i)
33         phi[i]=(phi[i]+phi[i-1])%mod;
34 }
35 
36 ll mul (ll a,ll b)
37 {
38     ll s=0;
39     while(b)
40     {
41         if(b&1LL) s=(s+a)%mod;
42         a=(a+a)%mod;
43         b>>=1LL;
44     }
45     return s;
46 }
47 
48 ll get_phi (ll x) {    return (x<=N)?phi[x]:m[x]; }
49 
50 ll S (ll x)
51 {
52     ll y=x+1;
53     if(y%2) x=x/2; else y=y/2;
54     return mul(x,y);
55 }
56 
57 ll solve (ll x)
58 {
59     if(x<=N) return phi[x];
60     ll i=2,l,ans=0;
61     if(m[x]) return m[x];
62     while(i<=x)
63     {
64         l=x/(x/i);
65         solve(x/i);
66         ans=(ans+(l-i+1)*get_phi(x/i));
67         i=l+1;
68     }
69     ans=((S(x)-ans)%mod+mod)%mod;
70     return m[x]=ans;
71 }
72 
73 int main()
74 {
75     scanf("%lld",&n);
76     init();
77     ll i=1,l,ans=0,a;
78      while(i<=n)
79     {
80         l=n/(n/i);
81         a=solve(n/i);
82         ans=(ans+mul((S(l)-S(i-1)+mod)%mod,a))%mod;
83         i=l+1;
84     }
85     ans=ans*2%mod;
86     ans=(ans-S(n)+mod)%mod;
87     printf("%lld",ans);
88     return 0;
89 }

1237

　---shzr