min_25筛 学习小记

2 篇文章 0 订阅

终于在9102年之前搞完了这个东西。。

关于min_25筛,一种常数和写法优于洲阁筛的神奇筛法,复杂度大概是 O ( n 3 4 log ⁡ n ) O\left(\frac{n^{\frac{3}{4}}}{\log \sqrt{n}}\right) O(logn n43)的。我们可以利用它求一些积性函数的前缀和,要求 p p p为质数时 f ( p ) f(p) f(p) f ( p k ) f(p^k) f(pk)比较好求,并且可以做到1e10这样子

举个栗子,假设我们要求 ∑ i = 1 n F ( i ) \sum\limits_{i=1}^n F(i) i=1nF(i) F ( x ) = x F(x)=x F(x)=x。考虑把 x x x分成质数、合数和 1 1 1三类,我们先讨论对 x x x为质数时求和

g ( i , j ) = ∑ k = 1 i F ( k ) [ k ∈ P 或 m n ( k ) > P j ] g(i,j)=\sum\limits_{k=1}^i {F(k)[k\in P或mn(k)>P_j]} g(i,j)=k=1iF(k)[kPmn(k)>Pj],其中 P P P是质数集合, m n ( x ) mn(x) mn(x)表示 x x x的最小质因子, P i P_i Pi表示第i大的质数
形象地,我们可以把 g ( i , j ) g(i,j) g(i,j)看成是筛素数筛了 j j j轮后剩余的 F ( ) F() F()的和,并且十分显然地 g ( n , ∣ P ∣ ) g(n,|P|) g(n,P)就是所有素数的 F ( ) F() F()的和。注意我们这里实际上是把剩余的所有数字都视作质数带入仅有质数满足的柿子在算了,因此这里算出来的并不是真正的答案

考虑怎么求 g ( i , j ) g(i,j) g(i,j),我们只需要考虑第j次删去了哪些数字,显然就是最小质因子恰好为 P j P_j Pj的数。因此 g ( i , j ) = g ( i , j − 1 ) − f ( P j ) ⋅ [ g ( ⌊ i P j ⌋ , j − 1 ) − g ( P j , j − 1 ) ] g(i,j)=g(i,j-1)-f(P_j)\cdot[g(\lfloor\frac{i}{P_j}\rfloor,j-1)-g(P_j,j-1)] g(i,j)=g(i,j1)f(Pj)[g(Pji,j1)g(Pj,j1)],注意还要加上那些本身就是质数的贡献,并且当 P j 2 > i {P_j}^2> i Pj2>i的时候是没有后面的贡献的

显然 m n ( n ) ≤ n mn(n)\le \sqrt n mn(n)n ,因此我们只需要线性筛出 n \sqrt n n 个质数。并且第一维只有 n \sqrt n n 个数有用

再考虑怎么求剩余部分。记 S ( i , j ) = ∑ k = 1 i F ( k ) [ k ∈ P 或 m n ( k ) ≥ P j ] S(i,j)=\sum\limits_{k=1}^i {F(k)[k\in P 或mn(k)\ge P_j]} S(i,j)=k=1iF(k)[kPmn(k)Pj],那么 S ( n , 1 ) + F ( 1 ) = ∑ i = 1 n F ( i ) S(n,1)+F(1)=\sum\limits_{i=1}^n F(i) S(n,1)+F(1)=i=1nF(i)

由F(x)为积性函数可知,我们只需要枚举一个最小质数 P k P_k Pk和指数 e e e,考虑加上 P k e {P_k}^e Pke的贡献即可,因此 S ( i , j ) = ( g ( i , ∣ P ∣ ) − g ( P j − 1 , j − 1 ) ) + ∑ k = j ∣ P ∣ ∑ e ≥ 1 且 P k e + 1 ≤ i F ( P k e ) S ( ⌊ i P k ⌋ , k + 1 ) + F ( P k e + 1 ) S(i,j)=\left(g(i,|P|)-g(P_{j-1},j-1)\right)+\sum\limits_{k=j}^{|P|}\sum\limits_{e\ge 1且{P_k}^{e+1}\le i}{F({P_k}^e)S(\lfloor\frac{i}{P_k}\rfloor,k+1)+F({P_k}^{e+1})} S(i,j)=(g(i,P)g(Pj1,j1))+k=jPe1Pke+1iF(Pke)S(Pki,k+1)+F(Pke+1),柿子的前半部分是质数的答案,这里查了很多个blog好像都有点问题,わからない

具体实现的话可以预处理所有可能的根号取值并离散,然后求出g。对于S的计算递归即可

模板

loj6053
某一天,你发现了一个神奇的函数 f ( x ) f(x) f(x),它满足很多神奇的性质:

f ( 1 ) = 1 f(1)=1 f(1)=1
f ( p c ) = p ⊕ c f(p^c)=p \oplus c f(pc)=pc p p p为质数, ⊕ \oplus 表示异或)。
f ( a b ) = f ( a ) ⋅ f ( b ) f(ab)=f(a)\cdot f(b) f(ab)=f(a)f(b) a a a b b b 互质)。

你看到这个函数之后十分高兴,于是就想要求出 ∑ i = 1 n f ( i ) \sum\limits_{i=1}^n f(i) i=1nf(i)

由于这个数比较大,你只需要输出 ∑ i = 1 n f ( i )   m o d   ( 1 0 9 + 7 ) \sum\limits_{i=1}^n f(i) \bmod (10^9+7) i=1nf(i)mod(109+7)
n ≤ 1 0 10 n\le 10^{10} n1010

Solution


由于除2外所有质数都是奇数,因此当 x ∈ P 且 x ≠ 2 x\in P且x\ne2 xPx̸=2 x ⊕ 1 = x − 1 x\oplus 1=x-1 x1=x1
我们定义 g ( x ) = x [ x ∈ P ] g(x)=x[x\in P] g(x)=x[xP] h ( x ) = [ x ∈ P ] h(x)=[x\in P] h(x)=[xP],然后令 f ( x ) = g ( x ) − h ( x ) f(x)=g(x)-h(x) f(x)=g(x)h(x)套进上面直接算就可以了

Code


#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
#define rep(i,st,ed) for (register LL i=st;i<=ed;++i)

typedef long long LL;
const int MOD=1000000007;
const LL ny2=500000004;
const int N=200005;

int id1[N],id2[N],B,m;

bool not_prime[N];

LL f[N],g[N],sum[N],prime[N],w[N],n;

void pre_work(int n) {
	rep(i,2,n) {
		if (!not_prime[i]) prime[++prime[0]]=i,sum[prime[0]]=(sum[prime[0]-1]+i)%MOD;
		for (int j=1;i*prime[j]<=n&&j<=prime[0];++j) {
			not_prime[i*prime[j]]=1;
			if (i%prime[j]==0) break;
		}
	}
}

LL solve(LL x,LL y) {
	if (x<=1||prime[y]>x) return 0;
	LL res=(y==1)*2;
	int k=(x<=B)?(id1[x]):(id2[n/x]);
	res=(res+(f[k]-g[k])-(sum[y-1]-(y-1))+MOD)%MOD;
	for (LL i=y;prime[i]*prime[i]<=x&&i<=prime[0];++i) {
		LL p1=prime[i],p2=p1*p1;
		for (LL e=1;p2<=x;++e,p1=p2,p2*=prime[i]) {
			res=(res+solve(x/p1,i+1)*(prime[i]^e)%MOD+(prime[i]^(e+1)))%MOD;
		}
	}
	return (res%MOD+MOD)%MOD;
}

int main(void) {
	freopen("data.in","r",stdin);
	scanf("%lld",&n); B=sqrt(n);
	pre_work(B);
	for (LL i=1,j;i<=n;i=j+1) {
		w[++m]=n/i; j=n/w[m];
		(w[m]<=B)?(id1[w[m]]=m):(id2[j]=m);
		g[m]=w[m]-1; f[m]=(1LL*(w[m]+2)%MOD)*((w[m]-1)%MOD)%MOD*ny2%MOD;
	}
	rep(j,1,prime[0]) {
		for (int i=1;i<=m&&prime[j]*prime[j]<=w[i];++i) {
			int k=((w[i]/prime[j])<=B)?(id1[w[i]/prime[j]]):(id2[n/(w[i]/prime[j])]);
			f[i]=(f[i]-prime[j]*(f[k]-sum[j-1])%MOD+MOD)%MOD;
			g[i]=((g[i]-(g[k]-(j-1)))%MOD+MOD)%MOD;
		}
	}
	printf("%lld\n", solve(n,1)+1);
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值