这个dp比较显然吧.
拆系数FFT的精度堪忧啊…
是时候改一下FFT的版了
1.所有的单位根都用欧拉公式求.
2.取整用round和llround,都是cmath中的.
3.拆系数可以增大精度.
多项式快速幂:
因为要mod
如果是NTT模数那最好,如果不是,那就拆系数FFT吧,反正也没有人想打CRTNTT
注意清零…和mod
AC Code:
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#define maxn 70005
#define M1 31623
#define M2 14122
#define mod 1000000007
#define ld long double
#define LL long long
using namespace std;
const double PI = 3.1415926535897932384626433832795;
LL n,k,lgn,len;
struct cplx
{
ld r,i;
cplx(ld r=0,ld i=0):r(r),i(i){}
cplx operator +(const cplx &B)const{ return cplx(r+B.r,i+B.i); }
cplx operator -(const cplx &B)const{ return cplx(r-B.r,i-B.i); }
cplx operator *(const cplx &B)const{ return cplx(r*B.r-i*B.i,i*B.r+r*B.i); }
cplx conj(){ return cplx(r,-i); }
}w[16][maxn]={1};
int r[maxn];
inline void FFT(cplx *A,int tp)
{
for(int i=0;i<len;i++) if(i < r[i]) swap(A[i],A[r[i]]);
for(int Len=2,k=0;Len<=len;Len<<=1)
{
int l = Len >> 1;
for(int st=0;st<len;st+=Len) for(int j=0;j<l;j++)
{
cplx tmp = (tp==1 ? w[k][j] : w[k][j].conj()) * A[st + j + l];
A[st + j + l] = A[st + j] - tmp;
A[st + j] = A[st + j] + tmp;
}
k++;
}
if(tp == -1) for(int i=0;i<len;i++) A[i].r/=len;
}
inline int Pow(int base,int k)
{
int ret = 1;
for(;k;k>>=1,base=1ll*base*base%mod) if(k&1) ret=1ll*ret*base%mod;
return ret;
}
inline LL si(cplx a){ return (LL)(a.r+0.5) %mod; }
void Move(int A[maxn],int B[maxn],int ret[maxn])
{
static cplx sta[2][2][maxn];
for(int i=0;i<len;i++)
if(i <= k)
{
sta[0][0][i] = A[i] / M1 , sta[0][1][i] = A[i] % M1;
sta[1][0][i] = B[i] / M1 , sta[1][1][i] = B[i] % M1;
}
else sta[0][0][i] = sta[0][1][i] = sta[1][0][i] = sta[1][1][i] = 0;
FFT(sta[0][0],1),FFT(sta[0][1],1),FFT(sta[1][0],1),FFT(sta[1][1],1);
static cplx rt[3][maxn];
for(int i=0;i<len;i++)
rt[0][i] = sta[0][1][i] * sta[1][1][i],
rt[1][i] = sta[0][1][i] * sta[1][0][i] + sta[1][1][i] * sta[0][0][i],
rt[2][i] = sta[0][0][i] * sta[1][0][i];
FFT(rt[0],-1),FFT(rt[1],-1),FFT(rt[2],-1);
for(int i=0;i<len;i++)
ret[i] = (llround(rt[0][i].r) % mod + 1ll * llround(rt[1][i].r) % mod * M1 %mod + 1ll * llround(rt[2][i].r) %mod * M2 %mod) % mod;
}
void Solve(int A[maxn],int B[maxn],int idx,int ret[maxn])
{
int p2=Pow(2,idx),sp2=1;
static int sta[2][maxn];
for(int i=0;i<len;i++,sp2=1ll*sp2*p2%mod)
sta[0][i] = 1ll * A[i] * sp2 % mod , sta[1][i] = B[i];
Move(sta[0],sta[1],ret);
}
int fac[maxn],inv[maxn],invf[maxn];
int dpret[maxn],baseret[maxn];
int main()
{
scanf("%I64d%I64d",&n,&k);
if(n > k){ printf("0\n");return 0; }
for(lgn=0;((k+1)<<1)>=(1<<lgn);lgn++);
len = 1<<lgn;
for(int i=1;i<len;i++) r[i] = (r[i>>1]>>1) | ((i&1) << (lgn-1));
for(int i=2,k=0;i<=len;i<<=1)
{
for(int j=0;j<i;j++)
w[k][j] = cplx(cos(j*PI/(i/2)),sin(j*PI/(i/2)));
k++;
}
fac[0] = fac[1] = inv[0] = inv[1] = invf[0] = invf[1] = 1;
for(int i=2;i<=k;i++)
fac[i] = 1ll * fac[i-1] * i % mod,
inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod,
invf[i] =1ll * invf[i-1]* inv[i] % mod;
for(int i=1;i<=k;i++) baseret[i] = invf[i];
dpret[0] = 1;
int idx = 1;
for(;n;n>>=1,Solve(baseret,baseret,idx,baseret),idx<<=1)
if(n&1)
Solve(dpret,baseret,idx,dpret);
int ans = 0;
for(int i=1;i<=k;i++)
ans = (ans + 1ll * dpret[i] * fac[k] % mod * invf[k-i] ) % mod;
printf("%d\n",(ans+mod)%mod);
}