【UOJ448】【集训队作业2018】人类的本质 min_25筛

题目大意

　　给你 \(n,m\)，求

\[\sum_{i=1}^n\sum_{x_1,x_2,\ldots,x_m=1}^i\operatorname{lcm}(\gcd(i,x_1),\gcd(i,x_2),\ldots,\gcd(i,x_m)) \]

　　对 \({10}^9+7\) 取模。

　　\(nm\leq {10}^9\)

题解

　　先推一下式子：

\[ans=\sum_{i=1}^n\sum_{x_1,x_2,\ldots,x_m=1}^i\operatorname{lcm}(\gcd(i,x_1),\gcd(i,x_2),\ldots,\gcd(i,x_m))\\ =\sum_{i=1}^n\sum_{x_1,x_2,\ldots,x_m=1}^i\frac{i}{\gcd(\frac{i}{\gcd(i,x_1)},\frac{i}{\gcd(i,x_2)},\ldots,\frac{i}{\gcd(i,x_m)})}\\ \]

　　然后枚举 \(y_j=\frac{i}{\gcd(i,x_j)}\)。对于一个 \(y_j\)，有多少个满足要求的 \(x_j\) 呢？

\[\frac{i}{\gcd(i,x_j)}=y_j\\ \gcd(i,x_j)=\frac{i}{y_j}\\ \gcd(y_j,\frac{x_j}{y_j})=1\\ \]

　　所以有 \(\varphi(y_j)\) 个。

\[ans=\sum_{i=1}^n\sum_{y_1,y_2,\ldots,y_m\mid i}\varphi(y_1)\varphi(y_2)\cdots\varphi(y_m)\frac{i}{\gcd(y_1,y_2,\ldots,y_m)} \]

　　记

\[f(n)=\sum_{y_1,y_2,\ldots,y_m\mid n}\varphi(y_1)\varphi(y_2)\cdots\varphi(y_m)\frac{n}{\gcd(y_1,y_2,\ldots,y_m)} \]

　　容易发现 \(f(n)\) 是积性函数。

\[f(p^k)=\sum_{i=0}^kp^{k-i}\sum_{j=i}^k\mu(\frac{j}{i})g(j)\\ =\sum_{j=0}^kg(j)\sum_{i=0}^jp^{k-i}\mu(\frac{j}{i})\\ g(j)=\begin{cases} p^{km}&,j=1\\ (p^k-p^{j-1})^m&,2\leq j\leq k \end{cases}\\ \]

\[f(p^k)=\sum_{x_1,x_2,\ldots,x_m=1}^{p^k}lcm(\gcd(p^k,x_1),\gcd(p^k,x_2),\ldots,\gcd(p^k,x_m))\\ =\sum_{x_1,x_2,\ldots,x_m=0}^{k}p^{\max(x_1,x_2,\ldots,x_m)}\prod_{j=1}^m\varphi(p^{k-x_j})\\ =\sum_{i=0}^kp^i(g(i)-g(i-1))\\ g(i)=\sum_{x_1,x_2,\dots,x_m=0}^i\prod_{j=1}^m\varphi(p^{k-x_j})\\ =\prod_{j=1}^m(\sum_{l=0}^i\varphi(p^{k-l}))\\= \begin{cases} p^{km}&,i=k\\ (p^k-p^{k-1-i})^m&,0\leq i<k \end{cases}\\ f(p^{k+1})=p^mf(p^k)-p^k(p^{(k+1)m}-(p^{k+1}-1)^m)+p^{k+1}(p^{(k+1)m}-(p^{k+1}-1)^m)\\ =p^mf(p^k)+p^k(p-1)(p^{(k+1)m}-(p^{k+1}-1)^m) \]

　　当 \(n\leq {10}^7\) 的时候可以线性筛算出 \(f(1)\sim f(n)\)。

　　时间复杂度：\(O(n(1+\frac{\log m}{\log n}))\)

　　当 \(n\) 大的时候：

\[f(p)=g(0)+pg(1)-pg(0)\\ =(1-p)(p-1)^m+p^{m+1}\\ =p^{m+1}-(p-1)^{m+1}\\ =-\sum_{i=0}^m{(-1)}^{m+1-i}p^i\binom{m+1}{i}\\ =\sum_{i=0}^m{(-1)}^{m-i}p^i\binom{m+1}{i} \]

　　就可以上 min_25 筛了。

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

代码

　　这个代码有些东西没有预处理，时间复杂度可能比上面写的高一些。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<functional>
#include<cmath>
#include<vector>
#include<assert.h>
//using namespace std;
using std::min;
using std::max;
using std::swap;
using std::sort;
using std::reverse;
using std::random_shuffle;
using std::lower_bound;
using std::upper_bound;
using std::unique;
using std::vector;
typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef std::pair<int,int> pii;
typedef std::pair<ll,ll> pll;
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
}
void open2(const char *s){
#ifdef DEBUG
	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;}
const ll p=1000000007;
ll fp(ll a,ll b)
{
	ll s=1;
	for(;b;b>>=1,a=a*a%p)
		if(b&1)
			s=s*a%p;
	return s;
}
int n,m;
namespace gao1
{
	const int N=10000010;
	int f[N];
	int pw[N];
	bool b[N];
	int pri[N];
	int c[N];
	int cnt;
	void gao()
	{
		f[1]=1;
		pw[1]=1;
		for(int i=2;i<=n;i++)
		{
			if(!b[i])
			{
				pri[++cnt]=i;
				pw[i]=fp(i,m);
			}
			for(int j=1;j<=cnt&&i*pri[j]<=n;j++)
			{
				b[i*pri[j]]=1;
				pw[i*pri[j]]=(ll)pw[i]*pw[pri[j]]%p;
				if(i%pri[j]==0)
					break;
			}
		}
		for(int i=2;i<=n;i++)
		{
			if(!b[i])
			{
				c[i]=i;
				f[i]=fp(i,m+1)-fp(i-1,m+1);
			}
			for(int j=1;j<=cnt&&pri[j]<=i&&i*pri[j]<=n;j++)
			{
				int v=i*pri[j];
				if(i%pri[j]==0)
				{
					c[v]=c[i]*pri[j];
					if(c[v]==v)
						f[v]=((ll)pw[pri[j]]*f[i]+(ll)i*(pri[j]-1)*(pw[v]-pw[v-1]))%p;
					else
						f[v]=(ll)f[c[v]]*f[v/c[v]]%p;
					break;
				}
				c[v]=pri[j];
				f[v]=(ll)f[i]*f[pri[j]]%p;
			}
		}
		ll ans=0;
		for(int i=1;i<=n;i++)
			ans+=f[i];
		ans=(ans%p+p)%p;
		printf("%lld\n",ans);
	}
}
namespace gao2
{
	const int N=100010;
	ll inv[N],fac[N],ifac[N];
	ll binom(int x,int y)
	{
		return x>=y&&y>=0?fac[x]*ifac[y]%p*ifac[x-y]%p:0;
	}
	int _;
	int b[N],pri[N],cnt;
	void init()
	{
		inv[1]=fac[0]=fac[1]=ifac[0]=ifac[1]=1;
		for(int i=2;i<=100000;i++)
		{
			inv[i]=-p/i*inv[p%i]%p;
			fac[i]=fac[i-1]*i%p;
			ifac[i]=ifac[i-1]*inv[i]%p;
		}
		for(int i=2;i<=_;i++)
		{
			if(!b[i])
				pri[++cnt]=i;
			for(int j=1;j<=cnt&&i*pri[j]<=_;j++)
			{
				b[i*pri[j]]=1;
				if(i%pri[j]==0)
					break;
			}
		}
	}
	ll pw[N],s[N];
	void init(int x)
	{
		pw[1]=1;
		for(int i=2;i<=x+2;i++)
		{
			if(!b[i])
				pw[i]=fp(i,x);
			for(int j=1;j<=cnt&&i*pri[j]<=n&&pri[j]<=i;j++)
			{
				pw[i*pri[j]]=pw[i]*pw[pri[j]]%p;
				if(i%pri[j]==0)
					break;
			}
		}
		for(int i=1;i<=x+2;i++)
			s[i]=(s[i-1]+pw[i])%p;
	}
	ll pre[N],suf[N];
	ll calc(int x,int y)
	{
		if(x<=y+2)
			return s[x];
		pre[0]=1;
		for(int i=1;i<=y+2;i++)
			pre[i]=pre[i-1]*(x-i)%p;
		suf[y+3]=1;
		for(int i=y+2;i>=1;i--)
			suf[i]=suf[i+1]*(x-i)%p;
		ll res=0;
		for(int i=1;i<=y+2;i++)
			res=(res+s[i]*pre[i-1]%p*suf[i+1]%p*ifac[i-1]%p*ifac[y+2-i]%p*((y+2-i)&1?-1:1))%p;
		return res;
	}
	ll s1[N],s2[N],ans1[N],ans2[N];
	void gao(int x)
	{
		init(x);
		calc(3,x);
		for(int i=1;i<=_;i++)
			s1[i]=(fp(i,x)+s1[i-1])%p;
		for(int i=1;n/i>_;i++)
			s2[i]=calc(n/i,x);
		for(int i=1;i<=cnt;i++)
		{
			ll v=fp(pri[i],x);
			int j;
			for(j=1;n/j/pri[i]>_&&n/j>=(ll)pri[i]*pri[i];j++)
				s2[j]=(s2[j]-v*(s2[j*pri[i]]-s1[pri[i]-1]))%p;
			for(;n/j>_&&n/j>=(ll)pri[i]*pri[i];j++)
				s2[j]=(s2[j]-v*(s1[n/j/pri[i]]-s1[pri[i]-1]))%p;
			for(j=_;j>=(ll)pri[i]*pri[i];j--)
				s1[j]=(s1[j]-v*(s1[j/pri[i]]-s1[pri[i]-1]))%p;
		}
	}
	ll g(int z,int x,int y)
	{
		if(z<0)
			return 0;
		return z==y?fp(x,(ll)y*m):fp(fp(x,y)-fp(x,y-z-1),m);
	}
	ll get(int x,int y)
	{
		ll res=0;
		for(int i=0;i<=y;i++)
			res=(res+fp(x,i)*(g(i,x,y)-g(i-1,x,y)))%p;
		return res;
	}
	ll solve(ll x,int y)
	{
		if(x<=1||x<pri[y])
			return 0;
		if(y>cnt)
			return (x<=_?ans1[x]:ans2[n/x])-ans1[_];
		ll res=(x<=_?ans1[x]:ans2[n/x])-ans1[pri[y]-1];
		for(int i=y;i<=cnt&&(ll)pri[i]*pri[i]<=x;i++)
		{
			ll x1=x/pri[i];
			for(int j=1;x1>=pri[i];j++,x1/=pri[i])
				res=(res+solve(x1,i+1)*get(pri[i],j)+get(pri[i],j+1))%p;
		}
		return res;
	}
	void gao()
	{
		_=sqrt(n)+0.5;
		init();
		for(int i=0;i<=m;i++)
		{
			gao(i);
			ll v=binom(m+1,i)*((m-i)&1?-1:1);
			for(int i=1;i<=_;i++)
				ans1[i]=(ans1[i]+v*s1[i])%p;
			for(int i=1;n/i>_;i++)
				ans2[i]=(ans2[i]+v*s2[i])%p;
		}
		ll ans=solve(n,1);
		ans++;
		ans=(ans%p+p)%p;
		printf("%lld\n",ans);
	}
}
int main()
{
	open("c");
	scanf("%d%d",&n,&m);
	if(n<=1e7)
		gao1::gao();
	else
		gao2::gao();
	return 0;
}