差别

这里会显示出您选择的修订版和当前版本之间的差别。

--- 2020-2021:teams:i_dont_know_png:potassium:sieve [2020/05/26 05:46]
potassium [实例]
+++ 2020-2021:teams:i_dont_know_png:potassium:sieve [2020/06/07 12:46] (当前版本)
nikkukun fix some typos
@@ 行 450: / 行 450: @@
 $$
-\begin{aligned}\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma(ij)&=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{x\mid i}\sum_{y\mid j}\frac{iy}{x}[(x,y)=1]\\&=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{x\mid i}\sum_{y\mid j}\frac{iy}{x}\sum_{d|x,d|y}\mu(d)\\&=\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{x\mid i,d\mid x}\sum_{y\mid j,d\mid y}\frac{iy}{x}\\&=\sum_{d=1}^{n}\mu(d)\sum_{d\mid x}\sum_{d\mid y}\frac yx\sum_{x\mid i}\sum_{y\mid j}i\\&=\sum_{d=1}^{n}\mu(d)\sum_{x=1}^{\lfloor\frac nd\rfloor}\sum_{y=1}^{\lfloor\frac nd\rfloor}\frac{yd}{xd}\sum_{i=1}^{\lfloor\frac{n}{xd}\rfloor}ixd\sum_{j=1}^{\lfloor\frac{n}{yd}\rfloor}\\&=\sum_{d=1}^{n}d\mu(d)\sum_{x=1}^{\lfloor\frac nd\rfloor}\sum_{i=1}^{\lfloor\frac{n}{xd}\rfloor}i\sum_{y=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac {n}{yd}\rfloor}y\\&=\sum_{d=1}^{n}d\mu(d)\sum_{x=1}^{\lfloor\frac nd\rfloor}\sum_{i=1}^{\lfloor\frac{n}{xd}\rfloor}i     \sum_{j=1}^{\lfloor\frac {n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac n{dj}\rfloor}y\\&=\sum_{d=1}^{n}d\mu(d)S_1(\lfloor\frac nd\rfloor)^2\\其中S_1(n)&=\frac 12(1+n)n\\\end{aligned}
+\begin{aligned}
+\sum_{i=1}^{n}\sum_{j=1}^{n}\sigma(ij)&=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{x\mid i}\sum_{y\mid j}\frac{iy}{x}[(x,y)=1]\\
+&=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{x\mid i}\sum_{y\mid j}\frac{iy}{x}\sum_{d|x,d|y}\mu(d)\\
+&=\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{x\mid i,d\mid x}\sum_{y\mid j,d\mid y}\frac{iy}{x}\\
+&=\sum_{d=1}^{n}\mu(d)\sum_{d\mid x}\sum_{d\mid y}\frac yx\sum_{x\mid i}\sum_{y\mid j}i\\
+&=\sum_{d=1}^{n}\mu(d)\sum_{x=1}^{\lfloor\frac nd\rfloor}\sum_{y=1}^{\lfloor\frac nd\rfloor}\frac{yd}{xd}\sum_{i=1}^{\lfloor\frac{n}{xd}\rfloor}ixd\sum_{j=1}^{\lfloor\frac{n}{yd}\rfloor}\\
+&=\sum_{d=1}^{n}d\mu(d)\sum_{x=1}^{\lfloor\frac nd\rfloor}\sum_{i=1}^{\lfloor\frac{n}{xd}\rfloor}i\sum_{y=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac {n}{yd}\rfloor}y\\
+&=\sum_{d=1}^{n}d\mu(d)\sum_{x=1}^{\lfloor\frac nd\rfloor}\sum_{i=1}^{\lfloor\frac{n}{xd}\rfloor}i     \sum_{j=1}^{\lfloor\frac {n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac n{dj}\rfloor}y\\
+&=\sum_{d=1}^{n}d\mu(d)(\sum_{i=1}^{\lfloor\frac nd\rfloor}S_1(\lfloor\frac n{id}\rfloor))^2\\
+其中S_1(n)&=\frac 12(1+n)n
+\\
+\end{aligned}
 $$
-于是可以对 $d$ 数论分块，只需要快速计算出 $f=\mathrm{id}\cdot \mu$ 的前缀和即可，令 $g=\mathrm{id}$ ，显有 $f\ast g=\epsilon$。
+于是可以对 $d$ 数论分块， $\sum_{i=1}^{n}\frac 12\lfloor\frac ni\rfloor(\lfloor\frac ni\rfloor+1)$ 线筛预处理一部分，剩余 $O(\sqrt n)$ 直接算（或者不处理直接暴力， $O(n^{\frac 34})$ ）；另外需要快速计算出 $f=\mathrm{id}\cdot \mu$ 的前缀和，令 $g=\mathrm{id}$ ，显有 $f\ast g=\epsilon$。
+从上述步骤的第六步到第七步的转换，验证了等式 $\sum_{i=1}^{n}\frac 12\lfloor\frac ni\rfloor(\lfloor\frac ni\rfloor+1)=\sum_{i=1}^{n}i\lfloor\frac ni\rfloor$ 。而后者可以视作在 $i$ 处统计了 $i$ 的所有因数之和（对于每一个因数 $i$ ，在他共 $\lfloor\frac ni\rfloor$ 个整数倍的位置都被统计一次），即为 $\sum_{i=1}^{n}\sigma(i)$ 。实现的时候存一个 $div_i$ 表示 $i$ 的最小质因数做贡献的一项 $(1+p_m+p_m^2+\cdots+p_m^{k_m})$ 的值即可。如果不转化为 $\sigma$ 的前缀和，也可以直接打表验证积性函数后记录 $p^k$ 下的取值进行线筛。
+<hidden 参考实现>
+<code:cpp>
+#include<cstdio>
+#include<algorithm>
+#include<queue>
+#include<map>
+#include<cstring>
+#include<cmath>
+#include<cstdlib>
+#include<set>
+#include<unordered_map>
+#include<vector>
+typedef long long ll;
+using namespace std;
+#define pii pair<int,int>
+#define pb push_back
+#define mp make_pair
+#define fi first
+#define se second
+#define N 3000000
+#define mod 1000000007
+char isnp[N+10];
+int pri[N+10],cnt,mu[N+10],sm[N+10];
+ll sigma[N+10],dv[N+10],sr[N+10];
+void sieve(int n){
+	int i,j;
+	isnp[0]=isnp[1]=1;
+	mu[1]=1;sigma[1]=1;
+	for(i=2;i<=n;i++){
+		if(!isnp[i])pri[cnt++]=i,mu[i]=-1,sigma[i]=dv[i]=1+i;
+		for(j=0;j<cnt;j++){
+			if(pri[j]*i>n)break;
+			isnp[i*pri[j]]=1;
+			if(i%pri[j]==0){
+				mu[i*pri[j]]=0;
+				dv[i*pri[j]]=dv[i]*pri[j]+1;
+				sigma[i*pri[j]]=sigma[i]/dv[i]*dv[i*pri[j]];
+				break;
+			}else{
+				mu[i*pri[j]]=-mu[i];
+				sigma[i*pri[j]]=sigma[i]*sigma[pri[j]];
+				dv[i*pri[j]]=pri[j]+1;
+			}
+		}
+	}
+	for(i=1;i<=n;i++){
+		sm[i]=(sm[i-1]+i*mu[i])%mod;
+		sr[i]=(sr[i-1]+sigma[i])%mod;
+	}
+}
+ll inv2;
+ll qp(ll a,ll p){
+	ll ans=1;
+	for(;p;p>>=1,a=a*a%mod)if(p&1)ans=ans*a%mod;
+	return ans;
+}
+ll s1(int l,int r){
+	return (l+r)*(r-l+1LL)%mod*inv2%mod;
+}
+map<int,ll>m;
+ll getRight(int n){
+	if(n<=N)return sr[n];
+	int i,r;
+	ll res=0;
+	for(i=1;i<=n;i=r+1){
+		r=n/(n/i);
+		res+=n/i*(n/i+1LL)%mod*(r-i+1)%mod*inv2%mod;
+	}
+	return res%mod;
+}
+ll getS(int n){
+	if(n<=N)return sm[n];
+	if(m.count(n))return m[n];
+	int i,r;
+	ll res=1;
+	for(i=2;i<=n;i=r+1){
+		r=n/(n/i);
+		(res-=getS(n/i)*s1(i,r)%mod)%=mod;
+	}
+	return m[n]=res;
+}
+ll get(int n){
+	ll tmp,res=0;
+	int i,r;
+	for(i=1;i<=n;i=r+1){
+		r=n/(n/i);
+		tmp=getRight(n/i);
+		(res+=tmp*tmp%mod*(getS(r)-getS(i-1))%mod)%=mod;
+	}
+	return (res+mod)%mod;
+}
+int main(){
+	int n;
+	inv2=qp(2,mod-2);
+	sieve(N);
+	scanf("%d",&n);
+	printf("%lld",get(n));
+	return 0;
+}
+</code>
+</hidden>
+\\
@@ 行 501: / 行 618: @@
 ===== min_25 筛 =====
+min_25 筛可以用来解决具有下列性质的积性函数 $f$ 在单点处的前缀和：
+  - 存在完全积性函数 $f'$，使得质数 $p$ 处有 $f'(p) = f(p)$；
+  - $f'(p)$ 的前缀和很好计算；
+  - $f'(p^k)$ 的值很好算（是多项式或可以很快求出）。
+为了方便，我们约定：
+  * 若没有特别说明，下文中的所有 $p$ 都取质数；
+  * $\mathrm{pri}_j$ 表示第 $j$ 小的质数；
+  * $\mathrm{minp}(n)$ 表示 $n$ 的最小质因子。
+==== 计算质数贡献 ====
+我们先考虑解决求所有分块点处质数的贡献，即 $G(m) = \sum _{i=1}^m [i \text{ is prime}] f'(i)$，其中 $m$ 是 $n$ 的数论分块点。
+设状态 $g(j, m)$ 表示 $\sum _{i=1}^m [i \text{ is prime or } \mathrm{minp}(i) > \mathrm{pri}_j] f'(i)$，则我们要求的就是 $G(m) = g(\sum _{j \land \mathrm{pri}_j \leq m} 1, m)$ 了。这个状态可以理解为，Eratosthenes 筛法中用 $\mathrm{pri}_j$ 筛掉它的倍数后，剩余没被筛掉部分的和。
+显然，边界条件是 $g(0, m) = \sum _{i=2}^m f'(i)$。请注意一定要去掉 $f'(1)$，因为它不会被筛掉，我们只要最后加上即可。
+考虑转移。若 $\mathrm{pri}_j > \sqrt m$，此时筛不掉任何数，故 $g(j, m) = g(j-1, m)$，故考虑 $\mathrm{pri}_j \leq \sqrt m$ 的情况。此时 $\mathrm{pri}_j$ 能筛掉的**合数**都能表示为 $\mathrm{pri}_j$ 乘一个不超过 $\left\lfloor \dfrac m{\mathrm{pri}_j} \right\rfloor$ 的、最小质因子**大于等于** $\mathrm{pri}_j$ 的数，即
+$$
+g(j, m) = g(j - 1, m) - f'(\mathrm{pri}_j) \times [g(j - 1, \left\lfloor \dfrac m{\mathrm{pri}_j} \right\rfloor) - \sum _{i=1}^{j-1} f'(\mathrm{pri}_i)]
+$$
+因为 $f'$ 是完全积性函数，所以即使分解的东西不互质也可以直接乘在一起。另外，式子最后的部分用于减去质数的贡献（参考 $g$ 的定义）。
+具体实现时，可以用滚动数组存放 $g$，滚掉它的第一维 $j$。可以证明这一部分的时间复杂度是 $O\left(\dfrac {n^{3/4}}{\log n} \right)$ 的。
+==== 计算合数贡献 ====
+设 $S(n, j) = \sum _{i=1}^m [\mathrm{minp}(i) \geq \mathrm{pri}_j] f(i)$，注意中间对 $\mathrm{minp}(i)$ 的要求和 $G$ 有所不同。根据定义，我们最后要求的就是 $S(n, 1) + f(1)$。
+对于 $S(n, j)$，质数的贡献为 $G(n) - \sum _{i=1}^{j-1} f'(\mathrm{pri}_i)$，而合数的贡献需要枚举其最小质因子及其次数，分为两个互质的部分。转移如下：
+$$
+S(n, j) = G(n) - \sum _{i=1}^{j-1} f'(\mathrm{pri}_i) + \sum _{i = j}^{\mathrm{pri}_i \leq n} \sum _{e=1} ^{\mathrm{pri}_i ^{e+1} \leq n} [f(\mathrm{pri}_i ^e) S(\left\lfloor \dfrac n{\mathrm{pri}_i^e} \right\rfloor, i + 1) + f(\mathrm{pri}_i^{e+1})]
+$$
+这部分不需要记忆化，直接递归的时间复杂度是 $O(n^{1 - \varepsilon})$ 的，而且跑 $n = 10^{10}$ 的数据也很飞快。
@@ 行 506: / 行 668: @@
 ==== 实例 ====
+=== Luogu 5325 ===
+**题意**：积性函数 $f(n)$ 的质数幂次取值 $f(p^k) = p^k(p^k - 1)$，求前缀和。
+**题解**：$f(p^k) = (p^2)^k - p^k$，拆成两部分计算质数处贡献，最后合并即可。
+<hidden 参考实现>
+<code:cpp>
+#include <bits/stdc++.h>
+using namespace std;
+#define fi first
+#define se second
+#define lch (o << 1)
+#define rch (o << 1 | 1)
+typedef double db;
+typedef long long ll;
+typedef unsigned int ui;
+typedef pair<int, int> pint;
+const int N = 2e5 + 5;
+const int MOD = 1e9 + 7;
+const int INV2 = 500000004;
+const int INV6 = 166666668;
+const int INF = 0x3f3f3f3f;
+const ll INF_LL = 0x3f3f3f3f3f3f3f3f;
+// 注意两倍空间
+int idx[N];
+ll val[N];
+ll _n, sqrtN;
+int nIdx;
+int GetId(ll x) {
+	if (x <= sqrtN)
+		return idx[x];
+	else
+		return idx[sqrtN + _n / x];
+}
+void InitId(ll n) {
+	// x =  n / i 是该分块点中能取到最大 / 最右侧的值
+	for (ll i = 1; i <= n; i++) {
+		ll x = n / i;
+		val[++nIdx] = x;
+		if (x <= sqrtN)
+			idx[x] = nIdx;
+		else
+			idx[sqrtN + n / x] = nIdx;
+		i = n / (n / i);
+	}
+}
+int pri[N];
+ll s1[N], s2[N];
+bool notPri[N];
+void InitPri() {
+	for (int i = 2; i < N; i++) {
+		if (!notPri[i])
+			pri[++pri[0]] = i;
+		for (int j = 1; j <= pri[0]; j++) {
+			int p = pri[j];
+			if (i * p >= N) break;
+			notPri[i * p] = 1;
+			if (i % p == 0) break;
+		}
+	}
+	for (int i = 1; i <= pri[0]; i++) {
+		s1[i] = (s1[i - 1] + 1LL * pri[i] * pri[i]) % MOD;
+		s2[i] = (s2[i - 1] + pri[i]) % MOD;
+	}
+}
+ll f1[N], f2[N];
+// Sieve1(n) = sum [i = 1 -> n] [i is isprime] * f(i)
+// 传入的 n 是原始值
+void Sieve1(ll n) {
+	// 预处理初值
+	// 需要注意这里一般不计算 f(1)
+	for (int i = 1; i <= nIdx; i++) {
+		ll x = val[i] % MOD;
+		f1[i] = x * (x + 1) % MOD * (2 * x + 1) % MOD * INV6 % MOD - 1;
+		f2[i] = (1 + x) * x % MOD * INV2 % MOD - 1;
+	}
+	for (int j = 1; j <= pri[0]; j++) {
+		ll p = pri[j];
+		if (p * p > n) break;
+		for (int i = 1; i <= nIdx; i++) {
+			ll x = val[i];
+			if (p * p > x) break;
+			// f(m) = f(idx[m / pri[j]]) - sum[j - 1];
+			f1[i] -= (p * p) % MOD * (f1[GetId(x / p)] - s1[j - 1]);
+			f2[i] -= p * (f2[GetId(x / p)] - s2[j - 1]);
+			f1[i] %= MOD;
+			f2[i] %= MOD;
+		}
+	}
+}
+// f(p ^ e)
+ll F(ll p, int e) {
+	ll t = 1;
+	while (e--) t *= p;
+	t %= MOD;
+	return t * (t - 1) % MOD;
+}
+// Sieve2(n, j) = sum [i = 1 -> n] [minp(i) >= pri[j]] * f(i)
+// 传入的 n 是原始值
+ll Sieve2(ll n, int j) {
+	if (n <= 1 || pri[j] > n) return 0;
+	int id = GetId(n);
+	// ret = Sieve1(n) - sum [i = 1 -> j - 1] f(i) + ...
+	ll ret = (f1[id] - s1[j - 1]) - (f2[id] - s2[j - 1]);
+	// ret = ... + 合数部分
+	for (int i = j; i <= pri[0]; i++) {
+		ll p = pri[i], t = p; // t = p ^ e
+		if (p * p > n) break;
+		for (int e = 1; p * t <= n; t *= p, e++) {
+			ll tmp = F(p, e) * Sieve2(n / t, i + 1) + F(p, e + 1);
+			ret = (ret + tmp) % MOD;
+		}
+	}
+	return ret;
+}
+int main() {
+	ios::sync_with_stdio(0);
+	ll n;
+	cin >> n;
+	sqrtN = floor(sqrt(n));
+	_n = n;
+	InitId(n);
+	InitPri();
+	Sieve1(n);
+	ll ans = (Sieve2(n, 1) + 1) % MOD;
+	cout << (ans + MOD) % MOD;
+	return 0;
+}
+</code>
+</hidden>
+\\

CVBB ACM Team

用户工具

站点工具

差别

页面工具