【51NOD1847】奇怪的数学题 min_25筛

题目描述

　　记\(sgcd(i,j)\)为\(i,j\)的次大公约数。

　　给你\(n\)，求

\[\sum_{i=1}^n\sum_{j=1}^n{sgcd(i,j)}^k \]

　　对\(2^{32}\)取模。

　　\(n\leq {10}^9,k\leq 50\)

题解

　　记\(f(n)\)为\(n\)的次大因数

　　显然\(sgcd(i,j)=f(gcd(i,j))\)

　　先推一波式子。

\[\begin{align} &\sum_{i=1}^n\sum_{j=1}^n{sgcd(i,j)}^k\\ =&\sum_{i=1}^n{f(i)}^k(2\sum_{j=1}^\frac{n}{i}\varphi(j)-1) \end{align} \]

　　后面那个\(\varphi(i)\)的前缀和可以杜教筛/min_25筛+数论分块解决，所以只用关心前面这部分。

　　观察min_25筛求质数的\(k\)次方的前缀和的过程，发现在求\(f_{n,j}=\sum_{i=2}^n[i\text{是质数或}i\text{的每个质因子都}>p_j]i^k\)的时候，减掉的那部分就是（最小质因子\(=p_j\)的数除以\(p_j\)后的值）的\(k\)次方和。所以直接把这些东西加起来就是我们要求的答案了。

　　还要加上质数的答案，即区间质数个数。

　　自然数幂求和要用第二类斯特林数那个做法。

\[\begin{align} S_m(n)&=\sum_{i=1}^ni^m\\ &=\sum_{i=1}^n\sum_{j=1}^m\begin{Bmatrix}m\\j\end{Bmatrix}i^\underline{j}\\ &=\sum_{j=1}^m\begin{Bmatrix}m\\j\end{Bmatrix}\sum_{i=1}^ni^\underline{j}\\ &=\sum_{j=1}^m\begin{Bmatrix}m\\j\end{Bmatrix}\frac{{(n+1)}^\underline{j+1}}{j+1} \end{align} \]

　　时间复杂度：\(O(\sqrt{n}k^2+\frac{n^\frac{3}{4}}{\log n})\)或\(O(\sqrt{n}k^2+\frac{n^\frac{3}{4}}{\log n}+n^\frac{2}{3})\)

代码

1

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<utility>
#include<cmath>
#include<functional>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
void sort(int &a,int &b)
{
	if(a>b)
		swap(a,b);
}
void open(const char *s)
{
#ifndef ONLINE_JUDGE
	char str[100];
	sprintf(str,"%s.in",s);
	freopen(str,"r",stdin);
	sprintf(str,"%s.out",s);
	freopen(str,"w",stdout);
#endif
}
int rd()
{
	int s=0,c,b=0;
	while(((c=getchar())<'0'||c>'9')&&c!='-');
	if(c=='-')
	{
		c=getchar();
		b=1;
	}
	do
	{
		s=s*10+c-'0';
	}
	while((c=getchar())>='0'&&c<='9');
	return b?-s:s;
}
void put(int x)
{
	if(!x)
	{
		putchar('0');
		return;
	}
	static int c[20];
	int t=0;
	while(x)
	{
		c[++t]=x%10;
		x/=10;
	}
	while(t)
		putchar(c[t--]+'0');
}
int upmin(int &a,int b)
{
	if(b<a)
	{
		a=b;
		return 1;
	}
	return 0;
}
int upmax(int &a,int b)
{
	if(b>a)
	{
		a=b;
		return 1;
	}
	return 0;
}
ll n;
int k;
int fp(int a,int b)
{
	int s=1;
	for(;b;b>>=1,a*=a)
		if(b&1)
			s*=a;
	return s;
}
int S[60][60];
int s[60];
int calc(ll x)
{
	s[0]=x;
	for(int i=1;i<=k;i++)
	{
		s[i]=1;
		ll v=(x+1)/(i+1)*(i+1);
		s[i]*=(x+1)/(i+1);
		for(ll j=x+1;j>=x-i+1;j--)
			s[i]*=(j==v?1:j);
		for(int j=0;j<i;j++)
			s[i]-=S[i][j]*s[j];
	}
	return s[k];
}
void init()
{
	S[0][0]=1;
	for(int i=1;i<=k;i++)
		for(int j=1;j<=i;j++)
			S[i][j]=S[i-1][j-1]-(i-1)*S[i-1][j];
}
namespace gao1
{
	const int m=100000;
	const int M=100010;
	int b[M],pri[M];
	int cnt;
	int f1[M],f2[M];
	int g1[M],g2[M];
	int ans1[M],ans2[M];
	void init()
	{
		for(int i=2;i<=m;i++)
		{
			if(!b[i])
				pri[++cnt]=i;
			for(int j=1;j<=cnt&&i*pri[j]<=m;j++)
			{
				b[i*pri[j]]=1;
				if(i%pri[j]==0)
					break;
			}
		}
	}
	void gao()
	{
		init();
		for(int i=1;i<=m;i++)
		{
			f1[i]=calc(i)-1;
			g1[i]=i-1;
		}
		for(int i=1;n/i>m;i++)
		{
			f2[i]=calc(n/i)-1;
			g2[i]=n/i-1;
		}
		for(int i=1;i<=cnt;i++)
		{
			int j;
			int v=fp(pri[i],k);
			for(j=1;n/pri[i]/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
			{
				ans2[j]+=f2[pri[i]*j]-f1[pri[i]-1];
				f2[j]-=v*(f2[pri[i]*j]-f1[pri[i]-1]);
				g2[j]-=g2[pri[i]*j]-g1[pri[i]-1];
			}
			for(;n/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
			{
				ans2[j]+=f1[n/pri[i]/j]-f1[pri[i]-1];
				f2[j]-=v*(f1[n/pri[i]/j]-f1[pri[i]-1]);
				g2[j]-=g1[n/pri[i]/j]-g1[pri[i]-1];
			}
			for(j=m;j>=2&&j>=(ll)pri[i]*pri[i];j--)
			{
				ans1[j]+=f1[j/pri[i]]-f1[pri[i]-1];
				f1[j]-=v*(f1[j/pri[i]]-f1[pri[i]-1]);
				g1[j]-=g1[j/pri[i]]-g1[pri[i]-1];
			}
		}
		for(int i=1;i<=m;i++)
			ans1[i]+=g1[i];
		for(int i=1;n/i>m;i++)
			ans2[i]+=g2[i];
	}
	int query(ll x)
	{
		return x<=m?ans1[x]:ans2[n/x];
	}
}
namespace gao2
{
	const int m=1000000;
	const int M=1000010;
	int pri[M],b[M],cnt;
	int phi[M],s[M];
	int b2[M],s2[M];
	void init()
	{
		phi[1]=1;
		for(int i=2;i<=m;i++)
		{
			if(!b[i])
			{
				pri[++cnt]=i;
				phi[i]=i-1;
			}
			for(int j=1;j<=cnt&&i*pri[j]<=m;j++)
			{
				b[i*pri[j]]=1;
				if(i%pri[j]==0)
				{
					phi[i*pri[j]]=phi[i]*pri[j];
					break;
				}
				phi[i*pri[j]]=phi[i]*phi[pri[j]];
			}
		}
		for(int i=1;i<=m;i++)
			s[i]=s[i-1]+phi[i];
	}
	void gao()
	{
		init();
	}
	int query(ll x)
	{
		if(x<=m)
			return s[x];
		if(b2[n/x])
			return s2[n/x];
		b2[n/x]=1;
		int &res=s2[n/x];
		res=x*(x+1)/2;
		for(ll i=2,j;i<=x;i=j+1)
		{
			j=x/(x/i);
			res-=query(x/i)*(j-i+1);
		}
		return res;
	}
}
int main()
{
	open("51nod1847");
	scanf("%lld%d",&n,&k);
	init();
	gao1::gao();
	gao2::gao();
	int ans=0;
	for(ll i=2,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		ans+=(gao1::query(j)-gao1::query(i-1))*(2*gao2::query(n/i)-1);
	}
	printf("%u\n",(unsigned)ans);
	return 0;
}

2

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<utility>
#include<cmath>
#include<functional>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
void sort(int &a,int &b)
{
	if(a>b)
		swap(a,b);
}
void open(const char *s)
{
#ifndef ONLINE_JUDGE
	char str[100];
	sprintf(str,"%s.in",s);
	freopen(str,"r",stdin);
	sprintf(str,"%s.out",s);
	freopen(str,"w",stdout);
#endif
}
int rd()
{
	int s=0,c,b=0;
	while(((c=getchar())<'0'||c>'9')&&c!='-');
	if(c=='-')
	{
		c=getchar();
		b=1;
	}
	do
	{
		s=s*10+c-'0';
	}
	while((c=getchar())>='0'&&c<='9');
	return b?-s:s;
}
void put(int x)
{
	if(!x)
	{
		putchar('0');
		return;
	}
	static int c[20];
	int t=0;
	while(x)
	{
		c[++t]=x%10;
		x/=10;
	}
	while(t)
		putchar(c[t--]+'0');
}
int upmin(int &a,int b)
{
	if(b<a)
	{
		a=b;
		return 1;
	}
	return 0;
}
int upmax(int &a,int b)
{
	if(b>a)
	{
		a=b;
		return 1;
	}
	return 0;
}
ll n;
int k;
int fp(int a,int b)
{
	int s=1;
	for(;b;b>>=1,a*=a)
		if(b&1)
			s*=a;
	return s;
}
int S[60][60];
int s[60];
int calc(ll x)
{
	s[0]=x;
	for(int i=1;i<=k;i++)
	{
		s[i]=1;
		ll v=(x+1)/(i+1)*(i+1);
		s[i]*=(x+1)/(i+1);
		for(ll j=x+1;j>=x-i+1;j--)
			s[i]*=(j==v?1:j);
		for(int j=0;j<i;j++)
			s[i]-=S[i][j]*s[j];
	}
	return s[k];
}
void init()
{
	S[0][0]=1;
	for(int i=1;i<=k;i++)
		for(int j=1;j<=i;j++)
			S[i][j]=(S[i-1][j-1]-(i-1)*S[i-1][j]);
}
namespace gao1
{
	const int m=100000;
	const int M=100010;
	int b[M],pri[M];
	int cnt;
	int f1[M],f2[M];
	int ans1[M],ans2[M];
//	vector<int> s1[M],s2[M];
//	vector<short> b1[M],b2[M];
	void init()
	{
		for(int i=2;i<=m;i++)
		{
			if(!b[i])
				pri[++cnt]=i;
			for(int j=1;j<=cnt&&i*pri[j]<=m;j++)
			{
				b[i*pri[j]]=1;
				if(i%pri[j]==0)
					break;
			}
		}
	}
	void gao()
	{
		init();
		for(int i=1;i<=m;i++)
			f1[i]=calc(i)-1;
		for(int i=1;n/i>m;i++)
			f2[i]=calc(n/i)-1;
		for(int i=1;i<=cnt;i++)
		{
			int j;
			int v=fp(pri[i],k);
			for(j=1;n/pri[i]/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
			{
				int x=f2[pri[i]*j]-f1[pri[i]-1];
				ans2[j]+=x;
				f2[j]-=v*x;
			}
			for(;n/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
			{
				int x=f1[n/pri[i]/j]-f1[pri[i]-1];;
				ans2[j]+=x;
				f2[j]-=v*x;
			}
			for(j=m;j>=(ll)pri[i]*pri[i];j--)
			{
				int x=f1[j/pri[i]]-f1[pri[i]-1];
				ans1[j]+=x;
				f1[j]-=v*x;
			}
		}
	}
	int query(ll x)
	{
		return x<=m?ans1[x]:ans2[n/x];
	}
}
namespace gao2
{
	const int m=100000;
	const int M=100010;
	int pri[M],b[M],cnt;
	int f1[M],f2[M],g1[M],g2[M];
	int s1[M],s2[M];
	int b1[M],b2[M];
	void init()
	{
		for(int i=2;i<=m;i++)
		{
			if(!b[i])
				pri[++cnt]=i;
			for(int j=1;j<=cnt&&i*pri[j]<=m;j++)
			{
				b[i*pri[j]]=1;
				if(i%pri[j]==0)
					break;
			}
		}
		pri[cnt+1]=m+1;
		pri[0]=1;
	}
	void gao()
	{
		init();
		for(int i=1;i<=m;i++)
		{
			f1[i]=(ll)(i+2)*(i-1)/2;
			g1[i]=i-1;
		}
		for(int i=1;n/i>m;i++)
		{
			f2[i]=(n/i&1?(n/i-1)/2*(n/i+2):(n/i+2)/2*(n/i-1));
			g2[i]=n/i-1;
		}
		for(int i=1;i<=cnt;i++)
		{
			int j;
			int x=f1[pri[i]-1];
			int y=g1[pri[i]-1];
			for(j=1;n/pri[i]/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
			{
				f2[j]-=pri[i]*(f2[pri[i]*j]-x);
				g2[j]-=g2[pri[i]*j]-y;
			}
			for(;n/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
			{
				f2[j]-=pri[i]*(f1[n/pri[i]/j]-x);
				g2[j]-=g1[n/pri[i]/j]-y;
			}
			for(j=m;j>=2&&j>=(ll)pri[i]*pri[i];j--)
			{
				f1[j]-=pri[i]*(f1[j/pri[i]]-x);
				g1[j]-=g1[j/pri[i]]-y;
			}
		}
		for(int i=1;i<=m;i++)
			gao1::ans1[i]+=g1[i];
		for(int i=1;n/i>m;i++)
			gao1::ans2[i]+=g2[i];
		for(int i=1;i<=m;i++)
			f1[i]-=g1[i];
		for(int i=1;n/i>m;i++)
			f2[i]-=g2[i];
		for(int i=1;i<=m;i++)
			s1[i]=f1[i];
		for(int i=1;n/i>m;i++)
			s2[i]=f2[i];
		for(int i=cnt;i>=1;i--)
		{
			int x=pri[i];
			for(int j=1;n/j>m&&n/j>=(ll)x*x;j++)
			{
				s2[j]-=x-1;
				ll x2=n/j/x;
				int s=x-1;
				for(;x2;x2/=x,s*=x)
					s2[j]+=((x2<=m?s1[x2]-f1[min((int)x2,pri[i])]:s2[n/x2]-f1[pri[i]])+1)*s;
			}
			for(int j=m;j>=(ll)x*x;j--)
			{
				s1[j]-=x-1;
				int x2=j/x;
				int s=x-1;
				for(;x2;x2/=x,s*=x)
					s1[j]+=(s1[x2]-f1[min(x2,pri[i])]+1)*s;
			}
		}
	}
	int query(ll x)
	{
		return (x<=m?s1[x]:s2[n/x])+1;
	}
}
int main()
{
	open("51nod1847");
	scanf("%lld%d",&n,&k);
	init();
	gao1::gao();
	gao2::gao();
	int ans=0;
	for(ll i=2,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		ans+=(gao1::query(j)-gao1::query(i-1))*(2*gao2::query(n/i)-1);
	}
	printf("%u\n",ans);
	return 0;
}