【XSY2721】求和 杜教筛

题目描述

  设(n=prod a_i^{p_i}),那么定义(f_d(n)=prod{(-1)^{p_i}[p_ileq d]})。特别的,(f_1(n)=mu(n))

  给你(n,k),求

[sum_{i=1}^nsum_{j=1}^nsum_{d=1}^kf_d(gcd(i,j)) ]

  (nleq {10}^{10},kleq 40)

题解

  先做一些简单的处理

[egin{align} ans&=sum_{i=1}^nsum_{j=1}^nsum_{d=1}^kf_d(gcd(i,j))\ &=sum_{i=1}^nsum_{d=1}^kf_d(i)(2(sum_{j=1}^{lfloorfrac{n}{i} floor}varphi(j))-1) end{align} ]

  后面(varphi)用杜教筛可以在(O(n^frac{2}{3}))内搞出来。

  设(lambda(n)=f_infty(n)=prod{(-1)}^{p_i})

  考虑容斥,有:

[f_d(n)=lambda(n)sum_{i^d|n}mu(i) ]

[egin{align} F_d(n)&=sum_{i=1}^nf_d(i)\ &=sum_{i=1}^nlambda(i)sum_{j^d|i}mu(j)\ &=sum_{i=1}^nmu(i)sum_{j=1}^{lfloorfrac{n}{i^{d+1}} floor}lambda(i^{d+1}j)\ &=sum_{i=1}^{lfloorsqrt[d+1]{n} floor}lambda^{d+1}(i)mu(i)Lambda(lfloorfrac{n}{i^{d+1}} floor) end{align} ]

  (nleq {10}^7)的部分预处理,其他的每次枚举。这部分每次枚举是(sqrt{n})的。总的是(O(n^frac{2}{3}))的。(和杜教筛的分析方法一样。)

[egin{align} sum_{j|i}lambda(j)&=[i是完全平方数]\ sum_{i=1}^nsum_{j|i}lambda(j)&=sqrt{n}\ sqrt{n}=sum_{i=1}^nsum_{j}[j|i]lambda(i)&=sum_{frac{i}{j}=1}^nsum_{j=1}^{lfloorfrac{n}{frac{i}{j}} floor}lambda(j) =sum_{i=1}^nLambda(lfloorfrac{n}{i} floor)\ Lambda(n)&=sqrt {n}-sum_{i=2}^nLambda(lfloorfrac{n}{i} floor) end{align} ]

  后面(Lambda(n))用杜教筛可以在(O(n^frac{2}{3}))内搞出来

  反正总的是(O(n^frac{2}{3}))的就对了。。。

  时间复杂度:(O(n^frac{2}{3}))

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int getsqrt(ll n)
{
	int l=1,r=1000000;
	while(l<r)
	{
		int mid=(l+r+1)>>1;
		if((ll)mid*mid>n)
			r=mid-1;
		else
			l=mid;
	}
	return l;
}
ll n;
ll _sqrt;
namespace hashset
{
	int getnum(ll x)
	{
		return n/x;
	}
}
using hashset::getnum;
int miu[10000010];
int phi[10000010];
int c[10000010];
int cs[10000010];
const int maxn=10000000;
int b[10000010];
int pri[1000010];
int cnt;
int d[10000010];
int e[10000010];
int f[10000010];
int k;
int vis[10000010];
int qp[10000010];
int qc[10000010];
void init()
{
	c[1]=phi[1]=miu[1]=f[1]=e[1]=1;
	d[1]=f[1]=0;
	int i,j;
	for(i=2;i<=maxn;i++)
	{
		if(!b[i])
		{
			miu[i]=-1;
			phi[i]=i-1;
			c[i]=-1;
			pri[++cnt]=i;
			d[i]=e[i]=1;
			f[i]=1;
		}
		for(j=1;j<=cnt&&i*pri[j]<=maxn;j++)
		{
			b[i*pri[j]]=1;
			c[i*pri[j]]=-c[i];
			f[i*pri[j]]=f[i]+1;
			if(i%pri[j]==0)
			{
				miu[i*pri[j]]=0;
				phi[i*pri[j]]=phi[i]*pri[j];
				d[i*pri[j]]=d[i]+1;
				e[i*pri[j]]=max(d[i*pri[j]],e[i]);
				break;
			}
			d[i*pri[j]]=1;
			e[i*pri[j]]=e[i];
			miu[i*pri[j]]=-miu[i];
			phi[i*pri[j]]=phi[i]*phi[pri[j]];
		}
	}
	for(i=1;i<=maxn;i++)
	{
		if(e[i]>k)
			f[i]=0;
		else
			f[i]=(f[i]&1?-1:1)*(k-e[i]+1);
		f[i]+=f[i-1];
//		miu[i]+=miu[i-1];
		phi[i]+=phi[i-1];
		cs[i]=cs[i-1]+c[i];
	}
}
int getphi(ll n)
{
	if(n<=maxn)
		return phi[n];
	int x=getnum(n);
	if(vis[x]&1)
		return qp[x];
	vis[x]|=1;
	ll i,j;
	int s=n*(n+1)>>1;
	for(i=2;i<=n;i=j+1)
	{
		j=n/(n/i);
		s-=(j-i+1)*getphi(n/i);
	}
	qp[x]=s;
	return s;
}
int getc(ll n)
{
	if(n<=maxn)
		return cs[n];
	int x=getnum(n);
	if(vis[x]&2)
		return qc[x];	
	vis[x]|=2;
	int s=getsqrt(n);
	ll i,j;
	for(i=2;i<=n;i=j+1)
	{
		j=n/(n/i);
		s-=(j-i+1)*getc(n/i);
	}
	qc[x]=s;
	return s;
}
ll pw[1000010];
int pw2[1000010];
int pw3[1000010];
int getfd(ll n)
{
	if(n<=maxn)
		return f[n];
	int i,j;
	for(i=1;(ll)i*i<=n;i++)
	{
		pw[i]=i;
		pw2[i]=pw3[i]=c[i];
	}
	int m=i-1;
	int s=0;
	for(j=1;j<=k;j++)
	{
		for(i=1;i<=m;i++)
		{
			pw[i]*=i;
			if(pw[i]>n)
				break;
			pw2[i]*=pw3[i];
			s+=miu[i]*pw2[i]*getc(n/pw[i]);
		}
	}
	return s;
}
int main()
{
#ifndef ONLINE_JUDGE
	freopen("b.in","r",stdin);
	freopen("b.out","w",stdout);
#endif
	scanf("%lld%d",&n,&k);
	_sqrt=getsqrt(n);
	init();
	int s=0;
	ll i,j;
	int now,last=0;
	int ans=0;
	for(i=1;i<=n;i=j+1)
	{
		j=n/(n/i);
		now=getfd(j);
		ans+=(now-last)*(2*getphi(n/i)-1);
		last=now;
	}
	ans&=(1<<30)-1;
	printf("%d
",ans);
	return 0;
}
原文地址:https://www.cnblogs.com/ywwyww/p/8513546.html