orz ei

这种式子感觉自己完全推不出来 看到就不敢搞

原图感觉比较复杂，考虑进行一些变形。

建一张有重边无自环，边带编号的无向图.

对于原图中两条边$(i,u)(i,v)$，在新图中连接$(u,v)$，编号为$i$。

考虑一张没有编号的无向图，可以发现一种编号合法当且仅当每一条边与对应的点不相邻。

考虑容斥有哪些点与边相邻。

可以发现我们建的图能被表示为若干$>1$的环。于是对每个环考虑。

令$f_{i,j}$表示一个大小为$i$的环，有$j$条边没有讨论匹配。

可以发现

$$ f_{i,j}={1\over 2}(j-1)!i[x^iy^j]e^{(\sum_{i\ge 1}x^ii)y}={i[x^i]H^j\over 2j}\\ f_{i,0}=1 $$

这里$H={x\over (1-x)^2}$。

除2是要去掉翻转相同的情况。

那么我们的答案就是

$$ n!\sum_{i=0}^{n}(-1)^{n-i}i![x^ny^i]e^F \\ F=\sum_{i\ge 2,0\le j\le i}{f_{i,j}x^iy^j\over i!}(i-1)! $$

注意到$-ln(1-x)=\sum_{i\ge 1}{x^i\over i}$（被ei教育得到）

那么就有

$$ F=-{1\over 2}ln(1-Hy)-x-{1\over 2}xy-ln(1-x) \\ $$

$$ ans=n!\sum_{i=0}^{n}(-1)^{n-i}i![x^ny^i]{1\over e^{x+{1\over 2}xy}(1-x)\sqrt{1-{x\over (1-x)^2}y}} $$

一波推导（接下来5行来自ei）

$$ ans=n!\sum_{i=0}^{n}i![x^ny^i]{1\over e^{-x+{1\over 2}xy}(1+x)\sqrt{1-{x\over (1+x)^2}y}}\\ =n!\sum_{i=0}^{n}[x^ny^i]{e^x\over (1+x)}(\sum_{j\ge 0}{(-{1\over 2}xy)^j\over j!})(\sum_{k\ge 0}\binom{-{1\over 2}}{k}(-{xy\over (1+x)^2})^k) \\ =n![x^n]{e^x\over (1+x)}\sum_{j\ge 0}\sum_{k\ge 0}(-{1\over 2})^j\binom{-{1\over 2}}{k}\binom{j+k}{k}k!(-{x\over (1+x)^2})^kx^{j} \\ =n![x^n]{e^x\over (1+x)}\sum_{k\ge 0}\binom{-{1\over 2}}{k}k!(-{x\over (1+x)^2})^k\sum_{j\ge 0}\binom{j+k}{k}(-{1\over 2}x)^j \\ =n![x^n]{e^x\over (1+x)(1+{1\over 2}x)}\sum_{k\ge 0}\binom{-{1\over 2}}{k}k!(-{x\over (1+x)^2(1+{1\over 2}x)})^k $$

令$P(x)=\sum_{k\ge 0}\binom{-{1\over 2}}{k}k!(-x)^k=\sum_{k\ge 0}{(2k-1)!!\over 2^k}x^k$，这里认为$(-1)!!=1$

可以发现

$$ P'(x)=\sum_{k\ge 1}{(2k-1)!!\over 2^k}{2k+1-1\over 2}x^{k-1} \\ x^2P'(x)+({1\over 2}x-1)P(x)+1=0 \\ $$

令$H(x)=(1+x)^2(1+{1\over 2}x)$，带入$x\to {x\over H(x)}$得

$$ x^2H(P\circ{x\over H})'(x)+({1\over 2}x-H)(H-xH')(P\circ{x\over H})(x)+H(H-xH')=0 $$

于是就可以愉快的$O(n)$解决了。

代码:

#include <bits/stdc++.h>

int P=998244353;
long long p=1ll*P*P;
typedef std::vector<int>poly;
inline int mul(const int &a,const int &b){return 1ll*a*b%P;}
inline int sub(int a,const int &b){a-=b;return (a<0)?a+P:a;}
inline int add(int a,const int &b){a+=b;return(a>=P)?a-P:a;}
int qsm(int a,int b){
    int ans=1;
    while(b){
        if(b&1)ans=mul(ans,a);
        a=mul(a,a);
        b>>=1;
    }
    return ans;
}
long long tem[100001];
poly operator*(poly a,poly b){
    int size=a.size()+b.size()-1;
    poly c(size);
    memset(tem,0,size<<3);
    for(int i=0;i<a.size();i++)
        for(int j=0;j<b.size();j++)
            if((tem[i+j]+=1ll*a[i]*b[j]%P)>=p)tem[i+j]-=p;
    for(int i=0;i<size;i++)c[i]=tem[i]%P;
    return c;
} 
poly operator+(poly a,poly b){
    if(a.size()<b.size())a.resize(b.size());
    for(int i=0;i<b.size();i++)a[i]=add(a[i],b[i]);
    return a;
}
poly operator-(poly a,poly b){
    if(a.size()<b.size())a.resize(b.size());
    for(int i=0;i<b.size();i++)a[i]=sub(a[i],b[i]);
    return a;
}
poly diff(poly a){
    for(int i=0;i<a.size()-1;i++)a[i]=mul(a[i+1],i+1);
    a.pop_back();
    return a;
}
poly lift(poly a,int n){
    a.resize(a.size()+n);
    for(int i=a.size()-1;i>=n;i--)a[i]=a[i-n];
    for(int i=0;i<n;i++)a[i]=0;
    return a;
}
poly operator*(poly a,int n){
    for(int i=0;i<a.size();i++)a[i]=mul(a[i],n);
    return a;
}
struct fract{
    poly a,b;
    fract operator*(fract c){return (fract){a*c.a,b*c.b};}
    fract operator*(int n){return (fract){a*n,b};}
    fract operator+(fract c){return (fract){a*c.b+b*c.a,b*c.b};}
    fract operator-(fract c){return (fract){a*c.b-b*c.a,b*c.b};}
    fract operator-(){return (fract){poly()-a,b};}
    fract diff(){return (fract){::diff(a)*b-::diff(b)*a,b*b};}
    fract sum(){return (fract){b,b-a};}
    void print(){
        for(int i=0;i<a.size();i++)printf("%d ",a[i]);putchar('\n');
        for(int i=0;i<b.size();i++)printf("%d ",b[i]);putchar('\n');
    }
};
int _mul[10000001],invmul[10000001],n,m,g[10000001],f[10000001];
void div(int *a,int v){
    for(int i=1;i<=n;i++)a[i]=sub(a[i],mul(v,a[i-1]));
}
int C(int n,int m){if(n<m||m<0)return 0;return mul(mul(_mul[n],invmul[m]),invmul[n-m]);}
int main(){
    scanf("%d%d",&n,&P);
    _mul[0]=_mul[1]=invmul[0]=invmul[1]=1;
    for(int i=2;i<=n;i++)_mul[i]=mul(_mul[i-1],i);
    for(int i=2;i<=n;i++)invmul[i]=mul(P-P/i,invmul[P%i]);
    for(int i=2;i<=n;i++)invmul[i]=mul(invmul[i-1],invmul[i]);
    for(int i=0;i<=n;i++)g[i]=invmul[i];
    div(g,1);
    div(g,(P+1)>>1);
    poly H;
    H.push_back(1);
    H.push_back(2+((P+1)>>1));
    H.push_back(2);
    H.push_back((P+1)>>1);
    poly a=H;
    a=lift(a,2);
    poly b=poly()-H;
    b[1]=add(b[1],(P+1)>>1);
    poly c=(H-lift(diff(H),1));
    b=b*c;
    c=c*H;
    f[0]=1;
    for(int i=1;i<=n;i++){
        int sum=0;
        for(int j=0;j<a.size();j++)
            if(i+1-j>=0)sum=add(sum,mul(mul(i+1-j,f[i+1-j]),a[j]));
        for(int j=0;j<b.size();j++)
            if(i-j>=0)sum=add(sum,mul(f[i-j],b[j]));
        if(c.size()>i)sum=add(sum,c[i]);
        f[i]=mul(sub(0,qsm(b[0],P-2)),sum);
    }
    int ans=0;
    for(int i=0;i<=n;i++)ans=add(ans,mul(f[i],g[n-i]));
    printf("%d\n",mul(ans,_mul[n]));
}