Codeforce 623 E. Transforming Sequence(FFT + 快速幂)

这个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  ,int k)
{
	int ret = 1;
	for(;k;k>>=1, =1ll* * %mod) if(k&1) ret=1ll*ret* %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], ret[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++)  ret[i] = invf[i];
	dpret[0] = 1;

	int idx = 1;
	for(;n;n>>=1,Solve( ret, ret,idx, ret),idx<<=1)
		if(n&1)
			Solve(dpret, ret,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);
}