杜教筛

Posted hehe54321

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了杜教筛相关的知识,希望对你有一定的参考价值。

现在要算$s(n)=sum_{i=1}^n{f(i)}$

那么构造函数h和g,使得$h=f*g$,要求h和g的前缀和都好算

$sum_{i=1}^n{h(i)}=sum_{i=1}^n{sum_{d|i}f(frac{i}{d})g(d)}=sum_{i=1}^n{sum_{d=1}^i[d|i]f(frac{i}{d})g(d)}$
$=sum_{d=1}^ng(d)sum_{i=1}^{{lfloor}frac{n}{d}{ floor}}f(i)=sum_{d=1}^ng(d)s({lfloor}frac{n}{d}{ floor})$

$g(1)s(n)=sum_{i=1}^n{h(i)}-sum_{d=2}^ng(d)s({lfloor}frac{n}{d}{ floor})$

由于h和g前缀和都好算,可以对后面s(n/d)做整除分块,递归去算(要记忆化,用个哈希表之类的)

并不知道为什么要这么搞。。

复杂度?直接贴过来吧。。

以下转自https://www.cnblogs.com/onioncyc/p/8461267.html

根据(x/a)/b=x/(ab),杜教筛只用到2√n个值,所以杜教筛的复杂度是:

$sum_{i=1}^{sqrt n}O(sqrt i)+sum_{i=1}^{sqrt n}O(sqrt{frac{n}{i}})$

计算复杂度需要用到积分:

$int_{0}^{a}x^n=frac{1}{n+1}a^{n+1}$

所以,前半部分的复杂度:

$sum_{i=1}^{sqrt n}O(sqrt i)=int_{0}^{sqrt n}x^frac{1}{2}approx sqrt n^{frac{1}{2}+1}=O(n^{frac{3}{4}})$

后半部分的复杂度:(把√n提出来,最后是n^(3/4))

$sum_{i=1}^{sqrt n}O(frac{1}{sqrt i})=int_{0}^{sqrt n}x^{-frac{1}{2}}approx sqrt n^{-frac{1}{2}+1}=O(n^{frac{1}{4}})$

最终,复杂度$O(n^{frac{3}{4}})$

有的时候,可以用线性筛先预处理出一部分s函数的值,然后再算

此时的复杂度:转自https://blog.csdn.net/skywalkert/article/details/50500009

假设预处理了前k个正整数的$phi(n)$,且$kgesqrt{n}$,则复杂度变为$T(n)=sum_{i=1}^{frac{n}{k}}{sqrt{frac{n}{i}}}=O(frac{n}{sqrt{k}})$,当$k=O(n^frac{2}{3})$时可以取到较好的复杂度$T(n)=O(n^frac{2}{3})$

具体一点,就是根据前面那个一样去积分化简一下得到T(n)约等于$frac{n}{k^{frac{1}{2}}}+k$,就好算了


欧拉函数之和 51Nod - 1239

这个题就是要求欧拉函数的前缀和

显然,$1*varphi=id$

然后代进去,得到$s(n)=n(n+1)/2-sum_{d=2}^ns({lfloor}frac{n}{d}{ floor})$

 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cstring>
 4 #include<vector>
 5 #include<cmath>
 6 #include<unordered_map>
 7 using namespace std;
 8 #define fi first
 9 #define se second
10 #define mp make_pair
11 #define pb push_back
12 typedef long long ll;
13 typedef unsigned long long ull;
14 typedef pair<int,int> pii;
15 #define md 1000000007
16 unordered_map<ll,ll> ma;
17 ll n,K;
18 ll ma2[5001000],prime[1001000],len;
19 bool nprime[5001000];
20 ll calc(ll n)
21 {
22     if(n<=K)    return ma2[n];
23     if(ma.count(n))    return ma[n];
24     ll i,j,ans=0;
25     for(i=2;i<=n;i=j+1)
26     {
27         j=min(n,n/(n/i));
28         //printf("a%lld %lld
",i,j);
29         ans=(ans+(j-i+1)%md*calc(n/i)%md)%md;
30     }
31     //printf("b%lld %lld
",n,ans);
32     return ma[n]=(__int128(n)*(n+1)/2-ans+md)%md;
33 }
34 int main()
35 {
36     ll i,j;
37     scanf("%lld",&n);K=pow(n,0.666667);
38     ma2[1]=1;
39     for(i=2;i<=K;i++)
40     {
41         if(!nprime[i])    prime[++len]=i,ma2[i]=i-1;
42         for(j=1;j<=len&&i*prime[j]<=K;j++)
43         {
44             nprime[i*prime[j]]=1;
45             if(i%prime[j]==0)    {ma2[i*prime[j]]=ma2[i]*prime[j];break;}
46             else    ma2[i*prime[j]]=ma2[i]*(prime[j]-1);
47         }
48     }
49     for(i=1;i<=K;i++)    ma2[i]=(ma2[i]+ma2[i-1])%md;
50     //for(i=1;i<=K;i++)    printf("%lld %lld
",i,ma2[i]);
51     printf("%lld",calc(n));
52     return 0;
53 }

莫比乌斯函数之和 51Nod - 1244

显然,$mu*1=epsilon$

$s(n)=1-sum_{d=2}^ns({lfloor}frac{n}{d}{ floor})$

 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cstring>
 4 #include<vector>
 5 #include<cmath>
 6 #include<unordered_map>
 7 using namespace std;
 8 #define fi first
 9 #define se second
10 #define mp make_pair
11 #define pb push_back
12 typedef long long ll;
13 typedef unsigned long long ull;
14 typedef pair<int,int> pii;
15 unordered_map<ll,ll> ma;
16 ll n,K;
17 ll ma2[5001000],prime[1001000],len;
18 bool nprime[5001000];
19 ll calc(ll n)
20 {
21     if(n<=K)    return ma2[n];
22     if(ma.count(n))    return ma[n];
23     ll i,j,ans=0;
24     for(i=2;i<=n;i=j+1)
25     {
26         j=min(n,n/(n/i));
27         //printf("a%lld %lld
",i,j);
28         ans=ans+(j-i+1)*calc(n/i);
29     }
30     //printf("b%lld %lld
",n,ans);
31     return ma[n]=1-ans;
32 }
33 ll solve(ll n)
34 {
35     ll i,j;
36     //ma.clear();
37     K=pow(n,0.666667);
38     for(i=1;i<=K;i++)    nprime[i]=0;
39     len=0;
40     ma2[1]=1;
41     for(i=2;i<=K;i++)
42     {
43         if(!nprime[i])    prime[++len]=i,ma2[i]=-1;
44         for(j=1;j<=len&&i*prime[j]<=K;j++)
45         {
46             nprime[i*prime[j]]=1;
47             if(i%prime[j]==0)    {ma2[i*prime[j]]=0;break;}
48             else    ma2[i*prime[j]]=-ma2[i];
49         }
50     }
51     for(i=1;i<=K;i++)    ma2[i]=ma2[i]+ma2[i-1];
52     return calc(n);
53 }
54 int main()
55 {
56     ll a,b;
57     scanf("%lld%lld",&a,&b);
58     printf("%lld",solve(b)-solve(a-1));
59     return 0;
60 }

bzoj3944

洛谷P2413

看起来好像是双倍经验?然而过不去。。卡常。。。

可以用一个小技巧避免掉哈希表:注意到所有要存到表里面的s(k)都满足k=⌊n/i⌋且k>n^(2/3),那么i<n^(1/3),可以根据i来存

还有,筛法稍微多往前筛一点,毕竟多组数据

还有,把计算两个的过程和到一个函数里面,常数会小。。。

算是A掉了吧。。。(无视下面那个9.8s)

技术分享图片

 1 #pragma GCC optimize(3)
 2 #include<cstdio>
 3 #include<algorithm>
 4 #include<cstring>
 5 #include<vector>
 6 #include<cmath>
 7 #include<map>
 8 using namespace std;
 9 #define fi first
10 #define se second
11 #define mp make_pair
12 #define pb push_back
13 typedef long long ll;
14 typedef unsigned long long ull;
15 typedef pair<ll,ll> pll;
16 pll m[3010];
17 bool vm[3010];
18 ll n,K=5000000;
19 ll ma[5001000],ma2[5001000],prime[2001000],len;
20 bool nprime[5001000];
21 pll calc(ll n,ll t)
22 {
23     if(n<=K) return pll(ma[n],ma2[n]);
24     if(vm[t])  return m[t];
25     ll i,j,ans1=1,ans2=n*(n+1)/2;pll tt;
26     for(i=2;i<=n;i=j+1)
27     {
28         j=min(n,n/(n/i));tt=calc(n/i,t*i);
29         ans1-=(j-i+1)*tt.fi;
30         ans2-=(j-i+1)*tt.se;
31     }
32     vm[t]=1;
33     return m[t]=pll(ans1,ans2);
34 }
35 int main()
36 {
37     ll a,b,i,j,T;pll tt;
38     ma[1]=1;ma2[1]=1;
39     for(i=2;i<=K;i++)
40     {
41         if(!nprime[i])    prime[++len]=i,ma[i]=-1,ma2[i]=i-1;
42         for(j=1;j<=len&&i*prime[j]<=K;j++)
43         {
44             nprime[i*prime[j]]=1;
45             if(i%prime[j]==0)    {ma[i*prime[j]]=0;ma2[i*prime[j]]=ma2[i]*prime[j];break;}
46             else    ma[i*prime[j]]=-ma[i],ma2[i*prime[j]]=ma2[i]*(prime[j]-1);
47         }
48     }
49     for(i=1;i<=K;i++)    ma[i]=ma[i]+ma[i-1],ma2[i]+=ma2[i-1];
50     scanf("%lld",&T);
51     while(T--)
52     {
53         scanf("%lld",&a);
54         memset(vm,0,sizeof(vm));tt=calc(a,1);
55         printf("%lld %lld
",tt.se,tt.fi);
56     }
57     return 0;
58 }

 

以上是关于杜教筛的主要内容,如果未能解决你的问题,请参考以下文章

bzoj 3944 Sum —— 杜教筛

杜教筛题表(已完成)

杜教筛学习笔记

hdu5608杜教筛

杜教筛

杜教筛