[min25筛学习小记]LOJ6053

###min25筛
解决一类积性函数求前缀和问题,主旨为模拟普通筛法过程。
假设我们要求 ∑ i = 1.. n F ( i ) \sum_{i=1..n}F(i) i=1..nF(i)
一般来说,能做的题,F§可表示为p相关的多项式,其中p为质数,下面将以 F ( p ) = p k F(p)=p^k F(p)=pk为例。
也即统计 ∑ i i k \sum_i i^k iik
为了方便做,我们先做出所有质数的f和。怎么做呢?
g ( n , i ) = ∑ j = 2.. n F ( j ) [ j 为 质 数 或 不 含 小 于 等 于 p i 的 质 因 子 ] g(n,i)=\sum_{j=2..n}F(j)[j为质数或不含小于等于p_i的质因子] g(n,i)=j=2..nF(j)[jpi]
其中 p i p_i pi表示从小到大第i个质数。g的含义就是,你用带log的筛法去做,现在暂时排除掉前i个质数的倍数的f和。
和原始筛法一样,注意到,n以内的合数,必定存在 ≤ n ≤\sqrt n n 的质因子,那么i最大取到根号的位置,这时就把所有质数的F和求出来了。
转移:
g ( x , i ) = g ( x , i − 1 ) − p i k ( g ( ⌊ x p i ⌋ , i − 1 ) − s u m [ i − 1 ] ) g(x,i)=g(x,i-1)-p_i^k(g(\lfloor \frac{x}{p_i}\rfloor,i-1)-sum[i-1]) g(x,i)=g(x,i1)pik(g(pix,i1)sum[i1]).
而如果 p i 2 > x p_i^2>x pi2>x,后面部分可以略去。
其中sum[i]表示的是前i个质数的F和,这个可以预处理。
这个式子的含义是,在前i-1个质数筛去合数结果的基础上,再用 p i p_i pi筛,你把要筛的数除去 p i p_i pi,显然剩下不能有小于 p i p_i pi的质因子。g状态比较多,用递推会比较快。

再考虑怎么把所有的数都算上。
S ( x , i ) = ∑ j = 2.. n F ( i ) [ j 的 最 小 质 因 子 大 于 等 于 p i ] S(x,i)=\sum_{j=2..n}F(i)[j的最小质因子大于等于p_i] S(x,i)=j=2..nF(i)[jpi]
和g相反,我们要从大的i推到小的i。
质数部分,由g求出来。对于合数,枚举每个数的最小质因子及其次幂来转移:
S ( x , i ) = ( g ( x , m x ) − s u m [ i − 1 ] ) + ∑ p j ≥ p i ∑ e = 1.. l o g p j x F ( p j k ) S ( ⌊ x p j e ⌋ , j + 1 ) + F ( p j k + 1 ) S(x,i)=(g(x,mx)-sum[i-1])+\sum_{p_j≥p_i}\sum_{e=1..log_{p_j}x}F(p_j^k)S(\lfloor \frac{x}{p_j^e}\rfloor,j+1)+F(p_j^{k+1}) S(x,i)=(g(x,mx)sum[i1])+pjpie=1..logpjxF(pjk)S(pjex,j+1)+F(pjk+1)
mx为能取到的最大的质数。
然后玄学地,S直接转移不用记忆化,跑得非常快,1s可以过1e10的n;而g则需要递推来写,而不是哈希,才能过1e10。
###LOJ6053
考虑直接套公式即可。
###代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<map>
#include<set>
using namespace std;
typedef long long ll;
typedef double db;
#define fo(i,j,k) for(i=j;i<=k;i++)
#define fd(i,j,k) for(i=j;i>=k;i--)
#define cmax(a,b) (a=(a>b)?a:b)
#define cmin(a,b) (a=(a<b)?a:b)
const int N=2e5+5,M=1e6+5,mo=1e9+7;
ll pd[M],pri[M],spr[M],sprf[M];
ll val[N],g0[N],g1[N],id[2][N],pt,sn,n,nn,tt,i,j,v,upb,kan;
void predo(ll n)
{
	ll i,j,t;
	fo(i,2,n)
	{
		if (!pd[i])
			pri[++pri[0]]=i;
		fo(j,1,pri[0])
		{
			if (1ll*i*pri[j]>n) break;
			t=i*pri[j];
			pd[t]=1;
			if (i%pri[j]==0) break;
		}
	}
	fo(i,1,pri[0]) spr[i]=(spr[i-1]+pri[i])%mo;
	fo(i,1,pri[0]) sprf[i]=(sprf[i-1]+pri[i]+((i==1)?1:-1))%mo;
}
void sieve_g()
{
	sn=trunc(sqrt(n));
	tt=0;
	i=1;
	while (i<=n)
	{
		v=n/i;
		j=n/v;
		g0[++tt]=(v-1)%mo;
		g1[tt]=(v%mo*(v%mo+1ll)/2ll-1ll)%mo;
		val[tt]=v;
		if (v<=sn) id[0][v]=tt;else id[1][i]=tt;
		i=j+1;
	}
	fo(upb,1,pri[0]) if (!(n/pri[upb]/pri[upb])) break;
	upb--;
	fo(j,1,upb)
	{
		while (tt&&!(val[tt]/pri[j]/pri[j])) tt--;
		fo(i,1,tt)
		{
			pt=val[i]/pri[j];
			if (pt<=sn) pt=id[0][pt];else pt=id[1][n/pt];
			g0[i]=(g0[i]-g0[pt]+(j-1))%mo;
			g1[i]=(g1[i]-(g1[pt]-spr[j-1])*pri[j])%mo;
		}
	}
}
ll S(ll n,ll x)
{
	if (n<pri[x]||n<=1) return 0;
	if (n<=sn) pt=id[0][n];else pt=id[1][nn/n];
	ll ret=(g1[pt]-g0[pt]+2-sprf[x-1])%mo,prod;
	ll i,k;
	fo(i,x,pri[0])
	{
		if (!(n/pri[i]/pri[i])) break;
		prod=pri[i];
		fo(k,1,200)
		{
			if (!(n/prod/pri[i])) break;
			ret=(ret+(kan=S(n/prod,i+1))*(pri[i]^k)+(pri[i]^(k+1)))%mo;
			prod*=pri[i];
		}
	}
	return ret;
}
int main()
{
	freopen("t4.in","r",stdin);
	//freopen("t4.out","w",stdout);
	scanf("%lld\n",&n);
	predo(1e6);
	sieve_g();
	nn=n;
	printf("%lld\n",(S(n,1)+1+mo)%mo);
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值