杜教筛
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$，就好算了
这个题就是要求欧拉函数的前缀和
显然，$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 }
显然，$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 }
看起来好像是双倍经验？然而过不去。。卡常。。。
可以用一个小技巧避免掉哈希表：注意到所有要存到表里面的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 }
以上是关于杜教筛的主要内容，如果未能解决你的问题，请参考以下文章