杜教筛
逃不过ww。被模拟赛人均爆切的莫比乌斯反演T1的1e8数据范围逼着爬来学杜教筛了。
杜教筛是用来筛积性函数前缀和的筛法。复杂度 (O(n^{frac 2 3}))
前置芝士
@莫比乌斯反演 笨蛋莫反还没搞清楚就被逼迫来学。
(h(i)) 是积性函数,计算 (sum_{i=1}^n f(i)) 。
套路式:
构造两个积性函数 (h,g) 使得 (h=f*g)
令 (s(n)=sum_{i=1}^n f(i))
只要 (h(i)) 的前缀和好求,那么对之后的柿子整除分块,求s(n)的时间复杂度是(O(n^{frac 2 3}))
构造 (g,h) 只能靠经验/kk 太草了
举例
-
求 (s(n)=sum_{i=1}^n mu(i))
[g(1)·s(n)=sum_{i=1}^n h(i)-sum_{d=2}^n g(d)s(lfloor frac n d floor) ]寻找一个积性函数 (g) 使得 $f·mu =h $ 且 (h) 前缀和非常好求。
(mu*I=epsilon) 显然非常好求。
[s(n)=1-sum_{d=2}^n s(lfloor frac n d floor)\ ]就这样……
-
求 (s(n)=sum_{i=1}^n phi(i))
[g(1)·s(n)=sum_{i=1}^n h(i) -sum_{d=2}^n g(d)s(lfloor frac n d floor) ]我们有:
[varphi * I =id ](g:I,h:id)
[s(n)=frac {n*(n+1)} 2-sum_{d=2}^n s(lfloor frac n d floor) ] -
求 (S(n)=sum_{i=1}^{n}icdot varphi(i))
[g(1)·s(n)=sum_{i=1}^n h(i) -sum_{d=2}^n g(d)s(lfloor frac n d floor) ][ f*g=h\ h(x)=sum_{d|x} f(x)g(frac x d)=sum_{d|x}( d·mu(d))·g(frac x d)\ 为了把d 搞掉,尝试让g为id\ =sum_{d|x} mu(d)·x\ =xsum_{d|x} mu(d)\ =x^2 ][g(1)·s(n)=sum_{i=1}^n h(i) -sum_{d=2}^n g(d)s(lfloor frac n d floor)\ s(n)=frac {n·(n+1)·(2n+1)} 6-sum_{d=2}^n d·s(lfloor frac n d floor) ]
三个题最后都是整除分块一下求就行了。
代码实现
map存一下杜教筛的前缀和,然后递归实现的形式很有意义(?),就这样。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int N=5e6+10;
bool vis[N];
int cnt,mu[N],phi[N],p[N];
ll summu[N],sumphi[N];
map<ll,ll>mpmu;
map<ll,ull>mpphi;
int bmin(ll x,ll y){ return (x>y)?y:x; }
void init(ll n){
mu[1]=1; vis[1]=1; phi[1]=1;
for(ll i=2;i<=n;i++){
if(!vis[i]) p[++cnt]=i,mu[i]=-1,phi[i]=i-1;
for(ll j=1;j<=cnt&&p[j]*i<=n;j++){
vis[p[j]*i]=1;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
mu[i*p[j]]=-mu[i];
phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
for(ll i=1;i<=n;i++){
summu[i]=summu[i-1]+mu[i];
sumphi[i]=sumphi[i-1]+phi[i];
}
return;
}
ll getsummu(ll n){
if(n<=N-10) return summu[n];
else if(mpmu.count(n)) return mpmu[n];
//杜教筛
ll ret=1;
for(ll l=2,r;l<=n;l=r+1){
r=bmin(n,n/(n/l));
ret-=getsummu(n/l)*(r-l+1);
}
return mpmu[n]=ret;
}
ull getsumphi(ll n){
if(n<=N-10) return sumphi[n];
else if(mpphi.count(n)) return mpphi[n];
//杜教筛
ull ret=1ll*n*(n+1)/2;
for(ll l=2,r;l<=n;l=r+1){
r=bmin(n,n/(n/l));
ret-=(ull)getsumphi(n/l)*(r-l+1);
}
return mpphi[n]=ret;
}
int main(){
int t; scanf("%d",&t);
while(t--){
ll n; scanf("%lld",&n);
init(N-10);
printf("%llu %lld
",getsumphi(n),getsummu(n));
}
return 0;
}
/*
sum_{d=1}^{n} mu(d) frac {(t-1)·(t-2)} 4
mu: s(n)=1-sum_{d=2}^n s(lfloor frac n d
floor)
*/
例题
LGP3768 简单的数学题
给定 (n) ,求
(1le nle 10^{10})
不知道为甚么代码过不去样例,气炸了。
为了抄题解代码爬去推式子。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll N=5e6+10;
int p,cnt,sum[N],phi[N],prim[N],inv2,inv6;
bool vis[N];
map<ll,int>mp;
int power(int a,int b){
int ret=1;
while(b){
if(b&1) ret=1ll*a*ret%p;
a=1ll*a*a%p; b>>=1;
}
return ret;
}
void init(ll n){
phi[1]=1; vis[1]=1;
for(ll i=2;i<=n;i++){
if(!vis[i]) vis[i]=1,phi[i]=i-1,prim[++cnt]=i;
for(int j=1;j<=cnt&&prim[j]*i<=n;j++){
vis[prim[j]*i]=1;
if(i%prim[j]==0){
phi[i*prim[j]]=1ll*phi[i]*prim[j]%p;
break;
}
phi[i*prim[j]]=1ll*phi[i]*(prim[j]-1)%p;
}
}
for(ll i=1;i<=n;i++) sum[i]=(sum[i-1]+1ll*i*i%p*phi[i]%p)%p;
return;
}
int calc(ll x){ return x%p*(x%p+1)%p*inv2%p; }
int pfh(ll x){ return x%p*(x%p+1)%p*(2*x%p+1)%p*inv6%p; }
int getphi(ll n){
if(n<=N-10) return sum[n];
if(mp.count(n)) return mp[n];
int ret=1ll*calc(n)*calc(n)%p;
for(ll l=2,r;l<=n;l=r+1){
r=min(n,n/(n/l));
ret=(ret-1ll*(pfh(r)-pfh(l-1))*getphi(n/l)%p)%p;
}
return mp[n]=ret;
}
int main(){
ll n;
scanf("%d%lld",&p,&n);
inv6=power(6,p-2);inv2=power(2,p-2);
init(N-10);
int ans=0;
for(ll l=1,r;l<=n;l=r+1){
r=min(n,n/(n/l));
int qwq=calc(n/l);
ans=(ans+1ll*qwq*qwq%p*(getphi(r)-getphi(l-1))%p)%p;
}
printf("%d
",(ans+p)%p);
return 0;
}
/*
998244353 2000
334905957
883968974
*/
LGP1587 [NOI2016]循环之美
给定n,m,k.
首先考虑在什么情况下可以化为纯循环小数。
当y是999...99的因数的时候
假设有 (m) 个 (k-1) 时,满足要求。(y| (k^m-1))
那么计算到多少个 (k) 的时候是可以判断不可能的呢。。。
我们观察十进制的,一个不合法的会变成(999...9000...0) 显然如果 (2|y) 或者 (5|y) (y) 就要变成这个了。
思考一下,只要 (y) 与 (k) 互质即可。
柿子变为
(G:)
(F:)
……原来自己推出来了一个,然后以为搞不出来就删掉了。然后其实是对的。我是zz。但是我懒得再搞。
(frac n d) 去掉先。
将原来的F与最后的F放在一起。
可以发现,后半部分是 (F(frac n t,t)) 这是可以递归的,杜教筛F。边界是 (F(1))
时间复杂度
(G: O(1)),(F: O(n^{frac 2 3}))
预处理 (O(klog k))
整除分块+杜教筛mu!
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e6+10,K=2e3+10;
bool vis[N];
int cnt,p[N],mu[N],sum[N],pre[K];
unordered_map<int,int>mp,f[K];
int gcd(int a,int b){ return (b==0)?a:gcd(b,a%b); }
void init(int n,int k){
for(int i=1;i<=k;i++) pre[i]=pre[i-1]+(gcd(i,k)==1);
mu[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i]) p[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&p[j]*i<=n;j++){
vis[p[j]*i]=1;
if(i%p[j]==0) break;
mu[p[j]*i]=-mu[i];
}
}
for(int i=1;i<=n;i++) sum[i]=sum[i-1]+mu[i];
return;
}
int getmu(int n){
if(n<=N-10) return sum[n];
if(mp.count(n)) return mp[n];
ll ret=1;
for(int l=2,r;l<=n;l=r+1){
r=n/(n/l);
ret-=(r-l+1)*getmu(n/l);
}
return mp[n]=ret;
}
int getf(int n,int k){
if(f[k].count(n)) return f[k][n];
if(n==0) return 0;
if(k==1) return f[k][n]=getmu(n);
ll ret=0;
for(int i=1;i*i<=k;i++){
if(k%i!=0) continue;
if(mu[i]) ret+=mu[i]*mu[i]*getf(n/i,i);
if(mu[k/i]&&i*i!=k) ret+=mu[k/i]*mu[k/i]*getf(n/(k/i),k/i);
}
return f[k][n]=ret;
}
int g(int n,int k){
return 1ll*pre[k]*(n/k)+pre[n%k];
}
int main(){
int n,m,k;
scanf("%d%d%d",&n,&m,&k);
init(N-10,k);
ll ans=0; int tmp=min(n,m);
for(int l=1,r;l<=tmp;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=1ll*(getf(r,k)-getf(l-1,k))*g(m/l,k)*(n/l);
}
printf("%lld
",ans);
return 0;
}