codeforces 623 E FFT

最新推荐文章于 2022-08-11 23:17:14 发布

stony_oi

最新推荐文章于 2022-08-11 23:17:14 发布

阅读量316

点赞数

分类专栏： FFT

本文链接：https://blog.csdn.net/stony_oi/article/details/71170834

版权

FFT 专栏收录该内容

3 篇文章 0 订阅

订阅专栏

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
using namespace std;
typedef long long LL;
#define R register
const int N=66005,P=1000000007,bt=16,mas=65535;
const long double _2Pi=2*acos(-1.0);
struct na {long double x,y;};
inline na operator + (const na &a,const na &b) {return (na){a.x+b.x,a.y+b.y};}
inline na operator - (const na &a,const na &b) {return (na){a.x-b.x,a.y-b.y};}
inline na operator * (const na &a,const na &b) {return (na){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
inline na conj (const na &a) {return (na){a.x,-a.y};}
LL K;int n,m,X,rev[N],fac[N],fr[N],f[N],g[N],ans=0;
na ha,w[2][N/2],A[N],B[N],C[N],D[N];

inline int pw(int a,int k){
	int rec=1;
	for (;k;a=(LL)a*a%P,k>>=1) (k&1)?rec=(LL)rec*a%P:0;
	return rec;
}

void FFT(na A[],int r){
	for (int i=0;i<m;i++) if (i<rev[i]) swap(A[i],A[rev[i]]);
	for (R int i=1,l=m/2;i<m;i*=2,l>>=1)
		for (R int j=0;j<m;j+=2*i)
			for (R int k=0,d=0;k<i;k++,d+=l)
				ha=A[i+j+k]*w[r][d],A[i+j+k]=A[j+k]-ha,A[j+k]=A[j+k]+ha;
	if (r) for (int i=0;i<m;i++) A[i].x/=m,A[i].y/=m;
}

void conv(int f[],int g[],int h[],int len)
{
	for (int i=0;i<=n;i++)
	{
		int a=(LL)f[i]*pw(2,(int)((LL)i*len%(P-1)))%P,b=g[i];
		A[i]=(na){a>>bt,a&mas};
		B[i]=(na){b>>bt,b&mas};//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	}
	for (int i=n+1;i<=m;i++) A[i]=B[i]=(na){0,0};
	FFT(A,0),FFT(B,0);
	A[m]=A[0];B[m]=B[0];
	for (int i=0;i<m;i++)
	{
		na a1=(A[i]+conj(A[m-i]))*(na){0.5,0},a2=(A[i]-conj(A[m-i]))*(na){0,-0.5};
		na b1=(B[i]+conj(B[m-i]))*(na){0.5,0},b2=(B[i]-conj(B[m-i]))*(na){0,-0.5};
		C[i]=a1*b1+(na){0,1}*a2*b2;
		D[i]=a1*b2+a2*b1;
	}
	FFT(C,1),FFT(D,1);
	for (int i=0;i<=n;i++)
	{
		int a=(LL)round(C[i].x)%P,b=(LL)round(C[i].y)%P,c=(LL)round(D[i].x)%P;
		h[i]=(((LL)a<<2*bt)%P+((LL)c<<bt)%P+b)%P;
	}
}

int main()
{
//	freopen("1.in","r",stdin);
	cin>>K>>n;
	if (K>n) return puts("0"),0;
	for (m=2,X=1;m<=2*n;m*=2,X++);
	for (int i=1;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<X-1);
	for (int i=0;i<m/2;i++){
		long double k=_2Pi*i/m;
		w[0][i]=(na){cos(k),sin(k)},w[1][i]=(na){cos(-k),sin(-k)};
	}
	fac[0]=1;fr[0]=1;
	for (int i=1;i<=m;i++) fac[i]=(LL)fac[i-1]*i%P,fr[i]=pw(fac[i],P-2);
	for (int i=1;i<=n;i++) f[i]=fr[i];
	g[0]=1;
	for (int x=1;K;K>>=1,x*=2)
	{
		if (K&1)
			conv(g,f,g,x);
		conv(f,f,f,x);
	}
	for (int i=0;i<=n;i++) ans=(ans+(LL)g[i]*fac[n]%P*fr[n-i]%P)%P;
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	cout<<ans<<endl;
	return 0;
}