Min_25's Sieve

2019-05-03 00:00:00 +0000 | DOlaBMOon

tags: math.

我们有一个函数 , 它满足如下性质:

  1. ;
  2. ;
  3. .

我们想要计算:

Step1.

我们先想要计算

既然, 是一个简单多项式, 那么我们就可以直接对所有的项分别计算和.

设 为 . 为 的最小素因子, 为第 个素数, 假装 , 是将 带入素数的多项式, 就是假装所有数都是素数的 .

但, 当然我们现在计算的是某一项的和, 所以 .

这是我们想要的.

比如, 如果我们要求 :

for(int tp=1;pri[tp]*pri[tp]<=n;++tp)
	for(int i=1;i<=cnt&&pri[tp]*pri[tp]<=w[i];++i)
		g[i]-=g[(t=w[i]/pri[tp])<lim?id[t]:id2[n/t]]-tp+1;

Step2.

有递归和不递归的写法, 递归的比较快, 但不递归的方法在时间复杂度不变的情况下, 还求出了很多别的东西, 有时候有用 (如 ).

递归方法

设 .

这是我们要求的.

然后就可以算啦.

SPOJ DIVCNTK

#include<iostream>

using namespace std;

#define ULL unsigned long long

const int M=3e5+10;

bool notp[M];
ULL pri[M];
int totp;
void Prework()
{
	pri[0]=1;
	for(int i=2;i<M;++i)
	{
		if(!notp[i])
			pri[++totp]=i;
		for(int j=1,k;j<=totp&&(k=pri[j]*i)<M;++j)
		{
			notp[k]=true;
			if(i%pri[j]==0)
				break;
		}
	}
}

int cnt,tp;
ULL w[M],n,c,g[M],id[M],id2[M],lim;

ULL S(ULL x,int y)
{
	if(x<=1||pri[y]>=n)
		return 0;
	ULL res=(x<lim?g[id[x]]:g[id2[n/x]])-(c+1)*y,p;
	for(int j=y+1;j<=totp&&(p=pri[j])*pri[j]<=x;++j)
		for(int k=1;p*pri[j]<=x;p*=pri[j],++k)
			res+=S(x/p,j)*(c*k+1)+c*(k+1)+1;
	return res;
}

void Work()
{
	scanf("%llu%llu",&n,&c);
	for(lim=1;lim*lim<=n;++lim);
	cnt=0;
	for(ULL i=1,k;i<=n;i=n/k+1)
	{
		w[++cnt]=k=n/i;
		g[cnt]=k-1;
		(k<lim?id[k]:id2[n/k])=cnt;
	}
	ULL t;
	for(tp=1;pri[tp]*pri[tp]<=n;++tp)
		for(int i=1;i<=cnt&&pri[tp]*pri[tp]<=w[i];++i)
			g[i]-=g[(t=w[i]/pri[tp])<lim?id[t]:id2[n/t]]-tp+1;
	for(int i=1;i<=cnt;++i)
		g[i]*=c+1;
	printf("%llu\n",S(n,0)+1);
}

signed main()
{
	Prework();
	int T;
	for(scanf("%d",&T);T--;Work());
	return 0;
}

不递归方法

这个看起来更正统一些:

设 . 这里的 是真的 哦.

这是我们要求的.

SPOJ DIVCNTK

#include<iostream>

using namespace std;

#define ULL unsigned long long

const int M=3e5+10;

bool notp[M];
ULL pri[M];
int totp;
void Prework()
{
	pri[0]=1;
	for(int i=2;i<M;++i)
	{
		if(!notp[i])
			pri[++totp]=i;
		for(int j=1,k;j<=totp&&(k=pri[j]*i)<M;++j)
		{
			notp[k]=true;
			if(i%pri[j]==0)
				break;
		}
	}
}

int cnt,tp;
ULL w[M],n,c,g[M],s[M],id[M],id2[M],lim;

void Work()
{
	scanf("%llu%llu",&n,&c);
	for(lim=1;lim*lim<=n;++lim);
	cnt=0;
	for(ULL i=1,k;i<=n;i=n/k+1)
	{
		w[++cnt]=k=n/i;
		g[cnt]=k-1;
		(k<lim?id[k]:id2[n/k])=cnt;
	}
	ULL t,d,p;
	for(tp=1;pri[tp]*pri[tp]<=n;++tp)
		for(int i=1;i<=cnt&&pri[tp]*pri[tp]<=w[i];++i)
			g[i]-=g[(t=w[i]/pri[tp])<lim?id[t]:id2[n/t]]-tp+1;
	for(int i=1;i<=cnt;++i)
		s[i]=(g[i]*=c+1);
	for(--tp;tp;--tp)
	{
		d=tp*(c+1);
		for(int i=1;i<=cnt&&(p=pri[tp])*pri[tp]<=w[i];++i)
			for(int k=1;p*pri[tp]<=w[i];p*=pri[tp],++k)
				s[i]+=(k*c+1)*(s[(t=w[i]/p)<lim?id[t]:id2[n/t]]-d)+(k+1)*c+1;
	}
	printf("%llu\n",s[1]+1);
}

signed main()
{
	Prework();
	int T;
	for(scanf("%d",&T);T--;Work());
	return 0;
}