题目描述

　　三倍经验题。

　　给你\(n,m\)，求

　　\(n\leq {10}^9,1\leq m\leq 500000\)

题解

　　当\(m=1\)时\(ans=\frac{n(n+1)}{2}\)

　　剩下的部分这篇博客有讲YWW's Blog

　　时间复杂度：\(O(m+\log n)\)

代码

#include<cstdio>

#include<cstring>

#include<algorithm>

using namespace std;

typedef long long ll;

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 pri[100010];

int b[1000010];

int cnt;

ll s[1000010];

ll fac[1000010];

ll ifac[1000010];

ll inv[1000010];

ll f1[1000010];

ll f2[1000010];

ll f[1000010];

ll pre[1000010];

ll suf[1000010];

ll getc(int x,int y)

{

	return fac[x]*ifac[y]%p*ifac[x-y]%p;

}

int main()

{

#ifndef ONLINE_JUDGE

	freopen("bzoj4126.in","r",stdin);

	freopen("bzoj4126.out","w",stdout);

#endif

	int n,m;

	scanf("%d%d",&n,&m);

	if(m==1)

	{

		printf("%lld\n",ll(n)*(n+1)/2%p);

		return 0;

	}

	n++;

	fac[0]=fac[1]=ifac[0]=ifac[1]=inv[1]=1;

	for(int i=2;i<=m+2;i++)

	{

		inv[i]=-p/i*inv[p%i]%p;

		ifac[i]=ifac[i-1]*inv[i]%p;

		fac[i]=fac[i-1]*i%p;

	}

	s[0]=0;

	s[1]=1;

	for(int i=2;i<=m+2;i++)

	{

		if(!b[i])

		{

			s[i]=fp(i,m);

			pri[++cnt]=i;

		}

		for(int j=1;j<=cnt&&i*pri[j]<=m+2;j++)

		{

			b[i*pri[j]]=1;

			s[i*pri[j]]=s[i]*s[pri[j]]%p;

			if(i%pri[j]==0)

				break;

		}

	}

	f1[0]=1;

	f2[0]=0;

	ll invm=fp(m,p-2);

	for(int i=1;i<=m+1;i++)

	{

		f1[i]=f1[i-1]*invm%p;

		f2[i]=(f2[i-1]+s[i-1])*invm%p;

	}

	ll v1=0,v2=0;

	for(int i=0;i<=m+1;i++)

	{

		v1=(v1+((m+1-i)&1?-1:1)*getc(m+1,i)*f1[i])%p;

		v2=(v2+((m+1-i)&1?-1:1)*getc(m+1,i)*f2[i])%p;

	}

	f[0]=-v2*fp(v1,p-2)%p;

	for(int i=1;i<=m+1;i++)

		f[i]=(f1[i]*f[0]+f2[i])%p;

	if(n<=m+1)

	{

		ll ans=fp(m,n)*f[n]-f[0];

		ans=(ans%p+p)%p;

		printf("%lld\n",ans);

		return 0;

	}

	for(int i=0;i<=m;i++)

	{

		pre[i]=n-i;

		if(i)

			pre[i]=pre[i-1]*pre[i]%p;

	}

	for(int i=m;i>=0;i--)

	{

		suf[i]=n-i;

		if(i!=m)

			suf[i]=suf[i+1]*suf[i]%p;

	}

	ll ans=0;

	for(int i=0;i<=m;i++)

	{

		ll v=1;

		if(i)

			v=v*pre[i-1]%p;

		if(i!=m)

			v=v*suf[i+1]%p;

		ans=(ans+f[i]*v%p*ifac[i]%p*ifac[m-i]%p*((m-i)&1?-1:1))%p;

	}

	ans=fp(m,n)*ans-f[0];

	ans=(ans%p+p)%p;

	printf("%lld\n",ans);

	return 0;

}