Min_25 筛与杜教筛

杜教筛

\(𝐹\)\(𝑓\) 的前缀和,\(𝐺\), \(𝐻\) 同理。

假设 \(𝑓 × 𝑔 = ℎ\) ,并且 \(𝐹, 𝐻\) 易求出,\(𝐺\) 难求出。

那么

\[H (𝑛) = \sum_{𝑖 \cdot 𝑗≤𝑛} 𝑓(𝑖) 𝑔(𝑗) = \sum_{𝑖≤𝑛} 𝑓(𝑖) 𝐺 (\frac {𝑛} {𝑖})\\ = f(1)\cdot 𝐺(𝑛) + \sum_{2≤𝑖≤𝑛} 𝑓(𝑖) 𝐺(\frac {𝑛} {𝑖})\]

有:

\[f(1)\cdot G(n)=H(n)-\sum_{2≤𝑖≤𝑛} 𝑓(𝑖) 𝐺(\frac {𝑛} {𝑖}) \]

整除分块,可以在 \(O(n^{\frac {3}{4}})\) 的时间内计算出 \(𝐺(n)\)

代码按照理解大概可以写成这样:

ll GetSum(int n){
	ll ans = H(n);
	for(ll l = 2, r; l <= n; l = r + 1) {
    	r = (n / (n / l)); 
 		ans -= (F(r) - F(l - 1)) * GetSum(n / l);
 	} return ans;
}

时间复杂度证明:

考虑到求 \(S(n)\) 时,需要求出 \(2 \times \sqrt n\)\(S(\lfloor \frac {n}{i} \rfloor)\)

由于 \(\lfloor \frac{\lfloor \frac {n}{a}\rfloor}{b} \rfloor =\lfloor \frac {n}{a b} \rfloor\) ,显然所需要求的 \(S(\lfloor \frac {n}{i} \rfloor)\) 数量只有 \(2 \times \sqrt n\) 个。

考虑到整除分块的复杂度,可得总复杂度为:

\[O(\sum_{i=1}^{\sqrt n} \sqrt i+ \sum_{i=1}^{\sqrt n} \sqrt \frac{n}{i})=O(n^{\frac {3}{4}}) \]

还可以进一步优化杜教筛,即先线性筛出前 \(m\) 个答案,之后再用杜教筛。

这个优化之后的复杂度是:

\[O(\sum_{i=1}^{\lfloor \frac{n}{m} \rfloor}\sqrt \frac{n}{i}) = O(\frac {n}{\sqrt m}) \]

\(m=n^{\frac {2}{3}}\) 时,总复杂度为 \(O(n^{\frac {2}{3}})\)

此时预处理复杂度与筛法复杂度相等,为杜教筛的最优复杂度。


\[\]

杜教筛小应用

比如说我们要求 \(g(𝑖)=𝜑(𝑖)\times 𝑖\) 的前缀和

那么观察 \(𝑔 * 𝑖𝑑 = 𝑖𝑑^2\) ,就令 \(𝑓 = 𝑖𝑑\)\(ℎ = 𝑖𝑑^2\) 即可。

比如说我们要求 \(g(𝑖)=𝜑 (𝑖) \times 𝑖^𝑘\) 的前缀和

那么观察 \(𝑔 * 𝑖𝑑^𝑘 = 𝑖𝑑^{𝑘+1}\) ,就令 \(𝑓 = 𝑖𝑑^𝑘\)\(ℎ = 𝑖𝑑^{𝑘+1}\) 即可。

比如说我们要求 \(g (𝑖) = 𝜇 (𝑖) \times 𝑖^𝑘\) 的前缀和

那么观察 \(𝑔 * 𝑖𝑑^𝑘 = 𝑖𝑑^𝑘∗ 𝜀 = 𝜀\),就令 \(𝑓 = 𝑖𝑑^𝑘\),ℎ = 𝜀 即可。


\[\]

例题:

【模板】杜教筛(Sum)

\(\sum_{i=0}^{n} \mu(i)\)

发现 \(\mu * I = e\)

有:

\[I(1) \cdot \sum_{i=1}^{n}\mu(i)=\sum_{i=1}^{n}e(i)-\sum_{i=2}^{n}I(i)\sum_{j=1}^{\lfloor \frac{n}{i} \rfloor}\mu(t) \]

\[\sum_{i=1}^{n}\mu(i)=1-\sum_{i=2}^{n}I(i)\sum_{j=1}^{\lfloor \frac{n}{i} \rfloor}\mu(t) \]

\[ans(n)=1-\sum_{i=2}^{n}I(i)\cdot ans(\lfloor \frac{n}{i} \rfloor) \]

\(\sum_{i=0}^{n} \varphi (i)\)

发现 \(\varphi * I = id\)

有:

\[I(1) \cdot \sum_{i=1}^{n}\varphi (i)=\sum_{i=1}^{n}id(i)-\sum_{i=2}^{n}I(i)\sum_{j=1}^{\lfloor \frac{n}{i} \rfloor}\varphi (t) \]

\[\sum_{i=1}^{n}\varphi (i)=\sum_{i=1}^{n}i-\sum_{i=2}^{n}I(i)\sum_{j=1}^{\lfloor \frac{n}{i} \rfloor}\varphi (t) \]

\[ans(n)=\frac{n * (n+1)}{2}-\sum_{i=2}^{n}I(i)ans(\lfloor \frac{n}{i} \rfloor) \]

点击查看代码

P3768 简单的数学题

求:

\[\sum_{i=1}^{n}\sum_{j=1}^{n} i\cdot j \cdot gcd(i,j) \]

枚举 \(gcd\)

\[\sum_{d=1}^{n}d \cdot \sum_{i=1}^{n}\sum_{j=1}^{n} i\cdot j \cdot [gcd(i,j)=d] \]

\[\sum_{d=1}^{n}d^3 \cdot \sum_{i=1}^{n/d}\sum_{j=1}^{n/d} i\cdot j \cdot [gcd(i,j)=1] \]

\[\sum_{d=1}^{n}d^3 \cdot \sum_{i=1}^{n/d}\sum_{j=1}^{n/d} i\cdot j \cdot \sum_{k|gcd(i,j)}\mu(k) \]

\[\sum_{d=1}^{n}d^3 \cdot \sum_{k=1}^{n/d}\mu(k)\cdot \sum_{k|i}^{n/d}\sum_{k|j}^{n/d} i\cdot j \]

\[\sum_{d=1}^{n}d^3 \cdot \sum_{k=1}^{n/d}\mu(k)*k^2\cdot \sum_{i=1}^{n/id}\sum_{j=1}^{n/id} i\cdot j \]

\[\sum_{d=1}^{n}d^3 \cdot \sum_{k=1}^{n/d}\mu(k)*k^2\cdot (\sum_{i=1}^{n/id} i)^2 \]


\[\]

[NOI2016] 循环之美

求:

\[\sum_{i=1}^{n}\sum_{j=1}^{m}[i\perp j][j\perp k] \]

\([i\perp j]\) ,得:

\[\sum_{i=1}^{n}\sum_{j=1}^{m}[i\perp j][j\perp k]\\ =\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{t|i,t|j}\mu(t)[j\perp k]\\ =\sum_{t|j}\mu(t)\sum_{t|i}^{n}\sum_{j=1}^{m}[j\perp k]\\ =\sum_{t}\mu(t)\lfloor\frac{n}{t}\rfloor \sum_{t|j}^{m}[j\perp k]\\ =\sum_{t}\mu(t)\lfloor\frac{n}{t}\rfloor \sum_{j=1}^{\lfloor\frac{m}{t}\rfloor}[tj\perp k]\\ =\sum_{t}\mu(t)\lfloor\frac{n}{t}\rfloor \sum_{j=1}^{\lfloor\frac{m}{t}\rfloor}[t\perp k][j\perp k]\\ =\sum_{t}\mu(t)\lfloor\frac{n}{t}\rfloor[t\perp k] \sum_{j=1}^{\lfloor\frac{m}{t}\rfloor}[j\perp k]\\ \]

\(g(x)=\sum_{i=1}^x\limits\ [i\perp k]\) ,发现 \(g(x)=\lfloor\dfrac{x}{k} \rfloor\varphi(k)+g(x\bmod k)\) ,可以 \(O(k)\) 预处理 ,\(O(1)\) 查询。

变为:

\[\sum_{t}\mu(t)\lfloor\frac{n}{t}\rfloor[t\perp k] \sum_{j=1}^{\lfloor\frac{m}{t}\rfloor}[j\perp k]=\sum_{t}g(\lfloor\frac{m}{t}\rfloor)\lfloor\frac{n}{t}\rfloor\mu(t)[t\perp k] \]

对于式子的前两项可以整除分块,需要快速求 \(\mu(t)[t\perp k]\) 的前缀和。

\(f(n)=\sum_{i=1}^n\limits \mu(i)[i\perp k]\) ,有:

\[f(n)=\sum_{i=1}^{n}[i\perp k]f(\lfloor\frac{n}{i}\rfloor)-\sum_{i=2}^ {n}[i\perp k]f(\lfloor\frac{n}{i}\rfloor) \]

发现 :

\[\sum_{i=1}^{n}[i\perp k]f(\lfloor\frac{n}{i}\rfloor)=\sum_{i=1}^{n}[i\perp k]\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\mu(j)[j\perp k]\\ =\sum_{i=1}^{n}\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\mu(j)[ij\perp k]\\ =\sum_{t=1}^n[t\perp k]\sum_{d|t}\mu(d)\\ =\sum_{t=1}^n[t\perp k][t==1]\\ =1 \]

进而:

\[f(n)=1-\sum_{i=2}^n[i\perp k]f(\lfloor\frac{n}{i}\rfloor) \]

发现 \(g(n)=\sum_{i=1}^n\limits [i\perp k]\) 前面已经可以 \(O(1)\) 求出,杜教筛即可。

点击查看代码
#include<bits/stdc++.h>
using namespace std;
int n,m,k;
const int lim=1e7;
int pri[10000005],top,mu[10000005],g[1005],f[10000005],phi[10000005];
bool vis[10000005];
int gcd(int x,int y){
	return y?gcd(y,x%y):x;
}
inline void init(){
	mu[1]=g[1]=f[1]=1;
	for(int i=2;i<=lim;i++){
		if(!vis[i])pri[++top]=i,mu[i]=-1,phi[i]=i-1;
		if(i<k)g[i]=g[i-1]+(gcd(i,k)==1);
		f[i]=f[i-1]+(gcd(i,k)==1)*mu[i];
		for(int j=1;j<=top&&i*pri[j]<=lim;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){
				phi[i*pri[j]]=phi[i]*pri[j];break;
			}
			mu[i*pri[j]]=-mu[i];phi[i*pri[j]]=phi[i]*phi[pri[j]];
		}
	}
}
inline long long G(int x){
	return 1ll*x/k*phi[k]+g[x%k];
}
map<int,long long> mp;
inline long long F(int x){
	if(x<=lim)return f[x];
	if(mp.count(x))return mp[x];
	long long res=1;
	for(int l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		res-=(G(r)-G(l-1))*F(x/l);
	}return mp[x]=res;
}
inline void solve(){
	long long res=0;
	for(int l=1,r;l<=n&&l<=m;l=r+1){
		r=min(n/(n/l),m/(m/l));
		res+=(F(r)-F(l-1))*G(m/l)*(n/l);
	}
	printf("%lld",res);
}
int main(){
	scanf("%d%d%d",&n,&m,&k);
	init();solve();


	return 0;
}



SCZ 2022-3-17 密码

求 :

\[\sum_{i=1}^n\sum_{j=1}^m \sigma_k(ij) \]

有:

\[\sigma_k (ij)=\sum_{x|i}\sum_{y|j}[x\perp y](\frac{i}{x}\cdot y)^k \]

进而:

\[\sum_{i=1}^n\sum_{j=1}^m \sigma_k(ij)=\sum_{i=1}^n\sum_{j=1}^m \sum_{x|i}\sum_{y|j}[x\perp y](\frac{i}{x}\cdot y)^k\\ =\sum_{i=1}^n\sum_{j=1}^m \sum_{x|i}\sum_{y|j}\sum_{d|x,d|y}\mu(d)(\frac{i}{x}\cdot y)^k\\ =\sum_{i=1}^n\sum_{x|i}(\frac{i}{x})^k\sum_{j=1}^m \sum_{y|j}(y)^k\sum_{d|x,d|y}\mu(d)\\ =\sum_{d}\mu(d)\sum_{i=1}^n\sum_{d|x,x|i}(\frac{i}{x})^k\sum_{j=1}^m \sum_{d|y,y|j}(y)^k\\ =\sum_{d}\mu(d)d^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sigma_k(i)\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sigma_k(j) \]

考虑构造 \(\sigma_k\)\(\mu\cdot id^k\) 的杜教筛。

发现 \(\sigma_k*\mu=id^k\) ,有:

\[\sum_{i=1}^n\sigma_k(i)=\sum_{i=1}^ni^k-\sum_{i=2}^n\mu(i)\times\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\sigma_k(j) \]

要筛 \(\mu\) 的前缀和:

\[\sum_{i=1}^{n}\mu(i)=1-\sum_{i=2}^{n}1\times\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\mu(j) \]

再看 \(\mu\cdot id^k\) ,发现 \((\mu\cdot id^k)*id^k=1\) ,因此:

\[\sum_{i=1}^{n}\mu(i)\cdot i^k=1-\sum_{i=2}^{n}i^k\times\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\mu(j)\cdot j^k \]

其中自然数 \(k\) 次幂前缀和是 \(k+1\) 次多项式,可以拉格朗日插值或高斯消元求出。其余部分杜教筛即可。

点击查看代码

Min_25 筛

对于满足一下条件的积性函数 \(𝑓\) ,如果其满足一下两条件,那么可以快速计算 \(𝐹\)

  • \(𝑓 (𝑝)\) 是关于 \(𝑝\) 的低阶多项式。

  • \(𝑓(𝑝^𝑘)\) 有快速计算的方法。

比如 \(𝜑 (𝑝) = 𝑝 − 1\)\(𝜑 (𝑝^𝑘) = (𝑝 − 1) 𝑝^{𝑘−1} = 𝜑 (𝑝^{𝑘−1})∗ 𝑝\)

时间复杂度为 \(O(\frac{n^{\frac{3}{4}}}{log\ n})\),空间复杂度为 \(O(\sqrt n)\)

step 1:

对于每一个 \(x=\lfloor \frac{n}{i} \rfloor\),求出:

\[\sum_{i=1}^{x}[i \in prime ]\cdot i^k \]

令:

\[G(n,i)=\sum_{j=1}^{n}[j \in prime | div(j)>pri_i]\cdot i^k \]

其中 \(div(j)\) 表示 \(j\) 的最小质因子。

\(G(n,i)\) 即为 \(≤n\) 且在埃氏筛第 \(i\) 轮后未被筛去的数的 \(k\) 次方和。

一个合数 \(x\ (x≤n)\) ,一定含有 \(≤\sqrt n\) 的质因子,因此若令 \(cnt\) 表示 \(≤n\) 的质数个数,\(G(n,cnt)\) 即为所求。

考虑递推求 \(G(n,i)\)

\(pri_i^2>n\) ,则第 \(i\) 轮不会删除任何数 ,有 \(G(n,i)=G(n,i-1)\)

\(pri_i^2≤n\),则第 \(i\) 轮会删除的数的 \(k\) 次方和即为 \(pri_i^k\times (G(\lfloor \frac {n}{pri_i} \rfloor,i-1)-\sum_{j=1}^{i-1}\limits pri_j^k)\)

因此,有:

\[G(n,i)=\begin{cases} \ \sum_{j=1}^{n}\limits j^k \\ G(n,i-1) \\ G(n,i-1)-pri_i^k\times (G(\lfloor \frac {n}{pri_i} \rfloor,i-1)-\sum_{j=1}^{i-1}\limits pri_j^k) \end{cases}\ \ \ \begin{matrix} (i=0) \\ (pri_i^2>n) \\ (pri_i^2≤n) \end{matrix} \]

在储存时,由于 \(\lfloor \frac{n}{i} \rfloor\) 的种类数不超过 \(2\times \sqrt n\)

因此可以开两个数组 \(ind1,ind2\) 将这 \(2\times \sqrt n\) 个数离散化,若 \(\lfloor \frac{n}{i} \rfloor ≤ \sqrt n\) 则将其存在 \(ind1[\lfloor \frac{n}{i} \rfloor]\) 处,否则将其存在 \(ind2[\lfloor \frac{n}{\lfloor \frac{n}{i} \rfloor} \rfloor]\) 处。

最终可以递推出 \(G(n,cnt)\) ,代码大致如下:

ll r=0;
for(ll l = 1; l <= n; l = r + 1) {
	r = n / (n / l); w[++tot]= n / l;
	g[tot] = n / l;
	g[tot] = (g[tot] * (g[tot] + 1) * inv2 - 1);
	if( n / l <= B ) ind1[n / l] = tot;
	else ind2[n / (n / l)] = tot;
}
for(int i = 1; i <= top; i++) {
	for(int j = 1; j <= tot && pri[i] * pri[i] <= w[j]; j++) {
		int k = w[j] / pri[i] <= B ? ind1[w[j] / pri[i]] : ind2[n / (w[j] / pri[i])];
		g[j] = (g[j] - pri[i] * (g[k] - sum[i - 1]));
	}
}

step 2:

定义 \(S(n,i)=\sum_{j=1}^{n}[div(j)>=pri_i]\cdot f(j)\),显然,最后答案就是 \(S(n,1)+f(1)\)

其中,对于 \(S(n,i)\) 质数的贡献为 \(\sum_{k}\limits G^k(n)\cdot [x^k]f(i)-\sum_{j=1}^{i-1}\limits pri_j\)

对于合数,我们枚举它的最小质因子,以及最小质因子的幂次,贡献为 \(\sum_{j=i}^{pri_j^e≤n}\limits f(pri_j^e)\cdot (S(\lfloor \frac{n}{pri_i^e} \rfloor,j+1)+[e!=1])\)

因此,有:

\[S(n,i)=\begin{cases} 0 \\ \sum_{k}\limits G^k(n)\cdot [x^k]f(i)-\sum_{j=1}^{i-1}\limits pri_j-\sum_{j=i}^{pri_j^e≤n}\limits f(pri_j^e)\cdot (S(\lfloor \frac{n}{pri_i^e} \rfloor,j+1)+[e!=1])\end{cases}\ \ \ \begin{matrix}(pri_i>n)\\(pri_i≤n)\end{matrix} \]

简计:

\[S(n,i)=\begin{cases}0\\g[n]-sum[i-1]-\sum_{j=i}^{pri_j^{e+1}≤n}\limits (f(pri_j^e)\cdot S(\lfloor \frac{n}{pri_i^e}\rfloor,j+1)+f(pri_j^{e+1}))\end{cases}\ \ \ \begin{matrix}(pri_i>n)\\ (pri_i≤n)\end{matrix} \]

ll S(ll x,int i){
	if( pri[i] > x ) return 0;
	int k = ( x <= B ? ind1[x] : ind2[n / x] );
	ll res = g[k] - sum[i-1];
	for(int j = i; j <= top && pri[j] * pri[j] <= x; j ++ ){
		ll P = pri[j];
		for(ll t = 1; P * pri[j] <= x; t ++ , P *= pri[j] ){
			ll Pw = P;
			res = res + Pw * S(x / P,j + 1);
			res = res + Pw * pri[j];
		}
	}return res;
}

时间复杂度 \(O(\frac{n^{\frac{3}{4}}}{log\ n})\)\(O(n^{1-ϵ})\) ,不需要记忆化。


\[\]

例题:

【模板】Min_25筛

模板题

点击查看代码 ```cpp ```

#6235. 区间素数个数

点击查看代码 ```cpp ```

SCZ 2022-2-18 取数

对于正整数 \(x =∏_{1≤i≤m}\limits p^{s_i}_i\) ,其中 \(p_i\) 为质数,\(p_i < p_{i+1}\)\(s_i > 0\)

令序列 \(A\)\(A_i = p^{s_i}_i + s^{p_i}_i\)

\(f(x) = ∑_P\limits ∏A_i[i ∈ P]\),其中 \(P\)\(i|1 ≤ i ≤ m\) 的非空子集,且 \(P\) 中不包含相邻的元素,即 \(i\)\(i + 1\) 不同时属于 \(P\)

\(∑^n_{i=2}\limits f(x)\)\(P\) 取模的结果。

发现对于单个数 \(x\) 可以通过 \(dp\) 求出 \(f(x)\) 的值,而且是积性函数,可以转换成矩阵做 \(Min_25\) ,也可以开结构体直接做。

点击查看代码
#include<bits/stdc++.h>
using namespace std;
#define ll long long
int lim;
ll n,md;
struct mat{
	ll tmp0,tmp1;
	mat(){
		tmp0=0;tmp1=0;
	}
	mat(ll _x,ll _y){
		tmp0=_x%md;tmp1=_y%md;
	}
	inline mat operator *(mat b){
		return mat((b.tmp0+b.tmp1)%md*tmp0%md,b.tmp0*tmp1%md);
	}
	inline mat operator +(mat b){
		return mat((tmp0+b.tmp0)%md,(tmp1+b.tmp1)%md);
	}
	inline mat operator -(mat b){
		return mat((tmp0-b.tmp0+md)%md,(tmp1-b.tmp1+md)%md);
	}
};
int pri[100005],top;
bool vis[100005];
int ind1[100005],ind2[100005],cnt;
ll w[200005],g0[200005],g1[200005],sum[100005];
inline void Init(){
	for(int i=2;i<=lim;i++){
		if(!vis[i])pri[++top]=i,sum[top]=(sum[top-1]+i)%md;
		for(int j=1;j<=top&&1ll*i*pri[j]<=lim;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0)break;
		}
	}
	for(ll l=1,r;l<=n;l=r+1){
		r=n/(n/l);w[++cnt]=n/l;
		g0[cnt]=(w[cnt]-1)%md;g1[cnt]=((w[cnt]+1)%md*(w[cnt]%md)/2-1)%md;
		if(w[cnt]<=lim)ind1[w[cnt]]=cnt;
		else ind2[n/w[cnt]]=cnt;
	}
	for(int i=1;i<=top;i++){
		for(int j=1;j<=cnt&&1ll*pri[i]*pri[i]<=w[j];j++){
			ll k=w[j]/pri[i]<=lim?ind1[w[j]/pri[i]]:ind2[n/(w[j]/pri[i])];
			g0[j]=((g0[j]-(g0[k]-i+1))%md+md)%md;
			g1[j]=((g1[j]-pri[i]*(g1[k]-sum[i-1]))%md+md)%md;
		}
	}
}
inline ll pwr(ll x,ll y){
	ll res=1;
	while(y){
		if(y&1)res=res*x%md;
		x=x*x%md;y>>=1;
	}
	return res;
}
mat Min_25(ll x,int i){
	if(x<=1||pri[i]>x)return mat();
	int k=x<=lim?ind1[x]:ind2[n/x];
	mat res=mat(g0[k],(g1[k]+g0[k])%md)-mat(i-1,(sum[i-1]+i-1)%md);
	for(int j=i;j<=top&&1ll*pri[j]*pri[j]<=x;j++){
		ll P=pri[j];
		for(int e=1;P*pri[j]<=x;e++,P*=pri[j]){
			res=res+mat(1,(P+pwr(e,pri[j]))%md)*Min_25(x/P,j+1);	
			res=res+mat(1,(P*pri[j]%md+pwr(e+1,pri[j]))%md);
		}
	}
	return res;
}
int main(){
	scanf("%lld%lld",&n,&md);lim=sqrt(n)+1; 
	Init();
	mat res=Min_25(n,1);
	printf("%lld\n",((res.tmp0+res.tmp1-n+1)%md+md)%md);

	return 0;
}


#6053. 简单的函数

点击查看代码 ```cpp ```

#188. 【UR #13】Sanrd

点击查看代码 ```
</details>
posted @ 2022-01-08 18:15  一粒夸克  阅读(251)  评论(0编辑  收藏  举报