笔记：杜教筛

杜教筛

逃不过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))

[sum_{i=1}^n h(i)\ =sum_{i=1}^n sum_{d|i} g(d)*f(frac i d)\ =sum_{xd=1}^n sum_{d|i} g(d)*f(x)\ =sum_{x=1}^{frac n d} sum_{d|i} g(d)*f(x)（顺序不太对劲）\ = sum_{d|i} g(d) sum_{x=1}^{frac n d}f(x)\ =sum_{d=1}^n g(d) sum_{i=1}^{frac n d} f(i)\ 提取 sum_{i=1}^{frac n d} f(i) ，\ =sum_{d=1}^n g(d) s(frac n d)\ 提出左边的第一项。\ =g(1)·s(n)+sum_{d=2}^n g(d) s(frac n d)\ ∴ g(1)·s(n)=sum_{i=1}^n h(i)-sum_{d=2}^n g(d) s(lfloor frac n d floor) ]

只要 (h(i)) 的前缀和好求，那么对之后的柿子整除分块，求s(n)的时间复杂度是(O(n^{frac 2 3}))

构造 (g,h) 只能靠经验/kk 太草了

举例

1. 求 (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)\ ]

   就这样……

2. 求 (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) ]

3. 求 (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存一下杜教筛的前缀和，然后递归实现的形式很有意义（？），就这样。

LGP4213 【模板】杜教筛（Sum）

#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) ，求

[({sum_{i=1}^n sum_{j=1}^n {i·j·gcd(i,j)}}) mod p ]

(1le nle 10^{10})

[{sum_{i=1}^n sum_{j=1}^n {i·j·gcd(i,j)}}\ sum_{d=1}^n sum_{d|i} sum_{d|j} [gcd(i,j)=d] i·j·d\ sum_{d=1}^n sum_{i=1}^{frac n d} sum_{j=1}^{frac n d} [gcd(i,j)=1] i·j·d^3\ sum_{d=1}^n sum_{i=1}^{frac n d} sum_{j=1}^{frac n d} ijd^3sum_{t|gcd(i,j)} mu(t) \ sum_{d=1}^nd^3 sum_{t=1}^{frac n d} mu(t) sum_{i=1}^{frac n {dt}} sum_{j=1}^{frac n {dt}} ijt^2 \ sum_{d=1}^nd^3 sum_{t=1}^{frac n d} mu(t) t^2calc(frac {frac n d} {t})^2\ 先整除分块frac n d,然后整除分块t，mu 那里用一下杜教筛。\ 不对，你的复杂度爆炸了。 ]

[sum_{d=1}^nd^3 sum_{t=1}^{frac n d} mu(t) t^2calc(frac {n} {td})^2\ 令T=td,\ sum_{d=1}^nd^3 sum_{frac T d=1}^{frac n d} mu(frac T d) (frac T d)^2calc(frac {n} {T})^2\ sum_{d=1}^nd^3 sum_{T=d}^{n} mu(frac T d) (frac T d)^2calc(frac {n} {T})^2\ 把T丢出去，\ sum_{T=1}^n sum_{d|T} d^3 mu(frac T d) (frac T d)^2calc(frac {n} {T})^2\ 好了，该莫比乌斯反演了。\ f(T)= sum_{d|T} d^3 mu(frac T d) (frac T d)^2calc(frac {n} {T})^2\ 令 F(n)=sum_{n|T} f(T)\ 我谔谔，所以F并不好求啊搞他干什么。sum_{d=1}^nd^3 sum_{t=1}^{frac n d} mu(t) t^2calc(frac {n} {td})^2\ 令T=td,\ sum_{d=1}^nd^3 sum_{frac T d=1}^{frac n d} mu(frac T d) (frac T d)^2calc(frac {n} {T})^2\ sum_{d=1}^nd^3 sum_{T=d}^{n} mu(frac T d) (frac T d)^2calc(frac {n} {T})^2\ 把T丢出去，\ sum_{T=1}^n T^2calc(frac {n} {T})^2sum_{d|T} mu(frac T d) d\ 令f(x)=mu(x),g(x)=x\ sum_{d|T} mu(frac T d) d=f*g=mu*id=varphi\ \ 所以\ sum_{T=1}^ncalc(frac {n} {T})^2 T^2 varphi(T)\ sum_{T=1}^n frac {n^2*(frac n T +1)^2} 4 varphi(T)\ 然后，整除分块你的 frac n T 就行了。phi 用杜教筛。 ]

不知道为甚么代码过不去样例，气炸了。

为了抄题解代码爬去推式子。

[sum_{T=1}^ncalc(frac {n} {T})^2 T^2 varphi(T)\ f(T)=T^2varphi (T)\ f*g=h\ g(1)S(n)=sum_{i=1}^n h(i)-sum_{d=2}^n g(d)s(lfloor frac n d floor)\ h(x)=sum_{d|x} g(d) f(frac x d)\ h(x)=sum_{d|x} g(d) phi(frac x d) frac {x^2} {d^2}\ 设g(d)=d^2\ h(x)=sum_{d|x} phi(frac x d) {x^2}\ h(x)=x^2sum_{d|x} phi(frac x d)\ phi*I=id 所以h(x)=x^3？\ g(1)S(n)=sum_{i=1}^n h(i)-sum_{d=2}^n g(d)s(lfloor frac n d floor)\ S(n)=sum_{i=1}^n i^3 -sum_{d=2}^n d^2 s(frac n d)\ 整除分块frac n d\ sum_{i=1}^n i^3=frac {n(n+1)^2} 4 ]

#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.

[sum_{i=1}^n sum_{j=1}^m [gcd(i,j)=1&&frac i j 在k进制下 可以化为纯循环小数] ]

首先考虑在什么情况下可以化为纯循环小数。

当y是999...99的因数的时候

假设有 (m) 个 (k-1) 时，满足要求。(y| (k^m-1))

那么计算到多少个 (k) 的时候是可以判断不可能的呢。。。

我们观察十进制的，一个不合法的会变成(999...9000...0) 显然如果 (2|y) 或者 (5|y) (y) 就要变成这个了。

思考一下，只要 (y) 与 (k) 互质即可。

柿子变为

[sum_{i=1}^n sum_{j=1}^m [gcd(i,j)=1&&gcd(j,k)=1]\ 人傻了，换一下枚举顺序。\ sum_{j=1}^m [gcd(j,k)=1]sum_{i=1}^n [gcd(i,j)=1]\ 先考虑后面那坨？\ sum_{j=1}^m [gcd(j,k)=1]sum_{i=1}^n [gcd(i,j)=1]\ sum_{j=1}^m [gcd(j,k)=1]sum_{i=1}^n sum_{d|gcd(j,i)} mu(d)\ sum_{d=1}^{min(n,m)} mu(d) sum_{d|j}^m[gcd(j,k)=1] sum_{d|i}^n 1\ sum_{d=1}^{min(n,m)} mu(d) sum_{d|j}^m [gcd(j,k)=1] frac n d\ sum_{d=1}^{min(n,m)} frac n dmu(d) sum_{d|j}^m [gcd(j,k)=1] \ 关于k：他是一个可爱的2e3，不用担心。\ 所以你刚才tm以为n,m,k同阶拆出一个不能莫反的玩意是zz行为。\ sum_{d=1}^{min(n,m)} frac n dmu(d) sum_{j=1}^{frac m d} [gcd(jd,k)=1] \ 为了让它们不相关，\ sum_{d=1}^{min(n,m)} frac n dmu(d) [gcd(d,k)=1]sum_{j=1}^{frac m d} [gcd(j,k)=1] \ 令 F(n)=sum_{d=1}^{min(n,m)}[gcd(d,k)=1]frac n dmu(d) ,\ 令 G(frac m d)=sum_{j=1}^{frac m d} [gcd(j,k)=1]\ ]

(G:)

[G(m)=sum_{j=1}^{m} [gcd(j,k)=1] \ =lfloor frac m k floor varphi(k)+pre[m\%k]\ pre[m\%k]=sum_{i=1}^{m\%k} [gcd(i,k)=1]\ 进行O(klog k) 预处理即可 O(1) 求出答案。 ]

(F：)

……原来自己推出来了一个，然后以为搞不出来就删掉了。然后其实是对的。我是zz。但是我懒得再搞。

(frac n d) 去掉先。

[F(n)=sum_{d=1}^{min(n,m)}mu(d) [gcd(d,k)=1] \ =sum_{d=1}^{min(n,m)}mu(d)sum_{t|gcd(d,k)}mu(t) \ =sum_{d=1}^{min(n,m)}mu(d)sum_{t|k}mu(t) [t|d]\ =sum_{t|k}mu(t)sum_{d=1}^{frac {min(n,m)} t}mu(td)\ =sum_{t|k}mu(t)sum_{d=1}^{frac {min(n,m)} t}mu(t)mu(d)[gcd(t,d)=1]\ =sum_{t|k}mu(t)^2sum_{d=1}^{frac {min(n,m)} t}mu(d)[gcd(t,d)=1]\ ]

将原来的F与最后的F放在一起。

[F(n,k)=sum_{d=1}^{min(n,m)}mu(d) [gcd(d,k)=1] \ F(n,k)=sum_{t|k}mu(t)^2sum_{d=1}^{frac {min(n,m)} t}mu(d)[gcd(d,t)=1]\ ]

可以发现，后半部分是 (F(frac n t,t)) 这是可以递归的，杜教筛F。边界是 (F(1))

时间复杂度

(G: O(1))，(F: O(n^{frac 2 3}))

预处理 (O(klog k))

[ans=sum_{d=1}^n mu(d)[gcd(d,k)=1] G(frac m d) lfloor frac n d floor ]

整除分块+杜教筛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;
}