BZOJ 4228 ：Tibbar的后花园（EGF+牛顿迭代）

传送门

题意：
求$n$个点的无向图个数，满足：
任意三个点能互相到达的点$a,b,c$。
满足$dis(a,b),dis(b,c),dis(a,c)$不能三个同时相等。

题解：

容易发现每个联通块不是链就是大小不为3的倍数的环，那么列出$EGF$直接上多项式$exp$就好了，口胡完了。。

多项式$exp$的做法：分治FFT
首先设:

$$g(f)=\ln f - A$$

这是一个关于幂级数的函数
假设已经知道要求的函数的前$n$项，为$f_0$
对$g(f)$在$f_0$处泰勒展开：

$$g(f)=\sum_{i=0}^{\infty}\frac{g^{(i)}(f_0)}{i!}(f-f_0)^i$$
那么显然：
$$0=g(f)\equiv g(f_0)+g'(f_0)(f-f_0) \pmod {x^{2n}}$$
移项可得：
$$f\equiv f_0-\frac{g(f_0)}{g'(f_0)} \pmod{x^{2n}}$$
$$f\equiv f_0-\frac{\ln(f_0)-A}{\frac{1}{f_0}} \pmod{x^{2n}}$$
$$f\equiv f_0(1-\ln (f_0)+A)\pmod{x^{2n}}$$
然后可以这样迭代下去求了.

实际上，只要符合函数$g(f)$容易求导和复合的话，牛顿迭代都可以轻松求出要求的函数，举个例子：
求

$$f=\sqrt{A}$$ .
同样：
$$f\equiv f_0-\frac{g(f_0)}{g'(f_0)} \pmod{x^{2n}}$$
$$f\equiv f_0-\frac{f_0^2-A}{2f_0} \pmod{x^{2n}}$$

是不是非常简单？？
然而屡次口胡完$exp$，今天第一次打,居然是到卡常题（话说题解好像不是生成函数），然后就过不了了。。代码还是挂出来吧。。

#include <bits/stdc++.h>
#define RG register 
#define LL long long
using namespace std;
const int N=1e6+50;
const int Mod=1004535809;
const int G=3;
inline int power(int a,int b) {
    int rs=1;
    for(;b;b>>=1,a=(LL)a*a%Mod) if(b&1) rs=(LL)rs*a%Mod;
    return rs;
}
int n,k,fac[N],inv_n[N],invfac[N];
int f[N*8],g[N*8],ig[N*8],h[N*8],tp[N*8],pos[N*8],w[N*8];
inline void dft(int *a) {
    for(int i=1;i<k;i++) pos[i]=(i&1)?((pos[i>>1]>>1)^(k>>1)):(pos[i>>1]>>1);
    for(int i=1;i<k;i++) if(pos[i]>i) swap(a[i],a[pos[i]]);
    for(int bl=1;bl<k;bl<<=1) {
        int tl=bl<<1,wn=power(G,(Mod-1)/tl); w[0]=1;
        for(int i=1;i<bl;i++) w[i]=(LL)w[i-1]*wn%Mod;
        for(int bg=0;bg<k;bg+=tl)
            for(int j=0;j<bl;j++) {
                int &t1=a[bg+j],&t2=a[bg+j+bl],t3=(LL)t2*w[j]%Mod;
                t2=(t1-t3<0?t1-t3+Mod:t1-t3);
                t1=(t1+t3>=Mod?t1+t3-Mod:t1+t3);
            }
    }
}
inline void calc_inverse(int *a,int *b,int len) {
    if(len==1) {b[0]=power(a[0],Mod-2); return;}
    if(len!=1) calc_inverse(a,b,len>>1);
    k=len<<1; 
    for(int i=0;i<len;i++) tp[i]=a[i];
    for(int i=len>>1;i<len;i++) b[i]=0;
    for(int i=len;i<k;i++) tp[i]=(b[i]=0);
    dft(tp); dft(b);
    for(int i=0;i<k;i++) b[i]=(2ll*b[i]%Mod -(LL)tp[i]*b[i]%Mod*b[i]%Mod +Mod)%Mod;
    dft(b); reverse(b+1,b+k);
    const int inv=power(k,Mod-2);
    for(int i=0;i<len;i++) b[i]=(LL)b[i]*inv%Mod;
    for(int i=len;i<k;i++) b[i]=0;
}
inline void calc_ln(int *a,int *b,int len){
    calc_inverse(a,ig,len);
    k=len<<1; b[len-1]=0;
    for(int i=0;i<len-1;i++) b[i]=(LL)a[i+1]*(i+1)%Mod;
    for(int i=len;i<k;i++) b[i]=(ig[i]=0);
    dft(b); dft(ig);
    for(int i=0;i<k;i++) b[i]=(LL)b[i]*ig[i]%Mod;
    dft(b); reverse(b+1,b+k);
    const int inv=power(k,Mod-2);
    for(int i=0;i<len;i++) b[i]=(LL)b[i]*inv%Mod;
    for(int i=len-1;i>=1;i--) b[i]=(LL)b[i-1]*inv_n[i]%Mod;
    for(int i=len;i<k;i++) b[i]=0;
    b[0]=0; return;
}
inline void calc_exp(int *a,int *b,int len) {
    if(len==1) {b[0]=1; return;}
    if(len!=1) calc_exp(a,b,len>>1);
    calc_ln(b,h,len);  k=len<<1;
    for(int i=0;i<len;i++) tp[i]=(a[i]-h[i]+Mod)%Mod;
    tp[0]++;
    for(int i=len;i<k;i++) tp[i]=(b[i]=0);
    dft(b); dft(tp);
    for(int i=0;i<k;i++) b[i]=(LL)b[i]*tp[i]%Mod;
    dft(b); reverse(b+1,b+k);
    const int inv=power(k,Mod-2);
    for(int i=0;i<len;i++) b[i]=(LL)b[i]*inv%Mod;
    for(int i=len;i<k;i++) b[i]=0;
}
int main() {
    scanf("%d",&n); fac[0]=1; inv_n[1]=1;
    for(int i=1;i<=n;i++) fac[i]=(LL)fac[i-1]*i%Mod;
    for(int i=2;i<=n;i++) inv_n[i]=(Mod-(LL)(Mod/i)*inv_n[Mod%i]%Mod);
    invfac[n]=power(fac[n],Mod-2);
    for(int i=n-1;i>=1;i--) invfac[i]=(LL)invfac[i+1]*(i+1)%Mod;
    f[1]=1; f[2]=inv_n[2];
    for(int i=3;i<=n;i++) {
        f[i]=(LL)fac[i]*inv_n[2]%Mod;
        if(i%3)f[i]=(f[i]+(LL)fac[i-1]*inv_n[2]%Mod)%Mod;
        f[i]=(LL)f[i]*invfac[i]%Mod;
    }
    int len;for(len=1;len<=n;len<<=1);
    calc_exp(f,g,len);
    printf("%d\n",(LL)g[n]*fac[n]%Mod);
}