min-25筛 简记

min-25 筛

part1

求 \(\sum\limits_{p\in prime}^{p\le n}g(p)\)

其中 \(g(x)=x^k\)

设 \(G_{n,k}=\sum\limits_{p\in prime\ or\ minp\gt prime_k}^{p\le n}g(p)\)

其中 \(minp\) 表示 \(p\) 的最小质因子.

设 \(sg(n)\) 表示前 \(n\) 个质数处 \(g(p)\) 之和

则有 \(G_{n,j}=G_{n,j-1}-prime_j^k(G_{\frac{n}{prime_j},j-1}-sg_j)\)

其边界为 \(G_{n,0}=\sum\limits_{i=2}^ni^k\)

part2

求 \(\sum\limits_{i=1}^{n}f(x)\)

其中 \(f(x)\) 为一积性函数且满足 \(f(p)=poly(p)\) \((p\in prime)\)

设 \(F_{n,k}=\sum\limits_{p\in prime\ or\ minp\gt prime_k}^{p\le n} f(p)\)

其中 \(minp\) 表示 \(p\) 的最小质因子.

设 \(sp(n)\) 表示前 \(n\) 个质数处 \(f(p)\) 之和

\(F_{n,k}=G_{n}-sp(k)+\sum\limits_{j=k+1}\sum\limits_{{prime_j}^e\le n}f({prime_j}^e)(F_{\frac{n}{{prime_j}^e},j}+[e\neq 1])\)

其中 \(G_{n}=\sum\limits_{p\in prime} poly(p)\)

放一个0-index的实现

void solve(){
	i64 N;
    cin>>N;
    i64 sqN=std::sqrt(N);
	auto sieve=[&](int n){
		vec<int> prime;
		vec<bool> flag(n+1);
		vec<Z> sum1(1),sum2(1);
		for(int i=2;i<=n;++i){
			if(flag[i]==false){
				prime.emplace_back(i);
				sum1.emplace_back(sum1.back()+Z(i));
				sum2.emplace_back(sum2.back()+Z(i)*i);
			}
			for(auto p:prime){
				if(p*i>n){
					break;
				}
				flag[p*i]=true;
				if(i%p==0){
					break;
				}
			}
		}
		return std::make_tuple(prime,sum1,sum2);
	};
	auto [prime,sum1,sum2]=sieve(sqN);
	int tot=1,num=prime.size();
	vec<i64> w(2*sqN+2);
	vec<Z> g1(2*sqN+2),g2(2*sqN+2);
	constexpr Z iv2=1_z/2,iv3=1_z/3;
	for(i64 l=1,r;l<=N;l=r+1){
		i64 x=(N/l);Z z=Z::get(x);
		r=(N/x),w[tot]=x,g1[tot]=z*(z+1)*iv2-1,g2[tot]=z*(z+iv2)*(z+1)*iv3-1,++tot;
	}
	auto get_id=[&](i64 x){
		return x<=sqN?(tot-x):(N/x);
	};
	for(int i=0;i<int(prime.size());++i){
		for(int j=1;i64(prime[i])*prime[i]<=w[j];++j){
			i64 k=get_id(w[j]/prime[i]);
			g1[j]-=(g1[k]-sum1[i])*prime[i];
			g2[j]-=(g2[k]-sum2[i])*prime[i]*prime[i];
		}
	}
	//mind>prime[y]
	auto get_sum=[&](auto&&self,i64 x,int y)->Z {
		//if x <= prime[y] return 0z;
		i64 k=get_id(x);
		Z res=(g2[k]-g1[k])-(sum2[y+1]-sum1[y+1]);
		for(int i=y+1;i<num&&prime[i]*prime[i]<=x;++i){
			i64 pj=prime[i],x_p=x/prime[i];
			int j=1;
			for(;prime[i]<x_p;++j,pj*=prime[i],x_p/=prime[i]){
				Z z=Z::get(pj);
				res+=z*(z-1)*(self(self,x_p,i)+(j!=1));
			}
			for(;pj<=x;++j,pj*=prime[i]){
				Z z=Z::get(pj);
				res+=z*(z-1)*(j!=1);
			}
		}
		return res;
	};
	cout<<get_sum(get_sum,N,-1)+1;
}