[LOJ572] Misaka Network 与求和 [杜教筛+Min25]

最新推荐文章于 2021-02-10 20:25:59 发布

FSYo

最新推荐文章于 2021-02-10 20:25:59 发布

阅读量378

点赞数 1

分类专栏：筛法

本文链接：https://blog.csdn.net/sslz_fsy/article/details/99192722

版权

筛法专栏收录该内容

9 篇文章 0 订阅

订阅专栏

传送门

$ans=\sum_{d=1}^n f(d)^k\sum_{i=1}^{n}\sum_{j=1}^n [(i,j)=d]=\sum_{d=1}^nf(d)^k(2*\sum_{i=1}^{n/d}\varphi(i)-1)=\sum_{d=1}^nf(d)^kg(n/d)$

首先 g 可以整除分块，然后杜教筛处理 g

考虑如何用 Min25 筛出 f，一个 d 有 f(d) ^ k 的贡献，当且仅当有一个 >= 它的质数在它前面

所有的方案数就是 >= 它的质数个数

于是 Min25 先处理出质数个数，然后加上合数贡献时算一下 f(d)

具体一点就是 $S(n, i)=\sum_{j=i}^{tot} \sum_{e=1}^{p_j^{e+1}<=n}S(n/p_j^e, j+1)+p_j^k*\sum_{i=p_j}^{n/p_j^e}[i\in P]$

#include<bits/stdc++.h>
#define N 1000050
using namespace std;
typedef long long ll;
typedef unsigned int ui;
ll n, m, k, S, tot, prim[N], w[N], id1[N], id2[N]; bool isp[N];
ui f[N], phi[N], fk[N];
ll Id(ll x){ return x <= S ? id1[x] : id2[n/x];}	
ui power(ui a, ll b){ ui ans = 1;
	for(;b;b>>=1){ if(b&1) ans *= a; a *= a;}
	return ans;
} 
void prework(){
	phi[1] = 1;
	for(int i = 2; i <= N-50; i++){
		if(!isp[i]) prim[++tot] = i, fk[tot] = power(i, k), phi[i] = i-1;
		for(int j = 1; j <= tot; j++){
			if(i * prim[j] > N-50) break;
			isp[i * prim[j]] = 1; if(i % prim[j] == 0){ phi[i * prim[j]] = phi[i] * prim[j]; break;}
			phi[i * prim[j]] = phi[i] * phi[prim[j]];
		}
	}
	for(int i = 2; i <= N-50; i++) phi[i] += phi[i-1];
}
ll Solve(ll x, int p){
	if(prim[p] > x) return 0; ll ans = 0;
	for(int i = p; i <= tot && prim[i] * prim[i] <= x; i++){
		for(ll e = 1, l = prim[i]; l * prim[i] <= x; e++, l *= prim[i]){
			ans += Solve(x / l, i + 1) + fk[i] * (f[Id(x / l)] - (i-1));
		}
	} return ans;
}
ui Phi[N];
ui Getphi(ll x){
	if(x <= N-50) return phi[x];
	int k = Id(x);
	if(Phi[k]) return Phi[k];
	ui ans = 1ll * x * (x + 1) / 2;
	for(ll l = 2, r; l <= x; l = r+1){
		ll v = x/l; r = x/v;
		ans -= (r-l+1) * Getphi(v);
	} return Phi[k] = ans;
} 
int main(){
	scanf("%lld%lld", &n, &k); m = 0; S = sqrt(n); prework();
	for(ll l = 1, r; l <= n; l = r+1){
		ll v = n/l; r = n/v; 
		if(v <= S) id1[v] = ++m; else id2[n/v] = ++m;
		w[m] = v; f[m] = w[m]-1;
	}
	for(int i = 1; i <= tot; i++){
		if(prim[i] * prim[i] > n) break;
		for(int j = 1; j <= m && prim[i] * prim[i] <= w[j]; j++){
			f[j] -= f[Id(w[j]/prim[i])] - (i - 1);
		}
	} 
	ui ans = 0;
	for(ll l = 1, r, lst = 0, nxt = 0; l <= n; l = r+1){
		ui v = n/l; r = n/v;
		ui v1 = 2 * Getphi(v) - 1;
		nxt = Solve(r, 1) + f[Id(r)];
		ans += (nxt - lst) * v1;
		lst = nxt;
	} cout << ans; return 0;
}