### codeforces gym 102114 L Lost In The Echo

posted on 2019-03-11 22:08:35 | under 比赛题解 |

## Step1:只有加减法或者只有乘除法

$$\sum_{i=1}^{n}sign_i x_i$$

$$\prod_{i=1}^n x_i^{sign_i}$$

## Step2:表达式树与递归

$$\sum_{i=1}^{m}sign_i f_{a_{i}}$$

$$\prod_{i=1}^{m}f_{a_i}^ {sign_i}$$

## Step3:规避分配律

$$F(z)=\sum_{k}\frac{f(k)}{k!}$$

$$G(z)=\sum_{k}\frac{g(k)}{k!}$$

$$F(z)=\frac{1}{2}(exp(2G(z))-2G(z))$$

$$G(z)=exp(2F(z))-exp(F(z))-F(z)$$

$$F(z)=\sum_{k \geq 2}\frac{(2^{k-1})G^{k}(z)}{k!}$$

$$G(z)=\sum_{k \geq 2}\frac{(2^k-1)F^{k}(z)}{k!}$$

(其实和多项式exp的组合意义挺像的，懂exp组合意义的人应该不难看懂这个柿子)

## Step4 ：求解生成函数-分治fft

$exp(F(x))=\sum_{k}\frac{F^{k}(z)}{k!}$

$$E(F(z))=exp(F(z))-F(z)$$

$$F(z)=E(G(z))$$

$$G(z)=E(2F(z))-E(F(z))$$

$$E’(F(z))=F'(z)E'(F(z))$$

$$E'(x)=exp(x)$$

$$E'(F(z))=F'(z)(E(F(z))+F(z))$$

wow，这是一个优雅的"自己卷自己“类型的递归式,虽然现在已经是$9102$年了，不过您可能还是不会写自己卷自己类型的分治fft,这里简单科普一下

### 如何处理自己卷自己

$$F(n)=\sum_{i+j=n}P(i)Q(j)$$

$P,Q$的特点是除非你计算完了$F(0) \dots F(n-1)$否则你不可能求出这个数组的值

$$F(n)=\sum_{i+j=n}P(i)Q(j)$$

$$F(n)=\sum_{i=l}^{mid-1}P(i)Q(n-i)$$

## 一些细节

$$F(z)=E(G(z))$$

$$G(z)=E(2F(z))-E(F(z))$$

$$E'(F(z))=F'(z)(E'(F(z))+F(z))$$

$$E'(2F(z))=F'(z)(E'(2F(z))+F(z))$$

$$E'(G(z))=G'(z)(E'(G(z))+G(z))$$

## 去重

$$F(z)=\sum_{k \geq 2}\frac{G^{k}(z)}{k!}$$

$$G(z)=\sum_{k \geq 2}\frac{(2^k-1)F^{k}(z)}{k!}$$

$$T(n)=O(nlogn)+2T(n/2)=O(nlog^2n)$$

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;const int N=131072+10;typedef long long ll;typedef double db;
const ll mod=1e9+7;const int P=32768;const int Msk=32767;const int PP=(ll)P*P%mod;
const db pi=acos(-1.0);
struct cmp
{
db r;db v;
friend cmp operator +(cmp a,cmp b){return (cmp){a.r+b.r,a.v+b.v};}
friend cmp operator -(cmp a,cmp b){return (cmp){a.r-b.r,a.v-b.v};}
friend cmp operator *(cmp a,cmp b){return (cmp){a.r*b.r-a.v*b.v,a.r*b.v+a.v*b.r};}
void operator /=(const db& b){r/=b;v/=b;}
}rt[2][22][N];int rv[22][N];ll fac[N];ll inv[N];
inline void pre()
{
for(int d=1;d<=18;d++)
for(int i=1;i<(1<<d);i++)rv[d][i]=(rv[d][i>>1]>>1)|((i&1)<<(d-1));
for(int d=1,t=1;d<=18;d++,t<<=1)
for(int i=0;i<(1<<d);i++)rt[0][d][i]=(cmp){cos(i*pi/t),sin(i*pi/t)};
for(int d=1,t=1;d<=18;d++,t<<=1)
for(int i=0;i<(1<<d);i++)rt[1][d][i]=(cmp){cos(i*pi/t),-sin(i*pi/t)};
fac[0]=1;for(int i=1;i<N;i++)fac[i]=fac[i-1]*i%mod;
inv[0]=inv[1]=1;for(int i=2;i<N;i++)inv[i]=(mod-mod/i)*inv[mod%i]%mod;
}
inline void fft(cmp* a,int len,int d,int o)
{
for(int i=0;i<len;i++)if(i<rv[d][i])swap(a[i],a[rv[d][i]]);cmp* st;
for(int k=1,j=1;k<len;k<<=1,j++)
for(int s=0,i;s<len;s+=(k<<1))
for(i=s,st=rt[o][j];i<s+k;++i,++st)
{cmp a1=a[i+k]*(*st);a[i+k]=a[i]-a1;a[i]=a[i]+a1;}
if(o){for(int i=0;i<len;i++)a[i]/=len;}
}
namespace poly
{
cmp tr[N],tr1[N],tr2[N],tr3[N],tr4[N],tr5[N],tr6[N];
ll m13[N],m14[N],m23[N],m24[N];
inline void dbdft(ll* a,int len,int d,cmp* op1,cmp* op2)
{
for(int i=0;i<len;i++)tr[i]=(cmp){(db)(a[i]>>15),(db)(a[i]&Msk)};fft(tr,len,d,0);
tr[len]=tr[0];for(cmp *p1=tr,*p2=tr+len,*p3=op1;p1!=tr+len;++p1,--p2,++p3)
(*p3)=(cmp){p1->r+p2->r,p1->v-p2->v}*(cmp){0.5,0};
for(cmp* p1=tr,*p2=tr+len,*p3=op2;p1!=tr+len;++p1,--p2,++p3)
(*p3)=(cmp){p1->r-p2->r,p1->v+p2->v}*(cmp){0,-0.5};
}
inline void dbidft(cmp* a,int len,int d,ll* op1,ll* op2)
{
fft(a,len,d,1);for(int i=0;i<len;i++)op1[i]=(ll)(a[i].r+0.5)%mod;
for(int i=0;i<len;i++)op2[i]=(ll)(a[i].v+0.5)%mod;
}
inline void mul(ll* a,ll* b,ll* c,int len,int d)
{
dbdft(a,len,d,tr1,tr2);dbdft(b,len,d,tr3,tr4);
for(int i=0;i<len;i++)tr5[i]=tr1[i]*tr3[i]+(cmp){0,1}*tr2[i]*tr4[i];
for(int i=0;i<len;i++)tr6[i]=tr1[i]*tr4[i]+(cmp){0,1}*tr2[i]*tr3[i];
dbidft(tr5,len,d,m13,m24);dbidft(tr6,len,d,m14,m23);
for(int i=0;i<len;i++)c[i]=(m13[i]*PP+(m14[i]+m23[i])*P+m24[i])%mod;
}
}
namespace subsolve
{
ll a[N];ll b[N];ll tr[N];
# define coresolve(type1,type2)\
for(int i=l;i<mid;i++)a[i-l]=type1;for(int i=(len>>1);i<(len<<1);i++)a[i]=0;\
for(int i=0;i<len;i++)b[i]=type2;for(int i=len;i<(len<<1);i++)b[i]=0;\
poly::mul(a,b,tr,len<<1,d+1);\
for(int i=mid;i<r;i++)(def[i]+=tr[i-l])%=mod;
inline void c(ll* df,ll* ef,ll* f,ll* def,int l,int r,int d)
{
int mid=(l+r)>>1;int len=(r-l)>>1;
if(r-1-(l<<1)>=0)
{
for(int i=l;i<mid;i++)a[i-l]=(ef[i]+f[i])%mod;for(int i=len;i<(len<<1);i++)a[i]=0;
for(int i=l;i<mid;i++)b[i-l]=df[i];for(int i=len;i<(len<<1);i++)b[i]=0;
poly::mul(a,b,tr,r-l,d);
for(int i=mid;i<r;i++)(def[i]+=tr[i-(l<<1)])%=mod;
}if(l==0)return;len<<=1;
coresolve((ef[i]+f[i])%mod,df[i]);coresolve(df[i],(ef[i]+f[i])%mod);
}
}
namespace solver1
{
ll f[N],g[N],tf[N],df[N],dg[N],dtf[N],ef[N],eg[N],etf[N],
def[N],deg[N],detf[N];
inline void solve(int l,int r,int d)
{
if(r-l==1)
{
if(l==0){df[0]=dg[0]=1,dtf[0]=2;return;}
if(l==1){f[1]=g[1]=1,tf[1]=2;}(def[l]+=df[0]*(ef[l]+f[l]))%=mod;
(deg[l]+=dg[0]*(eg[l]+g[l]))%=mod;(detf[l]+=dtf[0]*(etf[l]+tf[l]))%=mod;
ef[l+1]=def[l]*inv[l+1]%mod;eg[l+1]=deg[l]*inv[l+1]%mod;etf[l+1]=detf[l]*inv[l+1]%mod;
f[l+1]=eg[l+1];g[l+1]=(etf[l+1]+mod-ef[l+1])%mod;tf[l+1]=f[l+1]*2%mod;
df[l]=f[l+1]*(l+1)%mod;dg[l]=g[l+1]*(l+1)%mod;dtf[l]=tf[l+1]*(l+1)%mod;
return;
}int mid=(l+r)>>1;solve(l,mid,d-1);
subsolve::c(df,ef,f,def,l,r,d);subsolve::c(dg,eg,g,deg,l,r,d);
subsolve::c(dtf,etf,tf,detf,l,r,d);solve(mid,r,d-1);
}
}
namespace solver2
{
ll f[N],g[N],tf[N],tg[N],df[N],dtg[N],dtf[N],ef[N],etg[N],etf[N],
def[N],detg[N],detf[N];
inline void solve(int l,int r,int d)
{
if(r-l==1)
{
if(l==0){df[0]=1;dtg[0]=dtf[0]=2;return;}
if(l==1){f[1]=g[1]=1;tf[1]=tg[1]=2;}(def[l]+=df[0]*(ef[l]+f[l]))%=mod;
(detg[l]+=dtg[0]*(etg[l]+tg[l]))%=mod;(detf[l]+=dtf[0]*(etf[l]+tf[l]))%=mod;
ef[l+1]=def[l]*inv[l+1]%mod;etg[l+1]=detg[l]*inv[l+1]%mod;etf[l+1]=detf[l]*inv[l+1]%mod;
f[l+1]=etg[l+1]*inv[2]%mod;g[l+1]=(etf[l+1]+mod-ef[l+1])%mod;
tf[l+1]=f[l+1]*2%mod;tg[l+1]=g[l+1]*2%mod;
df[l]=f[l+1]*(l+1)%mod;dtf[l]=tf[l+1]*(l+1)%mod;dtg[l]=tg[l+1]*(l+1)%mod;
return;
}int mid=(l+r)>>1;solve(l,mid,d-1);
subsolve::c(df,ef,f,def,l,r,d);subsolve::c(dtg,etg,tg,detg,l,r,d);
subsolve::c(dtf,etf,tf,detf,l,r,d);solve(mid,r,d-1);
}
}ll ans[N];
int main()
{
pre();solver1::solve(0,65536,16);solver2::solve(0,65536,16);
for(int i=1;i<65536;i++)
{
ll ret1=(solver1::f[i]+solver1::g[i]-(i==1))%mod;
ll ret2=(solver2::f[i]+solver2::g[i]-(i==1))*2%mod;
ans[i]=(ret2+mod-ret1)*fac[i]%mod;
}int T;scanf("%d",&T);
for(int i=1,n;i<=T;i++)scanf("%d",&n),printf("%I64d\n",ans[n]);
}