2019 ICPC Asia Yinchuan Regional M Crazy Cake【ICPC银川M】【多项式】【生成函数】

题目大意

  圆周上均匀地放着n个点，你可以将若干点对之间的线段染黑，但必须保证线段不能在圆内相交。
  旋转同构，问方案数，多组数据。
  $T\le 10^5,n\le 10^6$

题解

  首先同构肯定是要考虑Burnside的，所以使用了Burnside后现在的问题变成：d个点不同构时的方案数（$d|n$)。
  但要注意到，如果$\frac{n}{d}\ge 2$，即循环节个数大于等于2个时的情况和$d=n$，即只有一个循环时的情况是不一样的，要分开考虑。

d=n

 $d=n$时相当是询问n个点不考虑同构时的方案数，由于此时最外围的n条边（i和i+1的边）对答案没有实质性影响，所以最后再乘以$2^n$，记为$f_n$为n个点只考虑内部连边时的方案数。
 此时考虑0号点的连边情况：
①0号点没有连边
 这种情况比较简单，因为0号点没有向其他点，所以把0号点删掉不影响答案，所以此时方案数为$2\ast f_{n-1}$,（*2是因为蓝色虚边可连可不连）。

②0号点有连边
 考虑0号点连出去的最小标号为 $i（1<i<n-1$），则0-i这条边将整个图形化成了两边，而两边各为一个子问题，所以此时的方案数为$2\cdot f_i\cdot f_{n-i+1}$（*2是因为红色实边可连可不连），但注意当i=2时并不需要*2，因为此时红色实边和虚边重合在一起了。.

综上所述，有$f_n=2f_{n-1}+\left(f_{n-1}+2{\sum }_{i=3}^{n-2}{f}_i\cdot f_{n-i+1}\right)(n\ge 3)$，特别地有$f_1=f_2=1$。
好那我们就可以开始大力FFT了
其实并不需要使用FFT，再推一下式子就可以了：
先整理一下：
$f_n=2f_{n-1}+2{\sum }_{i=2}^{n-2}{f}_i\cdot f_{n-i+1}\Rightarrow 5f_n=2{\sum }_{i=1}^n{f}_i\cdot f_{n-i+1}-f_{n-1}$

记$F(x)={\sum }_{n=0}^{\infty }{f}_{n+1}\cdot x^n$，则据上面的式子有$5F(x)=2F(x)\cdot F(x)-xF(x)$。
注意到$f_1=f_2=1$，则有$5F(x)=2F(x)\cdot F(x)-xF(x)+2x+3$。

把$F(x)$看成未知数，这是一个一元二次方程，解出可得$F(x)=\frac{x+5\pm \sqrt{x^2-6x+1}}{4}$。
根据x=0时有$f_1=1$，所以$F(x)=\frac{x+5-\sqrt{x^2-6x+1}}{4}$，则很容易就能算出$f_n.$

循环节个数大于等于2时

 首先有一点很清楚，相邻循环节是不可能相互连边的（除了两个端点外），否则就一定会出现边相交的情况。
 之所以这个情况和上面的情况不一样，是因为我们并不知道0号点是处在某个循环节的哪个位置（一开始我自以为0号点就是循环节的开头，WA到怀疑人生）。
 类似上面的处理，循环节周围的一圈边对答案没有本质性影响，最后再乘以$2^d$，记$g_d$为循环节长度为d，且不考虑外围边的方案数。
 还是按0号点的连边情况分类：
①0号点没有连边
和上面类似，0号点没有连边就可以将他删掉，退化成d-1的情况，所以此时方案数为$2\ast g_{d-1}$

②0号点向右边[2,d)连边
 假设0号点向右边连边的最小标号为$i(2\le i<d)$，则0-d及其对称边将整个图案划成了两个部分（如下图所示），红色框部分为子问题，蓝色框部分内部可以自由连边（即为$f_i$)，所以此时方案数为$2\cdot f_i\cdot g_{d-i+1}$（*2是因为蓝色实边可连可不连），但注意当i=2时并不需要*2，因为此时蓝色实边和虚边重合在一起了。.

③0号点没有向右边连边，但有向左边[2,d]连边
 假设0号点向左边连边的最大标号为$i(2\le i\le d)$，则0-d及其对称边将整个图案划成了两个部分（如下图所示），蓝色框部分为子问题，红色框部分内部可以自由连边（即为$f_{i+1}$)，所以此时方案数为$2\cdot f_{i+1}\cdot g_{d-i}$（*2是因为蓝色虚边可连可不连）

 但这里有个边界情况挺坑的，就是当i=d时。
 因为这里的条件是0号点没有向右边连边，即0号点的对称点（d号点）不能向左边连边，所以上图的红框部分应该修正为下图的橘框部分，即0不能包括d号点，此时方案数为$2f_d$,（*2是因为橘色虚边可连可不连）

综上所述，有$g_d=2\ast g_{d-1}$+$\left(f_2\cdot g_{d-1}+{\sum }_{i=3}^{d-1}2\cdot f_i\cdot g_{d-i+1}\right)$+$\left(2\cdot f_d+{\sum }_{i=2}^{d-1}2\cdot f_{i+1}\cdot g_{d-i}\right)$
特别地有$g_1=1,g_2=3$.

稍加整理后有：$5\cdot g_d=4{\sum }_{i=1}^d{f}_i\cdot g_{d-i+1}-g_{d-1}$

记$G(x)={\sum }_{n=0}^{\infty }{g}_{n+1}\cdot x^n$，则据上面的式子有$5G(x)=4F(x)\cdot G(x)-xG(x)$。
注意到$g_1=1,g_2=3$，则有$5G(x)=4F(x)\cdot G(x)-xG(x)+1$。

把$G(x)$看成未知数，可解得$G(x)=-\frac{1}{4F(x)-x-5}$，带入$F(x)=\frac{x+5-\sqrt{x^2-6x+1}}{4}$后有，$G(x)=\frac{1}{\sqrt{x^2-6x+1}}$，至此$g_d$也出来了。

题外话

 假设已知F(x)是一个m次多项式，即$F(x)={\sum }_{i=0}^m{f}_i\cdot x^i$，那么$G(x)=F^k(x)$的麦克劳林展开式怎么求呢？
 答：我会多项式快速幂！
 问：那k不是整数呢？比如k=1/2.
 答：我会FFT全家桶，先ln再exp！
 问：又或者m很小呢？比如m=2，你还用FFT吗？
 答：。。。

 这里介绍一种特别的方法：
 我们考虑对G(x)求导，则$G^{\prime }(x)=[F^k(x){]}^{\prime }=k{F}^{k-1}(x)\cdot F^{\prime }(x)$；
 两边乘以F(x)：$F(x)G^{\prime }(x)=k{F}^k(x)\cdot F^{\prime }(x)$;
 注意到$G(k)=F^k(x)$：$F(x)\cdot G^{\prime }(x)=kG(x)\cdot F^{\prime }(x)$.

 所以我们现在有这样一条式子：$F(x)\cdot G^{\prime }(x)=kG(x)\cdot F^{\prime }(x)$，感觉很对称有木有。
 系数对齐一下就好了：

 $[x^{n-1}]F(x)\cdot G^{\prime }(x)={\sum }_{i=0}^m{f}_i\cdot (n-i)g_{n-i}$
 $[x^{n-1}]kG(x)\cdot F^{\prime }(x)=k{\sum }_{i=1}^{m}i{f}_i\cdot g_{n-i}$

 $\left(F^{\prime }(x)={\sum }_{i=1}^{m}i{f}_i\cdot x^{i-1}，G^{\prime }(x)={\sum }_{i=1}^{\infty }i{g}_i\cdot x^{i-1}\right)$

 两边系数相等，即${\sum }_{i=0}^m{f}_i\cdot (n-i)g_{n-i}=k{\sum }_{i=1}^{m}i{f}_i\cdot g_{n-i}$；
 整理一下：$n{f}_0{g}_n=k{\sum }_{i=1}^{m}i{f}_i\cdot g_{n-i}-{\sum }_{i=1}^m{f}_i\cdot (n-i)g_{n-i}$

 即：$g_n=\frac{1}{f_0}\cdot \left(\frac{k+1}{n}{\sum }_{i=1}^{m}i{f}_i\cdot g_{n-i}-{\sum }_{i=1}^m{f}_i\cdot g_{n-i}\right)$

代码

#include<bits/stdc++.h>
#define maxn 1000050
#define modu 1000000007
using namespace std;
typedef long long LL;

const LL inv6=(modu+1)/6;

LL pw2[maxn],pw6[maxn],pwi6[maxn];
LL I[maxn];
LL f[maxn],g[maxn];
LL h[maxn];//full convex,not iso

LL c1[maxn]={modu-2,13,modu-13,2},c2[maxn];
LL c[maxn]={1,3},cc[maxn];
//(2*x^3-13*x^2+13*x-2+sqrt(x^2-6*x+1)*(2*x^2-7*x))/(12*x-2)

LL dp[maxn];

LL t[maxn];


// LL f2[maxn]={1,2,3},g2[maxn]={5,modu-14,modu-23,modu-64};

//f: ((x+5-sqrt(x^2-6*x+1))/4)
//sqrt(x^2-6x+1) : g[i]=(6i-9)*g[i-1]+(3-i)*g[i-2] ,g0=1,g1=-3;

//f2: f^2

//2 4 16 72 502 3604

bool noprime[maxn];
int p[maxn],pnum;
LL phi[maxn];

void init() {
    I[0]=I[1]=1;
    f[0]=f[1]=1;  g[0]=1,g[1]=modu-3;
    h[0]=1,h[1]=h[2]=2;
}

void predo()    {
    phi[1]=1;
    for (int i=2;i<maxn;++i)    {
        if (!noprime[i])
            phi[i]=i-1,p[pnum++]=i;
        for (int j=0;j<pnum&&i*p[j]<maxn;++j)   {
            noprime[i*p[j]]=1;
            if (i%p[j])
                phi[i*p[j]]=phi[i]*(p[j]-1);
            else    {
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
        }
    }
}

void add(LL &a,LL b)    {
    a=(a+b)%modu;
}

int main()  {
    init();
    predo();

    for (int i=pw2[0]=1;i<maxn;++i) pw2[i]=pw2[i-1]*2%modu;
    for (int i=pw6[0]=1;i<maxn;++i) pw6[i]=pw6[i-1]*6%modu;
    for (int i=pwi6[0]=1;i<maxn;++i) pwi6[i]=pwi6[i-1]*inv6%modu;

    for (int i=2;i<maxn;++i) I[i]=modu-modu/i*I[modu%i]%modu;
    for (int i=2;i<maxn;++i)  {
        g[i]=(6*i-9+modu)*g[i-1]+(3-i+modu)*g[i-2];
        g[i]%=modu;
        g[i]=g[i]*I[i]%modu;

        f[i]=(modu-g[i])*I[4]%modu;
    }
    for (int i=3;i<maxn;++i)    
        h[i]=f[i-1]*pw2[i]%modu;
    
    for (int i=2;i<maxn;++i)    {
        c[i]=(6*i-3)*c[i-1]+(1-i)*c[i-2];
        c[i]%=modu;
        c[i]=(c[i]+modu)*I[i]%modu;
    }

    for (int i=1;i<maxn;++i)
        cc[i]=c[i-1]*pw2[i]%modu;

    for (int i=1;i<maxn;++i)    {
        add(dp[i],h[i]);
        for (int j=i+i;j<maxn;j+=i)
            add(dp[j],cc[i]*phi[j/i]);
        dp[i]=dp[i]*I[i]%modu;
    }

    int T;
    scanf("%d",&T);
    while (T--) {
        int n;
        scanf("%d",&n);
        if (n<maxn)
            printf("%lld\n",dp[n]);
        else
            puts("0");
    }
    return 0;
}