codeforces gym 102114 L Lost In The Echo

shadowice1984 · 2019-03-11 22:08:35 · 个人记录

丧心病狂的分治fft………

出题人把模数开成10^9+7就是想让我去死……

前置芝士:快速傅里叶变换

哦，你不会快速傅里叶变换？出门左转你咕模板区，详细的题解多了去了

前置芝士:任意模数ntt(需要"三次变两次“优化)

蛤？9120年了你还在写三模ntt这种被历史淘汰的东西？赶紧去康康毛爷爷论文更新一下自己的多项式水平

本题题解

给出n个变量x_1.x_2,x_4 \dots x_n现在要在这些变量之间加上+,-,×,÷,()问最后可以得到多少种本质不同的表达式，两个表达式相同当且仅当他们化简之后是一样的，比如a+(b-c)/d和a-(b+c)/d其实是一样的,注意，你可以任意安排这些变量之间的位置，比如说x_3/x_1-x_2也是一个合法的表达式

那么我们先来考虑一些比较简单的情况，最后再一步一步拓展到我们需要求解的情况

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

那么考虑我们最后的表达式会长什么样，无论我们加的运算符和括号到底有多悬，最后我们总是可以把表达式写成这两种形式

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

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

其中sign_i不是正1就是-1,表示第i个变量对于答案的贡献

那么我们注意到除了sign值全部是-1的情况是无法构造的之外(因为我们只能在变量之间加运算符)其他情况总是可以构造出来的

因此如果仅仅允许加减法和乘除法，那么总的方案数就是十分简单的2^n-1

Step2:表达式树与递归

为了方便起见，我们设f_{n}代表一个含有n个不同变量的表达式

众所周知，计算代数表达式的规则是先小括号，再乘除，后加减

那么我们可以按照运算顺序将一个表达式表示为一颗二叉树，我们称之为表达式树

具体的构建方法是按照运算顺序找出最先需要计算的运算符,新建一个节点代表这次运算所产生的表达式，它的左儿子是运算符左侧的表达式，右儿子是运算符右侧的表达式

如此这般我们就搞出了一个二叉树出来，但是我们发现我们并不能统计所有不同的二叉表达式树来推知答案，因为加法和乘法具有交换律，两个表达式树不同可能仅仅是采用了不同的方式计算了加法和乘法

那么我们该怎么办呢？

答案是把二叉树转成多叉树，考虑我们计算这个表达式值的过程当中最后一定连续的进行了一段加减法或者是乘除法，我们将这些运算符拿出来，从而将表达式切成一些散段，不妨设一共切出了m段，第i段表达式含有a_{i}个变量，那么我们会发现最后的表达式其实可以写成这样的两种形式之一

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

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

其中sign_i表示第i个表达式对值的贡献，取值是±1

我们知道对于每一个f代表的表达式,里面含有的变量都是不同的，那么我们这个时候就可以把f当成一个子问题来看待了

此时我们就可以新建一个节点表示表达式树的根，然后将拆分出来的表达式递归下去建树,在父亲和第i个儿子之间连接一条类型为sign_i的边

如此这般我们就搞了一个多叉树出来，那么此时我们可以说两个表达式等价当且仅当他们的表达式树一样吗？

然而答案是否定的，我们会发现一种情况，对于表达式a-(b-c)/d和表达式a+(c-b)/d来讲，他们显然是完全一样的，但是表达式树却不同，出现这种问题主要是由于乘除法对加减法的分配律引起的，现在我们考虑如何规避这种情况

Step3:规避分配律

我们发现分配律好烦人啊，看一看我们能不能绕开它

我们发现分配律的恶心之处就在于对于一个表达式f_n我们可以很容易的构造出一个f'_n=-f_n,比如说表达式a-b和b-a

这样的话尽管f_n和f'_n是两个不同的表达式但是我们发现f'_n=-f_n,f_n=-f'n也就是说一旦他们作为表达式树上的一个孩子，并且头上顶了一个加号或者减号的时候就会导致我们算重复

那么怎么办呢？很简单我们认为f_n和f'_n是等价的，换句话说两种方案按照一种方案来计数

这样我们就无需担心分配律造成的影响了,如果这样的两种方案被视为一种方案，那么无论怎么变相反数，这种等价类都不可能变成另外一种等价类

那么我们怎么生成最后的答案呢?很简单，把求出来的方案数乘2就可以了

当然我们会漏掉一种情况那就是如果一个表达式只有加号那么我们将无法构造出一个和这个表达式互为相反数的表达式，我们需要减掉这种情况的方案数，这一部分之后再说如何计算

在去掉了分配律的影响之后，我们使用生成函数的方式来求解长度为n的表达式数目

我们在拼合两个表达式的时候可以重新分配标号产生一些新的方案，具体来讲如果我们把长度为n和m的表达式拼起来，那么有{n+m \choose n}种方式重新分配变量的标号

那么这种卷积里带组合数的计数问题我们通常上一个EGF(Exponential Generating Function,指数型生成函数)来解决，EGF会帮我们带上柿子中的组合数因此我们唯一要做的就是列式子即可了

不如设f(n)表示长度为n,.并且在求值时最后执行了加减法的表达式数目，g(n)表示长度为n，并且在求值的时候最后执行了乘除法的表达式数目，我们对于这两个数列构造他们的指数生成函数

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!}

解释一下组合意义就是，我们枚举这个加法类型的生成函数到底是由几个乘法类型的生成函数构成的，如果这个加法类型的生成函数恰好由k个表达式拼凑而成，那么给这些表达式添加正负号一共有2^k种方案，但是因为互为相反数的方案在同一个等价类中，因此仅仅额外产生2^{k-1}种拼接的方案，又因为各个表达式之间是无序的，因此需要除上k!

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

喔，在各类egf题的狂轰滥炸下你还不知道多项式exp的组合意义啊？，那你可能需要找个博客翻一翻了

第二个等式也是同理，无序的枚举所有拼凑情况之后，由于乘除法不受分配律的干扰，因此一共有2^k-1种方案

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

那么我们知道多项式exp的展开式其实是这个东西

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

我们发现这个东西和我们列出来的柿子很像，所以我们尝试定义另一种函数

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

那么回过头来看我们关于F和G的方程组，我们会发现它可以被这样重写

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

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

那么我们知道分治fft可以处理一种被称之为"我卷我自己"的递归式，也可以处理两个数组反复卷积，"左右互搏"类型的递归式

那么问题来了，现在要处理的递归式是"我exp我自己"……稍微冷静下让我们拿出一张草稿纸和削好的铅笔，进行一些机械的计算，我们就能处理这样的递归式

具体点来讲我们在分治fft的过程维护F,G,2F,E(F),E(G),E(2F)以及这些函数的导函数

重点就是利用关于E,E'之间的联系,对函数E求导，我们可以得到一个递归式

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,这里简单科普一下

如何处理自己卷自己

分治fft最重要的一点就是：把递归式改写成递推式

假设我们需要处理这样的一个数组F

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

假设我们现在需要求解数列的前$n$项，我们此时数列进行cdq分治，同时我们将要维护这样一个**循环不变式**（或者你把它理解成归纳假设也可以）那就是:**在函数$solve(l,r)$结束的时候$F(l) \dots F(r-1)$均计算完毕** 边界条件显然成立，在$solve(0,1)$结束的时候我们可以求出$F$的初始值现在我们需要维护这个归纳假设,怎么维护呢？假设我们正在进行$solve(l,r)$那么我们选择这个序列的中点$mid$,然后递归进行$solve(l,mid)

这样我们的手上就有了关于F(l) \dots F(mid-1)的全部信息了，我们尝试用这些信息去更新F(mid) \dots F(r-1)

我们知道

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

如果仅仅考虑l \leq i \leq mid-1 的部分时，我们会得到这样的柿子

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

此时你会发现一个非常精妙的事实，由于我们仅仅关心那些mid \leq n \leq r-1的n,因此在卷积柿子当中涉及的n-i这个下标是在(0, r-l)之间的，因此我们可以仅仅提取Q数组的前r-l项进行卷积，这样的话复杂度就被控制在了一个合理的范围内

同理你也可以使用相似的技巧处理出l \leq j \leq mid-1部分的贡献

最后是我们尚未求出P,Q数组的前r-l项的情况，这种情况直接把左半边的P,Q拿去卷积就好

如此这般我们就考虑完了左边对右边的贡献，此后计算右侧的F和左侧的P,Q再无关联，因此我们就可以递归分治右半边了，当右侧的递归也结束的时候，我们就计算出了区间当中F数组的值

一些细节

那么到此为止我们其实已经可以做这道题了，无非就是利用这五个方程进行递推而已

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))

~~大概维护那么十几个数组就可以做了.jpg~~

但是此时我们需要合理安排分治fft的归纳假设，

因为这里我们有一个非常致命的问题,当我们知道F(z)[z^n]的时候,我们只能推出F'(z)[z^{n-1}]

因此如果我们分治的时候打算用solve(l,r)求解F数组和G数组,那么我们就会尴尬的发现在solve函数结束之后我们不能得知关于l,r的全部信息，因为F'(z)[z^{r-1}]，G'(z)[z^{r-1}]我们是算不出来的,从而导致归纳失败

正确的做法是使用solve(l,r)求解F',G'数组，当我们分治到最底层的时候我们可以得知E'(F(z)),E'(G(z)),E'(2F(z))之类的数组的信息，此时我们用他们去推出下一项的E(F),E(G),E(2F)然后推出F,G，最后用下一项的F,G去求解这一项的F',G'这样我们的归纳假设才能成立

去重

我们在上面使用了指数生成函数求解了忽略负号时长度为n的表达式的方案数,我们知道将这个柿子乘2并不是答案,因为对于一个全部是加号的表达式我们无法找到一个和他互为相反数的表达式来配对

那么我们还需要再次设立生成函数，这回我们令F(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!}

当然还是自己列五个方程组然后用分治fft维护即可

看到这里你也知道这个做法的常数有多大了~~(没事,cf机子跑8000ms什么都跑的出来)~~

当然这里还是要装膜做样的分析下复杂度

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

蕴含的常数因子嘛……在每一层分治当中我们需要做24次fft……(滑稽)

上代码~

#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]);
}