数学记录

· · 算法·理论

Update on 2025.6.12：修缮博客，加了一些题。

有码版已经基本停止修缮，不过 code 没有问题。

前言：如果觉得有代码很烦，可以看这份，这一份不一定同步更新。

题单

PGF

ARC 里有很多好玩的数学题 ARC 记录。

有些题没有代码：Lust、ARC154F、ARC156E、ARC184E，Calc（但是有加强版的代码），生成树问题。

推导

第一类斯特林数

行

递推公式：{n\brack m}=(n-1){n-1\brack m}+{n-1\brack m-1}。

设第 i 行斯特林数的生成函数为：F_i(x)=\sum\limits_{j=0}^\infty{i\brack j}x^j，有

\begin{aligned} F_i(x)&=(i-1)F_{i-1}(x)+xF_{i-1}(x)\\ &=(x+i-1)F_{i-1}(x) \\ \therefore F_n(x)&=\prod_{i=0}^{n-1}(x+i)\\ &=x^{\overline{n}} \end{aligned}

所以可以选择暴力分治 NTT，也有简洁的单 \log。

考虑到：x^{\overline{2n}}=x^{\overline{n}}(x+n)^{\overline n}，于是开始考虑给定 f(x) 时，如何快速求 f(x+c)。

\begin{aligned} f(x)&=\sum_{i=0}^na_ix^i\\ \therefore f(x+c)&=\sum_{i=0}^na_i(x+c)^i\\ &=\sum_{i=0}^na_i\sum_{j=0}^ix^jc^{i-j}{i\choose j}\\ &=\sum_{j=0}^n\frac{x^j}{j!}\sum_{i=0}^n\frac{c^{i-j}}{(i-j)!}i!a_i\\ \end{aligned}

可以倒过来卷积。

(￣▽￣)"，代码鸽了,给个 O(n\log ^2 n) 的。

#include<bits/stdc++.h>

using namespace std;

namespace Fread {
    const int SIZE=1<<21;char buf[SIZE],*S,*T;
    inline char getchar() {if(S==T){T=(S=buf)+fread(buf,1,SIZE,stdin);if(S==T)return '\n';}return *S++;}
}
namespace Fwrite {
    const int SIZE=1<<21;
    char buf[SIZE],*S=buf,*T=buf+SIZE;
    inline void flush(){fwrite(buf,1,S-buf,stdout);S=buf;}
    inline void putchar(char c){*S++=c;if(S==T)flush();}
    struct POPOSSIBLE{~POPOSSIBLE(){flush();}}ztr;
}
#define getchar Fread :: getchar
#define putchar Fwrite :: putchar
namespace Fastio{
    struct Reader{
        template<typename T>
        Reader& operator >> (T& x) {
            char c=getchar();T f=1;
            while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}x=0;
            while(c>='0'&&c<='9'){x=x*10+(c-'0');c=getchar();}x*=f;
            return *this;
        }
        Reader(){}
    }cin;
    struct Writer{
        template<typename T>
        Writer& operator << (T x) {
            if(x==0){putchar('0');return *this;}
            if(x<0){putchar('-');x=-x;}
            static int sta[45];int top=0;
            while(x){sta[++top]=x%10;x/=10;}
            while(top){putchar(sta[top]+'0');--top;}
            return *this;
        }
        Writer& operator << (char c) {putchar(c);return *this;}
        Writer& operator << (const char* str){int cur=0;while(str[cur])putchar(str[cur++]);return *this;}
        Writer(){}
    }cout;
}
#define endl '\n'
#define cin Fastio :: cin
#define cout Fastio :: cout

const int mod=998244353;
struct modint {
    int val;
    static int norm(const int& x) { return x < 0 ? x + mod : x; }
    modint inv() const {
        int a = val, b = mod, u = 1, v = 0, t;
        while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v);
        return modint(u);
    }
    modint() : val(0) {}
    modint(const int& m) : val(norm(m)) {}
    modint(const long long& m) : val(norm(m % mod)) {}
    modint operator-() const { return modint(norm(-val)); }
    bool operator==(const modint& o) { return val == o.val; }
    bool operator<(const modint& o) { return val < o.val; }
    modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; }
    modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; }
    modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; }
    modint operator-(const modint& o) const { return modint(*this) -= o; }
    modint operator+(const modint& o) const { return modint(*this) += o; }
    modint operator*(const modint& o) const { return modint(*this) *= o; }
};

const modint g=3;
const modint invg=(mod+1)/3;
const int N=1e6+10;
int rev[N],n,A,B;
modint frac[N],inv[N];

void GetRev(int lim){
    int bit=log2(lim);
    for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<(bit-1)));
}

inline modint expow(modint x,int y){
    modint res=1;
    for(;y;y>>=1){if(y&1) res*=x;x*=x;}
    return res;
}

struct Poly{
    modint a[N];
    modint& operator [](const int &x){return a[x];}
    inline void Readin(const int len){for(int i=0;i<len;++i) cin>>a[i].val;}
    inline void out(const int len){for(int i=0;i<len;++i) cout<<a[i].val<<' ';putchar('\n');}
    inline void clear(const int len){for(int i=0;i<len;++i) a[i]=modint();}
    void NTT(int limit,int type){
        for(int i=0;i<limit;++i) if(i>rev[i]) swap(a[i],a[rev[i]]);
        for(int i=1;i<limit;i<<=1){
            modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1));
            for(int j=0;j<limit;j+=(i<<1)){
                modint G=1;
                for(int k=0;k<i;++k,G*=Gn){
                    modint x=a[j+k],y=G*a[i+j+k];
                    a[j+k]=x+y;
                    a[i+j+k]=x-y;
                }
            }
        }
        if(type==-1){
            modint inv=expow(limit,mod-2);
            for(int i=0;i<limit;++i) a[i]*=inv;
        }
    }
}f[18];

void mul(int dep,int l,int r){
    int mid=(l+r)>>1,lim=1;
    for(;lim<=(r-l+5);lim<<=1) ;
    f[dep].clear(lim);
    if(l==r){
        f[dep][0]=l;
        f[dep][1]=1;
        return ;
    }
    mul(dep-1,l,mid);
    for(int i=0;i<=mid-l+1;++i) f[dep][i]=f[dep-1][i];
    mul(dep-1,mid+1,r);GetRev(lim);

//  cout<<l<<' '<<r<<endl;
//  f[dep].out(10);
//  f[dep-1].out(10);

    f[dep].NTT(lim,1);f[dep-1].NTT(lim,1);
    for(int i=0;i<lim;++i) f[dep][i]*=f[dep-1][i];
    f[dep].NTT(lim,-1);

    for(int i=r-l+2;i<lim;++i) f[dep][i]=0;
    f[dep-1].clear(lim);

}

signed main(){
    return 0;
}

受到 Permutation（下面的例题）的启发，可以将上升幂翻转，做 P5850，也是 O(n\log n) 的。

列

常规生成函数，设单个轮换的生成函数为：G=\sum\limits_{i=1}^n \frac{x^i}{i}。

有 F=G^m，暴力快速幂或分治，O(n\log n) 或 O(n\log^2n)。

#include<bits/stdc++.h>

using namespace std;

namespace Fread {
    const int SIZE=1<<21;char buf[SIZE],*S,*T;
    inline char getchar() {if(S==T){T=(S=buf)+fread(buf,1,SIZE,stdin);if(S==T)return '\n';}return *S++;}
}
namespace Fwrite {
    const int SIZE=1<<21;
    char buf[SIZE],*S=buf,*T=buf+SIZE;
    inline void flush(){fwrite(buf,1,S-buf,stdout);S=buf;}
    inline void putchar(char c){*S++=c;if(S==T)flush();}
    struct POPOSSIBLE{~POPOSSIBLE(){flush();}}ztr;
}
#define getchar Fread :: getchar
#define putchar Fwrite :: putchar
namespace Fastio{
    struct Reader{
        template<typename T>
        Reader& operator >> (T& x) {
            char c=getchar();x=0;
            while(c<'0'||c>'9') c=getchar();
            while(c>='0'&&c<='9'){x=x*10+(c-'0');c=getchar();}
            return *this;
        }
        Reader(){}
    }cin;
    struct Writer{
        template<typename T>
        Writer& operator << (T x) {
            if(x==0){putchar('0');return *this;}
            if(x<0){putchar('-');x=-x;}
            static int sta[45];int top=0;
            while(x){sta[++top]=x%10;x/=10;}
            while(top){putchar(sta[top]+'0');--top;}
            return *this;
        }
        Writer& operator << (char c) {putchar(c);return *this;}
        Writer(){}
    }cout;
}
#define endl '\n'
#define cin Fastio :: cin
#define cout Fastio :: cout

const int g=3;
const int mod=167772161;
const int invg=(mod+1)/3;
struct modint {
    int val;
    static int norm(const int& x) { return x < 0 ? x + mod : x; }
    inline int Mod(const int &x){return x>=mod?x-mod:x;}
    modint inv() const {
        int a = val, b = mod, u = 1, v = 0, t;
        while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v);
        return modint(u);
    }
    modint() : val(0) {}
    modint(const int& m) : val(norm(m)) {}
    modint operator-() const { return modint(norm(-val)); }
    modint& operator+=(const modint& o) { return val = Mod(val+o.val),*this; }
    modint& operator-=(const modint& o) { return val = norm(val-o.val),*this; }
    modint& operator*=(const modint& o) { return val = static_cast<int>(1ll*val*o.val%mod),*this; }
    modint operator-(const modint& o) const { return modint(*this) -= o; }
    modint operator+(const modint& o) const { return modint(*this) += o; }
    modint operator*(const modint& o) const { return modint(*this) *= o; }
    friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; }
};

int modread() {
    char c=getchar();long long x=0;
    while(c<'0'||c>'9') c=getchar();
    while(c>='0'&&c<='9'){x=(x*10+(c-'0'))%mod;c=getchar();}
    return x;
 }

const int N=1e6+10;
int rev[N],lim,n,m,k;

inline modint expow(modint x,int y){
    modint res=1;
    for(;y;y>>=1){if(y&1) res*=x;x*=x;}
    return res;
}

inline void GetRev(const int lim){
    int bit=log2(lim)-1;
    for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<bit));
}

struct Init{
    const int maxx=5e5;
    modint inv[N],frac[N];
    Init(){
        inv[0]=frac[0]=1;
        for(int i=1;i<=maxx;++i) frac[i]=frac[i-1]*i;
        inv[maxx]=frac[maxx].inv();
        for(int i=maxx-1;i;--i) inv[i]=inv[i+1]*(i+1);
        for(int i=1;i<=maxx;++i) inv[i]*=frac[i-1];
    }
    modint operator ()(const int &c){return inv[c];}
}inv;

struct Poly{
    modint a[N];
    ~Poly(){}
    modint& operator [](const int x){return a[x];}
    inline void Readin(const int len){for(int i=0;i<len;++i) cin>>a[i].val;}
    inline void out(const int len){for(int i=0;i<len;++i) cout<<a[i].val<<' ';cout<<endl;}
    inline void integral(const int len){for(int i=len;~i;--i) a[i+1]=a[i]*inv(i+1);a[0]=0;}
    inline void derivation(const int len){for(int i=0;i<len;++i) a[i]=a[i+1]*(i+1);a[len]=0;}
    inline void NTT(const int limit,const int type){
        for(int i=0;i<limit;++i) if(i<rev[i]) swap(a[i],a[rev[i]]);
        for(int i=1;i<limit;i<<=1){
            modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1));
            for(int j=0;j<limit;j+=(i<<1)){
                modint G=1;
                for(int k=0;k<i;++k,G*=Gn){
                    modint x=a[j+k],y=a[i+j+k]*G;
                    a[j+k]=x+y;
                    a[i+j+k]=x-y;
                }
            }
        }
        if(type==1) return ;
        modint inv=expow(limit,mod-2);
        for(int i=0;i<limit;++i) a[i]*=inv;
    }
}x;

inline void getinv(Poly &y,Poly &x,const int len){
    static Poly s,tmp;
    int lim=2,bas=1;s[0]=x[0].inv();
    while(bas<(len<<1)){
        for(int i=0;i<bas;++i) tmp[i]=x[i];
        GetRev(lim);tmp.NTT(lim,1);s.NTT(lim,1);
        for(int i=0;i<lim;++i) s[i]=s[i]*2-tmp[i]*s[i]*s[i];
        s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0;
        bas<<=1;lim<<=1;
    }
    for(int i=0;i<len;++i) y[i]=s[i];
    for(int i=0;i<=bas;++i) s[i]=tmp[i]=0;
}

inline void ln(Poly &y,Poly &x,const int len){
    static Poly t1,t2;
    for(int i=0;i<len;++i) t1[i]=t2[i]=x[i];
    t1.derivation(len);getinv(t2,t2,len);
    for(lim=1;lim<(len<<1);lim<<=1) ;
    GetRev(lim);t1.NTT(lim,1);t2.NTT(lim,1);
    for(int i=0;i<lim;++i) t1[i]*=t2[i];
    t1.NTT(lim,-1);t1.integral(lim);
    for(int i=0;i<len;++i) y[i]=t1[i];
    for(int i=0;i<=lim;++i) t2[i]=t1[i]=0;
}

inline void exp(Poly &y,Poly &x,const int len){
    static Poly s,tmp,t;
    int lim=2,bas=1;s[0]=1;
    while(bas<(len<<1)){
        for(int i=0;i<bas;++i) tmp[i]=x[i];
        ln(t,s,bas);tmp[0]+=1;for(int i=0;i<bas;++i) tmp[i]-=t[i];
        GetRev(lim);tmp.NTT(lim,1);s.NTT(lim,1);
        for(int i=0;i<lim;++i) s[i]*=tmp[i];
        s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0;
        bas<<=1;lim<<=1;
    }
    for(int i=0;i<len;++i) y[i]=s[i];
    for(int i=0;i<=bas;++i) s[i]=tmp[i]=t[i]=0;
}

inline void expow(Poly &y,Poly &x,const int len,const int k){
    ln(x,x,len);for(int i=0;i<len;++i) y[i]=x[i]*k;exp(y,y,len);
}

int main(){
    cin>>n>>k;++n;
    for(int i=0;i<n-1;++i) x[i]=inv(i+1);
//  x.out(n);
    expow(x,x,n,k);modint tmp=(inv.frac[k]).inv();
    for(int i=n-1;i>=k;--i) x[i]=x[i-k]*inv.frac[i];
    for(int i=0;i<k;++i) x[i]=0;
    for(int i=k;i<n;++i) x[i]*=tmp;
    x.out(n);
    return 0;
}

第二类斯特林数

有一个不常用的递推式：${n \brace m}={n-1 \brace m-1}+m{n-1 \brace m}$ ，从组合意义上来讲是显然的。 但是这样太慢了，我们期望一个通项公式。 设当前有 $n$ 个元素， $f_i$ 为将其划分入 $i$ 个两两区分（无空集）的方案数，$g_i$ 为允许空集的方案数。 显然有： $$ \begin{aligned} g_i&=i^n\\ g_i&=\sum_{j=0}^i{i \choose j}f_j\\ \end{aligned} $$ 使用二项式反演： $$ \begin{aligned} f_i&=\sum_{j=0}^i{i \choose j}(-1)^{i-j}g_j\\ f_i&=\sum_{j=0}^i\frac{i!(-1)^{i-j}j^n}{j!(i-j)!}\\ \end{aligned} $$ 又因为第二类斯特林数要求互不区分，所以 ${n \brace i}=\frac{f_i}{i!}=\sum\limits_{j=0}^{i}\frac{(-1)^{i-j}j^n}{j!(i-j)!}$ ，我们得到了通项公式，而且可以拆分成卷积形式。 ### 第二类斯特林数·行 第二类斯特林数的同一行可以卷积 $O(n\log n)$ 做：${n \brace i}=\sum\limits_{j=0}^i\frac{j^n}{j!} · \frac{(-1)^{i-j}}{(i-j)!}$。 ```cpp #include<bits/stdc++.h> using namespace std; inline int read(){ int x=0,f=0; char c=getchar(); while(c<'0'||c>'9'){f|=(c=='-');c=getchar();} while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();} return f?-x:x; } const int mod=167772161; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; } }; const modint g=3; const int N=1e6+10; const modint invg=(mod+1)/3; modint frac[N],inv[N]; int rev[N],n; inline void GetRev(int lim){ int bit=log2(lim)-1; for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<bit)); } inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void NTT(int limit,int type){ for(int i=0;i<limit;++i) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=a[i+j+k]*G; a[j+k]=x+y; a[i+j+k]=x-y; } } } if(type==-1){ modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } } }x,y; int main(){ n=read();frac[0]=1;int lim; for(int i=1;i<N;++i) frac[i]=frac[i-1]*i; inv[N-1]=frac[N-1].inv(); for(int i=N-2;~i;--i) inv[i]=inv[i+1]*(i+1); for(int i=0;i<=n;++i) x[i]=(inv[i]*expow(i,n)); for(int i=0;i<=n;++i) y[i]=(inv[i]*((i&1)?-1:1)); for(lim=1;lim<=(n<<1);lim<<=1);GetRev(lim); x.NTT(lim,1);y.NTT(lim,1); for(int i=0;i<lim;++i) x[i]*=y[i]; x.NTT(lim,-1); for(int i=0;i<=n;++i) printf("%d ",x[i].val); return 0; } ``` ### 第二类斯特林数·列 同一列也简单（只会大常数 $O(n \log n)$）。 还是考虑让集合两两区分，单个集合可以放 $1$ 到 $\infty$ 个元素。（因为不允许空集）。 设 $m$ 个盒子的生成函数为 $F$ ，单个函数的生成函数为 $G$ ，则 $F=G^m$ 。 显然有 $G=\sum\limits_{i=1}^\infty \frac{x^i}{i!}=e^x-1$。所以 $F=(e^x-1)^m$ 。 所以可以直接多项式快速幂做，记得把系数还原回去。 $G$ 的常数项为 $0$ 所以需要位移。 ```cpp #include<bits/stdc++.h> using namespace std; const int mod=167772161; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } }; inline int read(){ int x=0,f=0; char c=getchar(); while(c<'0'||c>'9'){f|=(c=='-');c=getchar();} while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();} return f?-x:x; } const modint g=3; const modint invg=(mod+1)/3; const int N=1e6+3e5; modint inv[N],frac[N],Inv[N]; int rev[N],n,m,k; inline int modread(){ modint x=0; char c=getchar(); while(c<'0'||c>'9')c=getchar(); while(c>='0'&&c<='9'){x=(x*10+(int)c-'0');c=getchar();} return x.val; } void GetRev(int lim){ int bit=log2(lim); for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<(bit-1))); } inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } void Init(int len){ inv[1]=1;Inv[0]=frac[0]=1; for(int i=2;i<=len;++i) inv[i]=inv[mod%i]*(-mod/i); for(int i=1;i<=len;++i) frac[i]=frac[i-1]*i;Inv[len]=expow(frac[len],mod-2); for(int i=len-1;i;--i) Inv[i]=Inv[i+1]*(i+1); } struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void Readin(const int len){for(int i=0;i<len;++i) a[i]=read();} void out(const int len){for(int i=0;i<len;++i) printf("%d ",a[i].val);putchar('\n');} void derivation(const int len){for(int i=0;i<len;++i) a[i]=a[i+1]*(i+1);a[len-1]=0;} void integral(const int len){for(int i=len-1;i;--i) a[i]=a[i-1]*inv[i];a[0]=0;} void NTT(int limit,int type){ for(int i=0;i<limit;++i) if(i>rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; } } } if(type==-1){ modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } } }F; void getinv(Poly &y,Poly &x,const int len){ static Poly s,tmp; int lim=2,bas=1;s[0]=expow(x[0],mod-2); while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; GetRev(lim);s.NTT(lim,1);tmp.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=(s[i]*2)-s[i]*s[i]*tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; bas<<=1;lim<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=tmp[i]=0; } void ln(Poly &y,Poly &x,const int len){ static Poly t1,t2;int lim=1; for(int i=0;i<len;++i) t1[i]=x[i]; getinv(t2,t1,len);t1.derivation(len); for(;lim<(len<<1);lim<<=1) ;GetRev(lim); t1.NTT(lim,1);t2.NTT(lim,1);for(int i=0;i<lim;++i) t1[i]*=t2[i]; t1.NTT(lim,-1);t1.integral(len);for(int i=0;i<len;++i) y[i]=t1[i]; for(int i=0;i<=(len<<2);++i) t1[i]=t2[i]=0; } void exp(Poly &y,Poly &x,const int len){ static Poly s,tmp,t; int lim=2,bas=1;s[0]=1; while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; ln(t,s,bas);GetRev(lim);tmp[0]+=1; for(int i=0;i<bas;++i) tmp[i]-=t[i]; tmp.NTT(lim,1);s.NTT(lim,1); for(int i=0;i<lim;++i) s[i]*=tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; lim<<=1;bas<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=t[i]=tmp[i]=0; } inline void expow(Poly &x,int k,int len){ ln(x,x,len); for(int i=0;i<=len;++i) x[i]*=k; exp(x,x,len); } signed main(){ n=read()+1;k=read();Init(n<<3); for(int i=1;i<n;++i) F[i-1]=Inv[i];expow(F,k,n); for(int i=0;i<n;++i) if(i<k) printf("%d ",0); else printf("%d ",(F[i-k]*Inv[k]*frac[i]).val); return 0; } ``` # 刷题 ## [珍珠](https://www.luogu.com.cn/problem/P5401) 有点神秘，不过第一步转换想出来后就很好推。 用 $cnt_i$ 表示 $i$ 出现次数，则发现题目要求等价于 $\sum\limits_{i=1}^D \lfloor \frac{cnt_i}{2} \rfloor \geq m$ 成立，有如下转换： $$ \begin{aligned} \sum\limits_{i=1}^D \lfloor \frac{cnt_i}{2} \rfloor &\geq m \\ \sum\limits_{i=1}^D \frac{cnt_i-cnt_i\&1}{2} &\geq m \\ \sum\limits_{i=1}^D cnt_i-cnt_i\&1 &\geq 2m \\ \sum\limits_{i=1}^D cnt_i\&1 &\leq n-2m \end{aligned} $$ 也就是说，出现次数为奇数的数的个数小于等于 $n-2m$ 则合法。 先把 $n-2m \geq D$ 判掉，为 $D^{n}$ 。 记 $odd_i$ 表示刚好有 $i$ 个出现次数为奇数的方案数，发现恰好 $i$ 个这个限制很强，不是很好算，考虑通过二项式反演将限制弱化。 记 $f_i$ 表示我们钦定 $i$ 个数出现次数为奇数，剩下的数随便选的方案数，明显有如下关系式： $$ \begin{aligned} f_i&=\sum\limits_{j=i}^D {j \choose i} odd_j\\ odd_i&=\sum\limits_{j=i}^D{j \choose i}(-1)^{j-i}f_j\\ odd_i&=\frac{1}{i!}\sum\limits_{i=1}^D\frac{(-1)^{j-i}}{(j-i)!}f_j*j!\\ odd_i&=\frac{1}{i!}\sum\limits_{i=1}^D\frac{(-1)^{i-j}}{(-(i-j))!}f_j*j! \end{aligned} $$ 这是一个卷积形式，只要知道 $f$ 就可以 $O(D \log D)$ 算了。 目前就只需要算 $f$ 这个限制弱得多的东西了，可以想想 EGF 了（毕竟是有标号的）。 我们有： $$ \begin{aligned} e^x&=\langle 1,1,1,1,\dots \rangle\\ e^{-x}&=\langle 1,-1,1,-1,\dots \rangle\\ \frac{e^x-e^{-x}}{2}&=\langle 0,1,0,1,\dots \rangle。 \end{aligned} $$ 也就是选一个奇数的 EGF 了。 由定义可得，我们在算 $f_i$ 时要钦定 $i$ 个奇数，剩下的数随便选，这种选择策略的生成函数就是 $(\frac{e^x-e^{-x}}{2})^ie^{(d-i)x}$。 所以得到 $f_j$ 的计算式：$f_j={D \choose j}n![x^n](\frac{e^x-e^{-x}}{2})^ie^{(d-i)x}$ ，前面的组合数是因为我们有 ${D \choose j}$ 种钦定方式，$n!$ 是为了抵消 EGF 的 $\frac{1}{n!}$ 。 开始推式子： $$ \begin{aligned} f_j&={D \choose j}n![x^n](\frac{e^x-e^{-x}}{2})^ie^{(d-i)x}\\ f_j&=\frac{n!}{2^j}{D \choose j} [x^n](e^{(D-j)x}\sum_{k=0}^j{j \choose k} e^{kx}e^{(k-j)x}(-1)^{k-j})\ \ \ \ [二项式定理展开]\\ f_j&=\frac{n!}{2^j}{D \choose j} [x^n](\sum_{k=0}^j{j \choose k}e^{D+2(k-j)}(-1)^{k-j})\\ \end{aligned} $$ 由 $e^{ax}=\langle 1,ax,(ax)^2,\dots \rangle$ 将式子展开： $$ \begin{aligned} f_j&=\frac{n!}{2^j}{D \choose j} (\sum_{k=0}^j{j \choose k}\frac{(D+2(k-j))^n}{n!}(-1)^{k-j})\\ f_j&=\frac{D!}{2^j(D-j)!}\sum_{k=0}^j\frac{(D-2(j-k))^n}{k!(j-k)!}*(-1)^{j-k}\\ f_j&=\frac{D!}{2^j(D-j)!}\sum_{k=0}^j\frac{(D-2k)^n(-1)^{k}}{k!}*\frac{1}{(j-k)!}\\ \end{aligned} $$ 我们就又获得了一个卷积形式。 这里顺便提一下 $odd$ 的卷积计算，因为 $i-j$ 是负数，所以要位移 $D$ 位。 ```cpp #include<bits/stdc++.h> using namespace std; inline int read(){ int x=0,f=0; char c=getchar(); while(c<'0'||c>'9'){f|=(c=='-');c=getchar();} while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();} return f?-x:x; } const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } static constexpr int get_mod() { return mod; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; } }; const modint g=3; const int N=6e5+10; const modint invg=(mod+1)/3; modint frac[N],inv[N],ans; int rev[N],n,d,m; inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } void GetRev(const int lim){ int bit=log2(lim)-1; for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<bit)); } struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void NTT(const int limit,const int type){ for(int i=0;i<limit;++i) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y;a[i+j+k]=x-y; } } } if(~type) return ; modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } void clear(int len){for(int i=0;i<=len;++i) a[i]=0;} void out(int len){for(int i=0;i<=len;++i) printf("%d ",a[i].val);printf("\n");} }A,B,f,odd; void init(const int len){ frac[0]=1; for(int i=1;i<=len;++i) frac[i]=frac[i-1]*i;inv[len]=frac[len].inv(); for(int i=len-1;~i;--i) inv[i]=inv[i+1]*(i+1); } int main(){ init(200000); d=read();n=read();m=read();int lim; if(d<=n-2*m) return printf("%d\n",expow(d,n).val),0; if(n-2*m<0) return printf("%d\n",0),0; //计算f for(lim=1;lim<=(d<<1);lim<<=1) ;GetRev(lim); for(int i=0;i<=d;++i) A[i]=expow(d-i*2,n)*((i&1)?-1:1)*inv[i]; for(int i=0;i<=d;++i) B[i]=inv[i]; A.NTT(lim,1);B.NTT(lim,1); for(int i=0;i<=lim;++i) f[i]=A[i]*B[i]; f.NTT(lim,-1); A.clear(lim);B.clear(lim); for(int i=0;i<=d;++i) f[i]*=frac[d]*expow(expow(2,i),mod-2)*inv[d-i]*frac[i]; for(int i=d+1;i<=lim;++i) f[i]=0; //计算odd for(int i=0;i<=d;++i) A[d-i]=(i&1)?-inv[i]:inv[i]; // A.out(d);f.out(d); A.NTT(lim,1);f.NTT(lim,1); for(int i=0;i<lim;++i) B[i]=A[i]*f[i]; B.NTT(lim,-1); for(int i=0;i<=n-2*m;++i) ans+=B[i+d]*inv[i]; printf("%d\n",ans.val); return 0; } ``` ## [UVA12298](https://www.luogu.com.cn/problem/UVA12298) 发现只有四个花色，可以直接对每一个花色开一个 OGF ，$1$ 表示可选 ，$0$ 表示不可选，再直接卷起来就行了。 NTT 板题，只需用到多项式乘法，发现答案很大，所以要换大模数+__int128。 ```cpp #include<bits/stdc++.h> #define LL long long #define ll __int128 using namespace std; inline int read(){ int x=0; char c=getchar(); while(c<'0'||c>'9')c=getchar(); while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();} return x; } const ll mod=4179340454199820289ll; struct modint { ll val; static ll norm(const ll& x) { return x < 0 ? x + mod : x; } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint(const ll &m):val(norm(m%mod)){} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<ll>(val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << (long long)a.val; } }; const modint g=3; const int N=5e5+10; const modint invg=(mod+1)/3; int rev[N],flag[N],prime[N],cnt,n,m,k; void make_list(const int n){ for(int i=2;i<=n;++i){ if(!flag[i]) prime[++cnt]=i; for(int j=1;j<=cnt&&prime[j]<=n/i;++j){ flag[i*prime[j]]=1; if(!(i%prime[j])) break; } } } void GetRev(const int lim){ ll bit=log2(lim)-1; for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<bit)); } int Init(const int len){int lim;for(lim=1;lim<=len;lim<<=1);GetRev(lim);return lim;} inline modint expow(modint x,LL y){ modint res=1; for(;y;y>>=1){if(y&1ll) res*=x;x*=x;} return res; } struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} inline void out(const int len){for(int i=0;i<=len;++i) printf("%lld ",(long long)a[i].val);putchar('\n');} inline void NTT(const int limit,const int type){ for(int i=0;i<limit;++i) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; } } } if(~type) return ; modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } void clear(int len){for(int i=0;i<=len;++i) a[i]=0;} }x[5],y; inline int ys(const char c){return (c=='D')?3:((c=='C'?2:c=='H'));} signed main(){ make_list(100000); while((n=read())&&(m=read())){ k=read(); for(int i=1;i<=m;++i) if(flag[i]) for(int j=0;j<4;++j) x[j][i]=1; for(int i=1;i<=k;++i){ int X=0; char c=getchar(); while(c<'0'||c>'9')c=getchar(); while(c>='0'&&c<='9'){X=(X<<3)+(X<<1)+c-'0';c=getchar();} if(X<=m) x[ys(c)][X]=0; } int lim=Init(m<<2); for(int i=0;i<4;++i) x[i].NTT(lim,1); for(int i=0;i<lim;++i) y[i]=x[0][i]*x[1][i]*x[2][i]*x[3][i]; y.NTT(lim,-1);for(int i=n;i<=m;++i) printf("%lld\n",(LL)y[i].val);printf("\n"); for(int i=0;i<4;++i) x[i].clear(lim);y.clear(lim); } return 0; } ``` ## [The Child and Binary Tree](https://www.luogu.com.cn/problem/CF438E) 主要是用 $G$ 表示权值限制没想到。 我们设：若集合 $C$ 中有元素 $i$ ，则 $G[i]$ 为 $1$ ，否则为 $0$。 设 $F[i]$ 是权值为 $i$ 的神犇二叉树的个数，特殊地，$F[0]=1$ （方便）。 于是便有如下递推式（下标从 $0$ 开始也是为了方便 ）： $$ F[n]=\begin{cases} 1 & n=0\\ \sum\limits_{i=0}^n G[i]\sum\limits_{j=0}^{n-i}F[j]F[n-i-j] & otherwise\end{cases} $$ 发现这是一个三函数卷积的形式，那么该式子成立：$F=F^2G+1$ 。（$+1$ 是因为 $F[0]=1$ ）。 解一个一元二次方程，得到：$F=\frac{1\pm\sqrt{1-4G}}{2G}$ ，现在考虑去掉 $\pm$ 。 如下： $$ 令x趋于0\\ 在 F=\frac{1+\sqrt{1-4G}}{2G} 中，F趋于无穷。\\ 在 F=\frac{1-\sqrt{1-4G}}{2G} 中，F收敛。\\ $$ 所以要取负号。 再搞一搞分子有理化，做一下求逆和开根即可。 ```cpp #include<bits/stdc++.h> using namespace std; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } }; inline int read(){ int x=0,f=0; char c=getchar(); while(c<'0'||c>'9'){f|=(c=='-');c=getchar();} while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();} return f?-x:x; } const modint g=3; const modint I=911660635; const modint invg=(mod+1)/3; const modint inv2=(mod+1)/2; const int N=6e5+10; int rev[N],n,m; void GetRev(int lim){ int bit=log2(lim); for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<(bit-1))); } inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void Readin(const int len){for(int i=0;i<len;++i) a[i]=read();} void out(const int len){for(int i=0;i<len;++i) printf("%d ",a[i].val);putchar('\n');} void NTT(int limit,int type){ for(int i=0;i<limit;++i) if(i>rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; } } } if(type==-1){ modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } } }G,F; void getinv(Poly &y,Poly &x,const int len){ static Poly s,tmp; int lim=2,bas=1;s[0]=expow(x[0],mod-2); while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; GetRev(lim);s.NTT(lim,1);tmp.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=(s[i]*2)-s[i]*s[i]*tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; bas<<=1;lim<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=tmp[i]=0; } void Sqrt(Poly &y,Poly &x,const int len){ static Poly s,tmp,t; int lim=2,bas=1;s[0]=1; while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; getinv(t,s,bas);GetRev(lim); t.NTT(lim,1);tmp.NTT(lim,1);s.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=(s[i]*s[i]+tmp[i])*t[i]*inv2; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0;bas<<=1;lim<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=t[i]=tmp[i]=0; } signed main(){ n=read();m=read()+1; for(int i=1;i<=n;++i) G[read()]=1; for(int i=0;i<m;++i) G[i]=-G[i]*4;G[0]+=1; Sqrt(G,G,m);G[0]+=1;getinv(G,G,m); for(int i=0;i<m;++i) G[i]*=2; for(int i=1;i<m;++i) printf("%d\n",G[i].val); return 0; } ``` ## [WD与积木](https://www.luogu.com.cn/problem/P5162) 抽象题，不过不难。 看到这个经典的数据范围与题目形式，直接就想到了多项式工业（先试试再说）。 然后就试出来了。 考虑我们选当前这一层的生成函数为 $F(x)=\sum\limits_{i=1}^{\infty}\frac{1}{i!}x^i=e^x-1$。 然后考虑方案数为 $\sum\limits_{i=1}^\infty F(x)^i=\frac{1}{1-F(x)}=\frac{1}{2-e^x}$ 。 那么总和也就是加个权罢了：$\sum\limits_{i=1}^\infty iF(x)^i$ ，明显是几何级数的形式，为：$\frac{e^x-1}{(2-e^x)^2}$ 。 求逆后直接卷积，别忘了还原系数，否则会寄。 某人 $F$ 没清空我不说是谁。 ```cpp #include<bits/stdc++.h> using namespace std; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } }; inline int read(){ int x=0,f=0; char c=getchar(); while(c<'0'||c>'9'){f|=(c=='-');c=getchar();} while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();} return f?-x:x; } const modint g=3; const modint invg=(mod+1)/3; const modint inv2=(mod+1)/2; const int N=6e5+10; int rev[N],n,m,q[N],T; modint inv[N],frac[N]; void GetRev(int lim){ int bit=log2(lim); for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<(bit-1))); } inline int Init(const int m){int lim=1;for(;lim<=m;lim<<=1);GetRev(lim);return lim;} inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void Readin(const int len){for(int i=0;i<len;++i) a[i]=read();} void out(const int len){for(int i=0;i<len;++i) printf("%d ",a[i].val);putchar('\n');} void NTT(int limit,int type){ for(int i=0;i<limit;++i) if(i>rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; } } } if(type==-1){ modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } } }G,F,H; void getinv(Poly &y,Poly &x,const int len){ static Poly s,tmp; int lim=2,bas=1;s[0]=expow(x[0],mod-2); while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; GetRev(lim);s.NTT(lim,1);tmp.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=(s[i]*2)-s[i]*s[i]*tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; bas<<=1;lim<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=tmp[i]=0; } void init(const int n){ frac[0]=inv[0]=1; for(int i=1;i<=n;++i) frac[i]=frac[i-1]*i;inv[n]=expow(frac[n],mod-2); for(int i=n-1;i;--i) inv[i]=inv[i+1]*(i+1); } signed main(){ T=read();init(200000); for(int i=1;i<=T;++i) n=max(n,q[i]=read());++n; for(int i=1;i<n;++i) F[i]=-inv[i];F[0]+=1; getinv(H,F,n);G=H; for(int i=1;i<n;++i) F[i]=inv[i];F[0]=0; int lim=Init(n<<2); G.NTT(lim,1); F.NTT(lim,1); for(int i=0;i<lim;++i) G[i]=G[i]*G[i]*F[i]; G.NTT(lim,-1); for(int i=1;i<=T;++i) printf("%d\n",(G[q[i]]*expow(H[q[i]],mod-2)).val); return 0; } ``` ## [Calc加强版](https://www.luogu.com.cn/problem/P5850) [Calc简单版](https://www.luogu.com.cn/problem/P4463) 发现序列是有序的，很烦，故将其除以 $n!$ ，得到无序序列。 考虑生成函数，对值进行生成函数，$OGF$ 。 总的生成函数为 $\prod\limits_{i=1}^k(1+ix)$ ，直接做感觉死路一条，看一下能不能 $\ln+\exp$ 一套连招。 $\ln$ 后为 $\sum\limits_{i=1}^k\ln(1+ix)$，对其泰勒展开得到： $$ \sum_{i=1}^k\sum_{j=1}^\infty \frac{(-1)^{j+1}(ix)^{j}}{j}\\ \sum_{j=1}^\infty\frac{(-1)^{j+1}}{j}\sum_{i=1}^k(ix)^{j}\\ \sum_{j=1}^\infty\frac{(-1)^{(j+1)}}{j}x^j\sum_{i=1}^ki^j\\ $$ 然后后面是自然数幂和的形式，可以 $EGF$ 做。~~我今天才知道这玩意可以生成函数~~ 记 $g_i=\sum\limits_{j=1}^kj^i$ ，$G$ 为 $g$ 的生成函数，则有如下推导： $$ \begin{aligned} G(x)&=\sum_{i=0}^\infty(\sum_{j=0}^kj^i)\frac{x^i}{i!}\\ &=\sum_{j=0}^k\sum_{i=0}^\infty \frac{(jx)^i}{i!}\\ &=\sum_{j=0}^ke^{jx}=\frac{e^{(k+1)x}-1}{e^{x}-1} \end{aligned} $$ 如上，我们就可以通过求逆得到自然数幂和的值究竟是多少。 你问多项式常数项为 $0$ 怎么求逆？分子分母都有公因式 $x$ ，提取就好了。 然后这道题你就需要打一个求逆和 $exp$ 就能轻松水过了~~~。 某人忘记了最后要乘上 $n!$ 看样例才看出来，我不说是谁。 ~~不知道为啥第一行特判了才过~~ 知道了（数据好水）。 ```cpp #include<bits/stdc++.h> using namespace std; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } }; inline int read(){ int x=0,f=0; char c=getchar(); while(c<'0'||c>'9'){f|=(c=='-');c=getchar();} while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();} return f?-x:x; } const modint g=3; const int N=4e6+10; const modint I=911660635; const modint invg=(mod+1)/3; modint inv[N],Inv[N],frac[N]; int rev[N],n,m,k; inline int modread(){ modint x=0; char c=getchar(); while(c<'0'||c>'9')c=getchar(); while(c>='0'&&c<='9'){x=(x*10+(int)c-'0');c=getchar();} return x.val; } void GetRev(int lim){ int bit=log2(lim); for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<(bit-1))); } inline int Rev(const int len){int lim=1;for(;lim<=len;lim<<=1);GetRev(lim);return lim;} inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } void Init(int len){ inv[1]=1;Inv[0]=frac[0]=1; for(int i=2;i<=len;++i) inv[i]=inv[mod%i]*(-mod/i); for(int i=1;i<=len;++i) frac[i]=frac[i-1]*i;Inv[len]=expow(frac[len],mod-2); for(int i=len-1;i;--i) Inv[i]=Inv[i+1]*(i+1); } struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void Readin(const int len){for(int i=0;i<len;++i) a[i]=read();} void out(const int len){for(int i=0;i<len;++i) printf("%d ",a[i].val);putchar('\n');} void derivation(const int len){for(int i=0;i<len;++i) a[i]=a[i+1]*(i+1);a[len-1]=0;} void integral(const int len){for(int i=len-1;i;--i) a[i]=a[i-1]*inv[i];a[0]=0;} void NTT(int limit,int type){ for(int i=0;i<limit;++i) if(i>rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; } } } if(type==-1){ modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } } void clear(const int len){for(int i=0;i<=len;++i) a[i]=0;} }G,x,y; void getinv(Poly &y,Poly &x,const int len){ static Poly s,tmp; int lim=2,bas=1;s[0]=expow(x[0],mod-2); while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; GetRev(lim);s.NTT(lim,1);tmp.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=(s[i]*2)-s[i]*s[i]*tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; bas<<=1;lim<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=tmp[i]=0; } void ln(Poly &y,Poly &x,const int len){ static Poly t1,t2;int lim=1; for(int i=0;i<len;++i) t1[i]=x[i]; getinv(t2,t1,len);t1.derivation(len); for(;lim<(len<<1);lim<<=1) ;GetRev(lim); t1.NTT(lim,1);t2.NTT(lim,1);for(int i=0;i<lim;++i) t1[i]*=t2[i]; t1.NTT(lim,-1);t1.integral(len);for(int i=0;i<len;++i) y[i]=t1[i]; for(int i=0;i<=(len<<2);++i) t1[i]=t2[i]=0; } void exp(Poly &y,Poly &x,const int len){ static Poly s,tmp,t; int lim=2,bas=1;s[0]=1; while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; ln(t,s,bas);GetRev(lim);tmp[0]+=1; for(int i=0;i<bas;++i) tmp[i]-=t[i]; tmp.NTT(lim,1);s.NTT(lim,1); for(int i=0;i<lim;++i) s[i]*=tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; lim<<=1;bas<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=t[i]=tmp[i]=0; } signed main(){ Init(1000000); k=read();m=read()+1; for(int i=0;i<=m;++i) x[i]=Inv[i];x[0]-=1; for(int i=0;i<=m;++i) y[i]=expow(k+1,i)*Inv[i];y[0]-=1; for(int i=0;i<=m;++i) x[i]=x[i+1]; for(int i=0;i<=m;++i) y[i]=y[i+1]; getinv(x,x,m);int lim=Rev(m<<1);x.NTT(lim,1);y.NTT(lim,1); for(int i=0;i<lim;++i) G[i]=x[i]*y[i];G.NTT(lim,-1); for(int i=0;i<lim;++i) G[i]*=frac[i];G[0]=0; x.clear(lim);y.clear(lim); for(int i=1;i<m;++i) G[i]*=inv[i]*((i&1)?1:-1); exp(G,G,m);for(int i=1;i<m;++i) printf("%d\n",(G[i]*frac[i]).val); return 0; } ``` ## [城市规划](https://www.luogu.com.cn/problem/P4841) 考虑 EGF+NTT (毕竟模数都摆面前了，原根为 $3$ )。 发现连通这个限制比较强，考虑弱化或者转换形式，如下构造。 设 $g_i$ 表示点数为 $i$ 的**简单无向图**的数量，$f_i$ 为点数为 $i$ 的答案。 明显有 $g_i=2^{n \choose 2}$ 且有该转移方程：$g_n=\sum\limits_{i=1}^n{n-1 \choose i-1}f_{i}g_{n-i}$ ，可以考虑为枚举节点 $1$ 所在联通块的大小（所以 $n$ 和 $i$ 都要减 $1$）。 推一下： $$ \begin{aligned} g_n&=\sum\limits_{i=1}^n{n-1 \choose i-1}f_{i}g_{n-i}\\ \frac{g_n}{(n-1)!}&=\sum\limits_{i=1}^{n}\frac{f_i}{(i-1)!}·\frac{g_{n-i}}{(n-i)!}\\ \end{aligned} $$ 那么我们定义 $$ \begin{aligned} F&=\sum\limits_{i=1}^{\infty}\frac{f_i}{(i-1)!}x^i\\ G&=\sum_{i=0}^{\infty}\frac{2^{i\choose 2}}{i!}x^i\\ H&=\sum_{i=1}^{\infty}\frac{2^{i\choose 2}}{(i-1)!}x^i \end{aligned} $$ 有 $H=F*G$ ，则 $F=H*G^{-1}$，求逆即可。 某人直接把指数对 $mod$ 取模，我不说是谁。 ```cpp #include<bits/stdc++.h> using namespace std; inline int read(){ int x=0,f=0; char c=getchar(); while(c<'0'||c>'9'){f|=(c=='-');c=getchar();} while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();} return f?-x:x; } const int MOD=1004535809; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + MOD : x; } modint inv() const { int a = val, b = MOD, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % MOD)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % MOD, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % MOD), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; } }; const modint g=3; const int N=1e6+10; const modint invg=(MOD+1)/3; const modint inv2=(MOD+1)/2; modint inv[N],frac[N]; int rev[N],n; inline void GetRev(const int lim){ int bit=log2(lim)-1; for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<bit)); } inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void NTT(const int limit,const int type){ for(int i=0;i<limit;++i) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(MOD-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; } } } if(type==1) return ; modint inv=expow(limit,MOD-2); for(int i=0;i<limit;++i) a[i]*=inv; } void out(const int len){for(int i=0;i<=len;++i) cout<<a[i]<<' ';cout<<endl;} }F,G,H; void getinv(Poly &x,const int len){ static Poly tmp,s; int bas=1,lim=2;s[0]=expow(x[0],MOD-2); while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; GetRev(lim);tmp.NTT(lim,1);s.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=s[i]*2-s[i]*s[i]*tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; lim<<=1;bas<<=1; } for(int i=0;i<len;++i) x[i]=s[i]; for(int i=0;i<(len<<2);++i) s[i]=tmp[i]=0; } void init(const int len){ frac[0]=inv[0]=1; for(int i=1;i<=len;++i) frac[i]=frac[i-1]*i;inv[len]=frac[len].inv(); for(int i=len-1;i;--i) inv[i]=inv[i+1]*(i+1); } int main(){ n=read();init(n<<2);int lim; for(int i=1;i<=n;++i) H[i]=expow(2,(1ll*i*(i-1)/2)%(MOD-1))*inv[i-1]; for(int i=0;i<=n;++i) G[i]=expow(2,(1ll*i*(i-1)/2)%(MOD-1))*inv[i]; getinv(G,n+1);for(lim=1;lim<=(n<<1);lim<<=1) ;GetRev(lim); G.NTT(lim,1);H.NTT(lim,1);for(int i=0;i<lim;++i) F[i]=G[i]*H[i];F.NTT(lim,-1); printf("%d\n",(F[n]*frac[n-1]).val); return 0; } ``` ## [Quo Vadis](https://www.luogu.com.cn/problem/P6384) 有结论（归纳可证）消元后矩阵为： $$ A_{i,j}= \begin{cases} ij\phi(i) & i\mid j\\ 0 & i\nmid j \end{cases} $$ 操作二、四是送的，看操作三。 记 $sum(i)=\sum\limits_{j=1}^ij^2

消元前：

\begin{aligned} Ans&=\sum_{i=1}^n\sum_{i=1}^nij\gcd(i,j)\\ &=\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor\frac{n}{d} \rfloor}ij[i\perp j]\\ &=\sum_{i=1}^n d^3\sum_{k=1}^{\lfloor\frac{n}{d} \rfloor}\mu(k)k^2sum(\lfloor \frac{n}{kd}\rfloor)^2\\ &=\sum_{T=1}^nT^2\phi(T)sum(\lfloor\frac{n}{T}\rfloor)^2\\ \end{aligned}

考虑进行杜教筛，非线性求出。

记 S(i)=\sum\limits_{j=1}^ij^2\phi(j)，f(i)=i^2\phi(i)，g_i=i^2。

则：

\begin{aligned} (f*g)(n)&=\sum_{d|n}d^2\phi(d)\frac{n^2}{d^2}\\ &=n^3\\ \end{aligned} \\ \begin{aligned} g(1)S(n)&=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)S(\lfloor \frac{n}{i} \rfloor)\\ S(n)&=(\frac{n(n+1)}{2})^2-\sum_{i=2}^ni^2S(\lfloor \frac{n}{i} \rfloor)\\ \end{aligned}

消元后呢？如下：

\begin{aligned} Ans&=\sum_{i=1}^n\sum_{j=1}^nij\phi(i)[i\mid j]\\ &=\sum_{j=1}^nj\sum_{d\mid j}d\phi(d)\\ \end{aligned}

后面一坨是积性函数，线筛，在 p^k 处取值 \frac{p^{2k+1}+1}{p+1}。

#include<bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>

using namespace std;
using namespace __gnu_pbds;
typedef long long ll;

inline ll read(){
    ll x=0;
    char c=getchar();
    while(c<'0'||c>'9') c=getchar();
    while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();}
    return x;
}

const int mod=998244353;
struct modint {
    int val;
    static int norm(const int& x) { return x < 0 ? x + mod : x; }
    modint inv() const {
        int a = val, b = mod, u = 1, v = 0, t;
        while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v);
        return modint(u);
    }
    modint() : val(0) {}
    modint(const int m) : val(norm(m)) {}
    modint(const long long m) : val(norm(m % mod)) {}
    modint operator-() const { return modint(norm(-val)); }
    bool operator==(const modint& o) { return val == o.val; }
    bool operator<(const modint& o) { return val < o.val; }
    modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; }
    modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; }
    modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; }
    modint operator-(const modint& o) const { return modint(*this) -= o; }
    modint operator+(const modint& o) const { return modint(*this) += o; }
    modint operator*(const modint& o) const { return modint(*this) *= o; }
    friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; }
};

const int N=5e6+10;
const int lenth=5e6;
const int Max=5e5+10;
const modint inv2=(mod+1)/2;
const modint inv6=(mod+1)/6;
inline modint s2(const modint x){return x*(x+1)*(x*2+1)*inv6;}//二次前缀和 
inline modint s3(const modint x){return x*x*(x+1)*(x+1)*inv2*inv2;}
gp_hash_table < ll , modint > mp;
int cnt,prime[Max],flag[N],low[N],opt,T;
modint phi[N],p[N],f[N],g[N],inv[Max];
ll n,m;

void predo(const int n){
    flag[1]=1;
    for(int i=2;i<=n;++i){
        if(!flag[i]){
            prime[++cnt]=i;
            phi[i]=i-1;
            low[i]=i;
            p[i]=1ll*i*i%mod-i+1;
            inv[cnt]=modint(i+1).inv();
        }
        for(int j=1;j<=cnt&&prime[j]*i<=n;++j){
            int tmp=i*prime[j];
            flag[tmp]=1;
            if(i%prime[j]) {
                low[tmp]=prime[j];
                p[tmp]=p[i]*p[prime[j]];
                phi[tmp]=phi[i]*phi[prime[j]];
            }else {
                low[tmp]=low[i]*prime[j];
                if(i==low[i]) p[tmp]=((p[i]*(prime[j]+1)-1)*prime[j]*prime[j]+1)*inv[j];
                else p[tmp]=p[i/low[i]]*p[low[i]*prime[j]];
                phi[tmp]=phi[i]*prime[j];
                break;
            }
        }
        f[i]=phi[i];
    }
    f[1]=phi[1]=1;g[0]=1;p[1]=1;
    for(int i=1;i<=n;++i) g[i]=g[i-1]*phi[i]*i*i;
    for(int i=1;i<=n;++i) phi[i]=phi[i-1]+phi[i]*i*i;
    for(int i=1;i<=n;++i) p[i]=p[i-1]+p[i]*i;
}

inline modint S(const ll n){
    if(n<=lenth) return phi[n];
    if(mp[n].val) return mp[n];
    modint res=s3(n);
    for(ll i=2,j;i<=n;i=j+1){
        j=n/(n/i);
        res-=(s2(j)-s2(i-1))*S(n/i);
    }
    return mp[n]=res;
}

inline int PreThree(const ll n){
    modint res=0;
    for(ll i=1,j=0;i<=n;i=j+1){
        j=n/(n/i);
        res+=(S(j)-S(i-1))*s3(n/i);
    }
    return res.val;
} 
inline int SufThree(const int n){return p[n].val;}
ll gcd(const ll x,const ll y){return y?gcd(y,x%y):x;}

int main(){
    predo(5000000);
    T=read();
    bool flag=0;
    while(T--){
        opt=read();
        if(opt==1) flag=1;
        else if(opt==2) {
            n=read();m=read();
            printf("%d\n",(flag?((m%n)?0:(f[n]*(n%mod)*(m%mod)).val):(modint(n)*m*gcd(n,m)).val));
        }else if(opt==3) {
            n=read();
            printf("%d\n",flag?SufThree(n):PreThree(n));
        }else {
            n=read();
            printf("%d\n",g[n].val);
        }
    }
    return 0;
}

刷题计划

\begin{bmatrix} a_{k+2}\\ a_{k+1}\\ b_{k+2}\\ b_{k+1}\\ c_{k+2}\\ c_{k+1}\\ w^{k+1}\\ z^{k+1}\\ (k+1)^2\\ k+1\\ 1\\ \end{bmatrix} = \begin{bmatrix} p&q&1&0&1&0&0&0&r&t&1\\ 1&0&0&0&0&0&0&0&0&0&0\\ 1&0&u&v&1&0&1&0&0&0&0\\ 0&0&1&0&0&0&0&0&0&0&0\\ 1&0&1&0&x&y&0&1&0&1&2\\ 0&0&0&0&1&0&0&0&0&0&0\\ 0&0&0&0&0&0&w&0&0&0&0\\ 0&0&0&0&0&0&0&z&0&0&0\\ 0&0&0&0&0&0&0&0&1&2&1\\ 0&0&0&0&0&0&0&0&0&1&1\\ 0&0&0&0&0&0&0&0&0&0&1\\ \end{bmatrix} * \begin{bmatrix} a_{k+1}\\ a_{k}\\ b_{k+1}\\ b_{k}\\ c_{k+1}\\ c_{k}\\ w^{k}\\ z^{k}\\ k^2\\ k\\ 1\\ \end{bmatrix} \\ 得到如下 \\ \begin{bmatrix} a_{n}\\ a_{n-1}\\ b_{n}\\ b_{n-1}\\ c_{n}\\ c_{n-1}\\ w^{n-1}\\ z^{n-1}\\ (n-1)^2\\ n-1\\ 1\\ \end{bmatrix} = \begin{bmatrix} p&q&1&0&1&0&0&0&r&t&1\\ 1&0&0&0&0&0&0&0&0&0&0\\ 1&0&u&v&1&0&1&0&0&0&0\\ 0&0&1&0&0&0&0&0&0&0&0\\ 1&0&1&0&x&y&0&1&0&1&2\\ 0&0&0&0&1&0&0&0&0&0&0\\ 0&0&0&0&0&0&w&0&0&0&0\\ 0&0&0&0&0&0&0&z&0&0&0\\ 0&0&0&0&0&0&0&0&1&2&1\\ 0&0&0&0&0&0&0&0&0&1&1\\ 0&0&0&0&0&0&0&0&0&0&1\\ \end{bmatrix} ^{n-2} * \begin{bmatrix} 3\\ 1\\ 3\\ 1\\ 3\\ 1\\ w\\ z\\ 1\\ 1\\ 1\\ \end{bmatrix} 
#include<bits/stdc++.h>
#define int __int128

using namespace std;

namespace Fread {
    const int SIZE=1<<21;
    char buf[SIZE],*S,*T;
    inline char getchar() {
        if(S==T){
            T=(S=buf)+fread(buf,1,SIZE,stdin);
            if(S==T)return '\n';
        }
        return *S++;
    }
}
namespace Fwrite {
    const int SIZE=1<<21;
    char buf[SIZE],*S=buf,*T=buf+SIZE;
    inline void flush(){
        fwrite(buf,1,S-buf,stdout);
        S=buf;
    }
    inline void putchar(char c){
        *S++=c;
        if(S==T)flush();
    }
    struct POPOSSIBLE{
        ~POPOSSIBLE(){flush();}
    }ztr;
}
#define getchar Fread :: getchar
#define putchar Fwrite :: putchar
namespace Fastio{
    struct Reader{
        template<typename T>
        Reader& operator >> (T& x) {
            char c=getchar();
            T f=1;
            while(c<'0'||c>'9'){
                if(c=='-')f=-1;
                c=getchar();
            }
            x=0;
            while(c>='0'&&c<='9'){
                x=x*10+(c-'0');
                c=getchar();
            }
            x*=f;
            return *this;
        }
        Reader(){}
    }cin;
}
#define cin Fastio :: cin

int mod;
struct modint {
    int val;
    static int norm(const int& x) { return x < 0 ? x + mod : x; }
    modint() : val(0) {}

    modint(const int &m) : val(norm(m%mod)) {}
    modint(const signed &m) : val(norm(m%mod)) {}
    modint(const long long& m) : val(norm(m % mod)) {}
    modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; }
    modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; }
    modint operator+(const modint& o) const { return modint(*this) += o; }
    modint operator*(const modint& o) const { return modint(*this) *= o; }
    friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << (long long)a.val; }
};
const int len=11;
const int Max=12;
const int inf=1e9;
const int N=5e4+10;

struct Matrix{
    modint a[Max][Max];
    Matrix(){for(int i=0;i<Max;++i) for(int j=0;j<Max;++j) a[i][j]=0ll;}
    void init(){for(int i=1;i<=len;++i) a[i][i]=1ll;}
    inline modint* operator [](const int &x){return a[x];}
    inline Matrix operator *(const Matrix &t)const{
        Matrix c=Matrix();
        for(int i=1;i<=len;++i) 
            for(int k=1;k<=len;++k) 
                for(int j=1;j<=len;++j) 
                    c[i][j]+=a[i][k]*t.a[k][j];
        return c;
    }
    void out(){for(int i=1;i<=len;++i) {for(int j=1;j<=len;++j) cout<<a[i][j]<<' ';cout<<endl;}}
}f,g;
int n,p,q,r,t,u,v,w,x,y,z;

inline Matrix expow(Matrix x,int y){
    Matrix res;res.init();
    for(;y;y>>=1){
        if(y&1) res=res*x;
        x=x*x;
    }
    return res;
}

void initG(){
    g[1][1]=3;
    g[2][1]=1;
    g[3][1]=3;
    g[4][1]=1;
    g[5][1]=3;
    g[6][1]=1;
    g[7][1]=w;
    g[8][1]=z;
    g[9][1]=1;
    g[10][1]=1;
    g[11][1]=1;
}

void initF(){
    f[1][1]=p;f[1][2]=q;f[1][3]=1;f[1][5]=1;f[1][9]=r;f[1][10]=t;f[1][11]=1;
    f[2][1]=1;
    f[3][1]=1;f[3][3]=u;f[3][4]=v;f[3][5]=1;f[3][7]=1;
    f[4][3]=1;
    f[5][1]=1;f[5][3]=1;f[5][5]=x;f[5][6]=y;f[5][8]=1;f[5][10]=1;f[5][11]=2;
    f[6][5]=1;
    f[7][7]=w;
    f[8][8]=z;
    f[9][9]=1;f[9][10]=2;f[9][11]=1;
    f[10][10]=1;f[10][11]=1;
    f[11][11]=1;
}

signed main(){
    cin>>n>>mod>>p>>q>>r>>t>>u>>v>>w>>x>>y>>z;n-=2;
    initG();initF();
    f=expow(f,n)*g;
    cout<<"nodgd "<<f[1][1]<<endl;
    cout<<"Ciocio "<<f[3][1]<<endl;
    cout<<"Nicole "<<f[5][1]<<endl;
    return 0;
}

一个人的数论题

就习惯而言，将原题面的 d 称作 k。

设 F(n)=\sum\limits_{i=1}^ni^k，G(n)=\sum\limits_{i=1}^ni^k[\gcd(i,n)==1]。

答案为 G(n)。

有 F(n)=\sum\limits_{d|n}d^kG(\frac{n}{d}) 表示每次枚举出与 n 的最大公约数刚好为 d 的数的 i 次方和。

由莫比乌斯反演可知：G(n)=\sum\limits_{d|n}\mu(d)d^kF(\frac{n}{d})。

有如下式： $$ \begin{aligned} G(n)&=\sum\limits_{d|n}\mu(d)d^k\sum\limits_{i=0}^{k+1}(\frac{n}d)^iP_i\\ &=\sum_{i=0}^{k+1}P_in^i\sum_{d|n}\mu(d)d^{k-i}\\ \end{aligned} $$ 注意到 $\sum_{d|n}\mu(d)d^{k-i}$ 是积性函数，用质因数分解的形式搞定。 $$ \begin{aligned} G(n)&=\sum_{i=0}^{k+1}P_in^i\sum_{d|n}\mu(d)d^{k-i}\\ &=\sum_{i=0}^{k+1}P_in^i\prod_{c=1}^w\sum_{j=0}^{\alpha_c}\mu(p_c^j)p_c^{j*(k-i)}\\ &=\sum_{i=0}^{k+1}P_in^i\prod_{c=1}^w(1-p_c^{k-i}) \end{aligned} $$ 知道 $P_i$ 的话就能 $O(kw)$ 算了，可以选择高斯消元 $O(k^3)$做。 ```cpp #include<bits/stdc++.h> using namespace std; namespace Fread { const int SIZE=1<<21; char buf[SIZE],*S,*T; inline char getchar() { if(S==T){ T=(S=buf)+fread(buf,1,SIZE,stdin); if(S==T)return '\n'; } return *S++; } } namespace Fwrite { const int SIZE=1<<21; char buf[SIZE],*S=buf,*T=buf+SIZE; inline void flush(){ fwrite(buf,1,S-buf,stdout); S=buf; } inline void putchar(char c){ *S++=c; if(S==T)flush(); } struct POPOSSIBLE{ ~POPOSSIBLE(){flush();} }ztr; } #define getchar Fread :: getchar #define putchar Fwrite :: putchar namespace Fastio{ struct Reader{ template<typename T> Reader& operator >> (T& x) { char c=getchar(); T f=1; while(c<'0'||c>'9'){ if(c=='-')f=-1; c=getchar(); } x=0; while(c>='0'&&c<='9'){ x=x*10+(c-'0'); c=getchar(); } x*=f; return *this; } Reader& operator >> (char& c) { c=getchar(); while(c==' '||c=='\n')c=getchar(); return *this; } Reader& operator >> (char* str) { int len=0; char c=getchar(); while(c==' '||c=='\n')c=getchar(); while(c!=' '&&c!='\n'&&c!='\r'){ str[len++]=c; c=getchar(); } str[len]='\0'; return *this; } Reader(){} }cin; struct Writer{ template<typename T> Writer& operator << (T x) { if(x==0){putchar('0');return *this;} if(x<0){putchar('-');x=-x;} static int sta[45]; int top=0; while(x){sta[++top]=x%10;x/=10;} while(top){putchar(sta[top]+'0');--top;} return *this; } Writer& operator << (char c) { putchar(c); return *this; } Writer& operator << (char* str) { int cur=0; while(str[cur])putchar(str[cur++]); return *this; } Writer& operator << (const char* str) { int cur=0; while(str[cur])putchar(str[cur++]); return *this; } Writer(){} }cout; } #define cin Fastio :: cin #define cout Fastio :: cout const int mod=1e9+7; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator !(){return !val;} modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint& operator/=(const modint& o) { return *this *= o.inv(); } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } modint operator/(const modint& o) const { return modint(*this) /= o; } friend std::istream& operator>>(std::istream& is, modint& a) {return is >> a.val, a.val = norm(a.val % mod), is;} friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; } }; const int N=1010; modint p[N],P[N],y[N],t[N],ans,tot; int n,k,w; struct Vector{ modint v[N]; modint& operator [](const int &x){return v[x];} Vector& operator -=(Vector &t){for(int i=1;i<=(n<<1);++i) v[i]-=t[i];return *this;} Vector& operator *=(const modint &x){for(int i=1;i<=(n<<1);++i) v[i]*=x;return *this;} Vector operator *(const modint &x){return Vector(*this)*=x;} }c[N],d[N]; inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } bool Gauss(Vector *c,const int n){ Vector tmp; for(int i=1;i<=n;++i){ if(!c[i][i]){ int j=i; for(j=i+1;j<=n&&!c[j][i];++j) ; if(j>n) return false; swap(c[i],c[j]); } for(int j=1;j<=n;++j){ if(i==j||!c[j][i]) continue; modint delta=c[j][i]/c[i][i]; tmp=c[i]*delta;c[j]-=tmp; } } return true; } int main(){ cin>>k>>w;n=k+2;tot=1; for(int i=1;i<=w;++i) cin>>p[i]>>t[i]; for(int i=1;i<=w;++i) tot*=expow(p[i],t[i].val); for(int i=1;i<=n;++i){ c[i][1]=1; for(int j=2;j<=n;++j) c[i][j]=c[i][j-1]*i; y[i]=y[i-1]+expow(i,k); c[i][i+n]=1; } assert(Gauss(c,n)); for(int i=1;i<=n;++i) { modint tmp=c[i][i].inv(); for(int j=1;j<=n;++j) d[i][j]=c[i][j+n]*tmp; } for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) P[i-1]+=y[j]*d[i][j]; for(int i=0;i<n;++i){ modint res=P[i]*expow(tot,i); for(int c=1;c<=w;++c) res*=(-expow(p[c],mod-1+k-i)+1); ans+=res; } cout<<ans.val<<'\n'; return 0; } ``` ## [魔力环](https://www.luogu.com.cn/problem/P4916) 考虑 Burnside 引理的套路，记 $f(i)$ 为不考虑旋转的情况下的答案。 这里有一个额外的条件，那就是环的个数必定整除 $m$，否则无解（显然）。 计算式即为： $$ \begin{aligned} Ans&=\frac 1n\sum_{d=1}^nf(\gcd(d,n))[m|\gcd(d,n)]\\ &=\frac 1n\sum_{d|\gcd(n,m)}f(d)\phi(\frac {\gcd(n,m)}d)\\ \end{aligned} $$ 考虑计算 $f(i)$ ，相当于将环上 $n$ 个球中的 $m$ 个染成黑色，使得没有那一段黑色颜色段长度超过 $k$。 环上的计数是不方便的，于是考虑破链成环，枚举最后一个的后面和第一个的前面有多少球，剩下的序列计数，于是记 $g(n,m)$ 表示 $n$ 个球染 $m$ 个的答案，考虑对其容斥，即每次枚举钦定多少个颜色段长度超过 $k$，这个可以提前把 $n$ 减去 $i(k+1)$，再插板。 用 $h(n,m)$ 表示插板方案，则有如下计算式： $$ \begin{aligned} h(n,m)&={n+m-1\choose m-1}\\ g(n,m)&=\sum_{i=0}^{\lfloor \frac{n}{k+1} \rfloor} (-1)^i{m \choose i}h(n-i(k+1),m)\\ \end{aligned} $$ 不难发现计算 $f$ 是单次线性的。 时间复杂度：$O(\sigma(\gcd(n,m))+n)$。 ```cpp #include<bits/stdc++.h> using namespace std; namespace Fread { const int SIZE=1<<21; char buf[SIZE],*S,*T; inline char getchar() { if(S==T){ T=(S=buf)+fread(buf,1,SIZE,stdin); if(S==T)return '\n'; } return *S++; } } namespace Fwrite { const int SIZE=1<<21; char buf[SIZE],*S=buf,*T=buf+SIZE; inline void flush(){ fwrite(buf,1,S-buf,stdout); S=buf; } inline void putchar(char c){ *S++=c; if(S==T)flush(); } struct POPOSSIBLE{ ~POPOSSIBLE(){flush();} }ztr; } #define getchar Fread :: getchar #define putchar Fwrite :: putchar namespace Fastio{ struct Reader{ template<typename T> Reader& operator >> (T& x) { char c=getchar(); while(c<'0'||c>'9') c=getchar(); x=0; while(c>='0'&&c<='9'){ x=(x*10+(c-'0')); c=getchar(); } return *this; } }cin; struct Writer{ template<typename T> Writer& operator << (T x) { if(x==0){putchar('0');return *this;} if(x<0){putchar('-');x=-x;} static int sta[45]; int top=0; while(x){sta[++top]=x%10;x/=10;} while(top){putchar(sta[top]+'0');--top;} return *this; } Writer& operator << (char c) { putchar(c); return *this; } }cout; } //#define cin Fastio :: cin //#define cout Fastio :: cout const int MOD=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + MOD : x; } modint inv() const { int a = val, b = MOD, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % MOD)) {} modint operator-() const { return modint(norm(-val)); } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % MOD, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % MOD), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; } friend std::string tostring(const modint& a) { return std::to_string(a.val); } }; template < typename T > inline modint expow(const modint& a, const T& b) { modint x = a, res = 1; for (T p = b; p; x *= x, p >>= 1) if (p & 1) res *= x; return res; } const int N=2e5+10; int flag[N],phi[N],prime[N],cnt,n,m,k,d; modint frac[N],inv[N]; void init(const int len){ phi[1]=1;frac[0]=frac[1]=1; for(int i=2;i<=len;++i){ if(!flag[i]){ prime[++cnt]=i; phi[i]=i-1; } for(int j=1;j<=cnt&&i*prime[j]<=len;++j){ int tmp=i*prime[j]; flag[tmp]=1; if(i%prime[j]) phi[tmp]=phi[i]*phi[prime[j]]; else { phi[tmp]=phi[i]*prime[j]; break; } } frac[i]=frac[i-1]*i; } inv[len]=frac[len].inv(); for(int i=len-1;~i;--i) inv[i]=inv[i+1]*(i+1); } inline modint C(const int &n,const int &m){return (n>=m)?(frac[n]*inv[m]*inv[n-m]):0;} inline modint put(const int &n,const int &m){if(m<1) return 0;return C(n+m-1,m-1);} inline modint g(const int &n,const int &m){ modint res=0; for(int i=0;i<=m&&i<=(n/(k+1));++i) res+=C(m,i)*put(n-i*(k+1),m)*((i&1)?-1:1); return res; } inline modint f(const int &n,const int &m){ modint res=0; if(m<=k) return C(n,m); if(k==1) C(n-m+1,m)-C(n-m-1,m-2); for(int i=0;i<=min(k,n);++i) res+=g(m-i,n-m-1)*(i+1); return res; } int main(){ modint ans=0; cin>>n>>m>>k;init(n<<1);d=__gcd(n,m); if(m==n) return printf("%d\n",k>=n),0; if(!m) return printf("1\n"),0; for(int i=1;i<=d;++i) if(d%i==0) ans+=f(n/i,m/i)*phi[i]; cout<<ans*((modint(n)).inv())<<endl; return 0; } ``` ## [按位或](https://www.luogu.com.cn/problem/P3175) 考虑 $\min-\max$ 容斥（用在概率上好像是套路了）。~~dottle是不是讲过~~ 约定： $S$ 和 $T$ 都是存储位置集合，$S$ 代表该题的全集。 $\min(T)$ 表示集合中第一次有位置 $1$ 的时间（变成 $1$ 的时间的最小值）。 $\max(T)$ 所有位置都变成 $1$ 的时间（变成 $1$ 的时间的最大值）。 $\min-\max$ 容斥： $$ \min(S)=\sum\limits_{T\subset S}(-1)^{|T|-1}\max(T)\\ \max(S)=\sum\limits_{T\subset S}(-1)^{|T|-1}\min(T)\\ $$ 对于期望亦成立： $$ E(\min(S))=\sum\limits_{T\subset S}(-1)^{|T|-1}E(\max(T))\\ E(\max(S))=\sum\limits_{T\subset S}(-1)^{|T|-1}E(\min(T))\\ $$ 发现 $n\le 20$ 所以可以暴力容斥，现在只需要求出 $E(\min(T))$。 对于满足几何分布的随机离散变量，若其被选中的概率是 $P$ 则选中其的期望步数为 $\frac 1P$。 故现在只需求出对于一个集合 $T$，其元素被选中的概率是多少，问题有点困难，考虑正难则反。 问题即转化为：对于一个集合 $T$，选出一个被其完全包含的数的概率是多少，发现恰好是原概率数组经由 FMT 后得到的数组。~~就是 or 卷积的 FWT 数组~~。 理由是考虑 FWT 的式子：$fwt_i=\sum\limits_{j\subset i}a_j$。 ```cpp #include<bits/stdc++.h> using namespace std; const int N=1e6+3e5; const double eps=1e-8; double f[N],ans; int n,p; int main(){ ios::sync_with_stdio(0); cin.tie(0);cout.tie(0); cin>>p;n=1;for(int i=1;i<=p;++i) n<<=1; for(int i=0;i<n;++i) cin>>f[i]; for(int i=1;i<n;i<<=1) for(int j=0;j<n;j+=(i<<1)) for(int k=0;k<i;++k) f[i+j+k]+=f[j+k]; --n; for(int i=1;i<=n;++i){ if((1-f[i^n])<eps) return cout<<"INF\n",0; double tmp=1/(1-f[i^n]);ans+=((__builtin_popcount(i)&1)?1:-1)*tmp; } cout<<fixed<<setprecision(10)<<ans<<endl; return 0; } ``` ## 喂鸽子 是集训队作业的题，改了一下背景放在了模拟赛。 题意：有 $n$ 个初值为 $0$ 变量，你每次会等概率选择一个变量使其加一，问期望多少次使得每一个数大于等于 $k$，模 $998244353$。 $$n\leq 50,k\leq1000$$。 $\min-\max$ 容斥萌萌题。 然后可以开始考虑设 $F_i$ 表示有 $i$ 个变量，第一个大于等于 $k$ 的数期望出现时间。 有 $Ans=\sum\limits_{i=1}^nF_i*\frac ni*{n \choose i}*(-1)^{i+1}$，考虑怎么求 $F_i$。 设 $G_{i,j}$ 表示有 $i$ 个变量，操作了 $j$ 次还没有数大于等于 $k$ 的方案数，有如下： $$ F_i=i\sum\limits_{j=1}^{(i-1)(k-1)}\frac{{j+k-1\choose k-1}*G_{i-1,j}*(j+k)}{i^{j+k}} $$ 该式子的意义是先确定哪一个数达到 $k$（乘 $i$ ），再枚举的 $j$ 的意义为给其他数加了多少个。贡献就是合法方案除以总方案的意思。 问题落在了求 $G_{i,j}$，可以考虑生成函数。 设指数型生成函数 $g=\sum\limits_{i=0}^{k-1}\frac {x^i}{i!}$，则有 $G_{i,j}=(g^{i})[x^j]*j!$，暴力 NTT。 来考虑一下组合意义：在有标号的情况下分配操作次数得到答案。 于是此题有了 $O(n^2k\log nk)$ 的暴力做法。 ```cpp #include<bits/stdc++.h> using namespace std; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; } }; const modint g=3; const modint invg=(mod+1)/3; const int N=2e5+10; int rev[N],n,k,len; modint G[51][N],frac[N],inv[N]; void GetRev(int lim){int bit=log2(lim);for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<(bit-1)));} inline modint expow(modint x,int y){modint res=1;for(;y;y>>=1){if(y&1) res*=x;x*=x;}return res;} struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void out(const int len){for(int i=0;i<=len;++i) printf("%d ",a[i].val);putchar('\n');} void NTT(int limit,int type){ for(int i=0;i<limit;++i) if(i>rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; } } } if(type==-1){ modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } } }x,y; inline void init(const int len){ frac[0]=1;for(int i=1;i<=len;++i) frac[i]=frac[i-1]*i; inv[len]=frac[len].inv();for(int i=len-1;~i;--i) inv[i]=inv[i+1]*(i+1); } inline modint C(const int n,const int m){return frac[n]*inv[m]*inv[n-m];} inline modint F(const int i){ if(i==1) return k; modint Ans=0; for(int j=0;j<=(i-1)*(k-1);++j) Ans+=expow(i,1ll*(mod-2)*(j+k)%(mod-1))*C(j+k-1,k-1)*G[i-1][j]*(j+k); return Ans*i; } signed main(){ cin>>n>>k;len=n*k;init(len<<1);int lim; for(lim=1;lim<=len;lim<<=1) ;GetRev(lim); for(int i=0;i<=k-1;++i) x[i]=inv[i];x.NTT(lim,1);y[0]=1; for(int i=1;i<=n;++i){ y.NTT(lim,1); for(int i=0;i<lim;++i) y[i]*=x[i]; y.NTT(lim,-1);for(int i=len;i<lim;++i) y[i]=0; // y.out(len); for(int j=0;j<=len;++j) G[i][j]=y[j]*frac[j]; } modint ans=0; for(int i=1;i<=n;++i) ans+=F(i)*n*expow(i,mod-2)*C(n,i)*((i&1)?1:-1); cout<<ans.val<<endl; return 0; } ``` ## [Card](https://www.luogu.com.cn/problem/CF1278F) 和 [加强](https://www.luogu.com.cn/problem/P6031) 记 $p=\frac{1}{m}$，也就是抽到王牌的概率。 则答案等价于 $\sum\limits_{i=1}^n{n\choose i}(1-p)^{n-i}p^ii^k$。 $$ \begin{aligned} Ans&=\sum\limits_{i=0}^n{n\choose i}(1-p)^{n-i}p^ii^k\\ &=\sum\limits_{i=0}^n{n\choose i}(1-p)^{n-i}p^i\sum_{j=0}^k{i\choose j}j!{k\brace j}\\ &=\sum_{j=0}^k{k\brace j}j!\sum_{i=0}^n{n\choose i}{i\choose j}(1-p)^{n-i}p^i\\ &=\sum_{j=0}^k{k\brace j}j!\sum_{i=0}^n\frac{n!}{(n-i)!j!(i-j)!}(1-p)^{n-i}p^i\\ \end{aligned} $$ 后面部分： $$ \begin{aligned} &=\sum_{i=0}^n\frac{n!}{(n-i)!j!(i-j)!}(1-p)^{n-i}p^i\\ &=\frac{n!}{j!(n-j)!}\sum_{i=0}^n\frac{(n-j)!}{(n-i)!(i-j)!}(1-p)^{n-i}p^i\\ &={n\choose j}\sum_{i=0}^n{n-j\choose i-j}(1-p)^{n-i}p^i\\ &={n\choose j}p^j\sum_{i=0}^n{n-j\choose i-j}(1-p)^{n-i}p^{i-j}\\ &={n\choose j}p^j\sum_{i=0}^{n-j}{n-j\choose i}(1-p)^{n-i-j}p^i\\ &={n\choose j}p^j \end{aligned} $$ 代入原式得： $$ \begin{aligned} \sum_{j=0}^k{k\brace j}j!{n\choose j}p^j \end{aligned} $$ $$ \begin{aligned} \because {k\brace j}&=\sum_{i=0}^j\frac{i^k(-1)^{j-i}}{i!(j-i)!}\\ \therefore Ans&=\sum_{j=0}^k{n\choose j}j!p^j\sum_{i=0}^j\frac{i^k(-1)^{j-i}}{i!(j-i)!}\\ &=\sum_{j=0}^kp^j\sum_{i=0}^j\frac{n!}{(n-j)!i!(j-i)!}i^k(-1)^{j-i}\\ &=\sum_{j=0}^kp^j\sum_{i=0}^j{n\choose i}{n-i\choose j-i}i^k(-1)^{j-i}\\ &=\sum_{i=0}^k{n\choose i}i^k\sum_{j=i}^k{n-i\choose j-i}(-1)^{j-i}p^j\\ &=\sum_{i=0}^k{n\choose i}i^kp^i\sum_{j=0}^{k-i}{n-i\choose j}(-1)^jp^j\\ \end{aligned} $$ 记 $f(i)=\sum\limits_{j=0}^j{n+i-k\choose j}(-1)^jp^j$，$F(i)={n+i-k\choose i}(-p)^i$。 则有： $$ \begin{aligned} Ans&=\sum_{i=0}^k{n\choose i}i^kp^if(k-i)\\ f(i+1)-f(i)&=F(i+1)+\sum_{j=0}^i({n+i-k+1\choose j}-{n+i-k\choose j})(-p)^j\\ &=F(i+1)+[n+i=k-1]+\sum_{j=1}^i{n+i-k\choose j-1}(-p)^j\\ &=F(i+1)+[n+i=k-1]+pF(i)-pf(i)\\ f(i)&=F(i)+pF(i-1)-(p-1)f(i-1)+[n+i=k]\\ \end{aligned} $$ 由此实现 $O(k)$ 递推，$i^k$ 可以线性筛出来，毛估很卡常。 暴力： ```cpp #include<bits/stdc++.h> using namespace std; const int N=5010; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } }S[N][N],ans,c,p,tmp,t2; int n,m,k; int main(){ cin>>n>>m>>k;S[0][0]=c=tmp=t2=1;p=modint(m).inv(); for(int i=1;i<=k;++i) for(int j=1;j<=k;++j) S[i][j]=S[i-1][j-1]+S[i-1][j]*j; for(int i=1;i<=k;++i) ans+=(c*=modint(i).inv()*(n-i+1))*(tmp*=p)*(t2*=i)*S[k][i]; cout<<ans.val<<endl; return 0; } ``` 没线筛也过~(￣▽￣)~*，完全不卡常。 ```cpp #include<bits/stdc++.h> using namespace std; const int N=1e7+2; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} bool operator !=(const modint &o) {return val != o.val;} modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; } }frac[N],inv[N],ans,p,f[N],g[N],c,P,pw[N]; int n,m,k; inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1)res*=x;x*=x;} return res; } inline void init(int len){ frac[0]=inv[0]=1; for(int i=1;i <=len;++i) frac[i]=frac[i-1]*i; inv[len]=frac[len].inv();for(int i=len-1;i;--i) inv[i]=inv[i+1]*(i+1); if(n>=k){g[0]=1;for(int i=1;i<=len;++i) g[i]=g[i-1]*(frac[i-1]*inv[i]*(n+i-k))*p;} } inline modint C(const int n,const int m){return (n<m||m<0)?0:frac[n]*inv[m]*inv[n-m];} inline modint F(const int i){return g[i]*((i&1)?-1:1);} int main(){ cin>>n>>m>>k;p=modint(m).inv();init(10000000);f[0]=(n>=k);c=1;P=1; for(int i=1;i<=k;++i) f[i]=F(i)+p*F(i-1)-(p-1)*f[i-1]+(n+i==k); for(int i=1;i<=k;++i) ans+=expow(i,k)*(c*=inv[i]*frac[i-1]*(n-i+1))*f[k-i]*(P*=p); cout<<ans<<endl; return 0; } ``` ## [Lust](https://www.luogu.com.cn/problem/CF891E) 简单题，晚自习花十几分钟就切了。 发现每次答案的增加量刚好是 $\prod\limits_{i=1}^na_i$ 的减少量，于是设最终序列第 $i$ 个位置被减了 $b_i$ 次，则答案为 $\prod\limits_{i=1}^na_i-\prod\limits_{i=1}^n(a_i-b_i)$，于是问题转化到求后面部分。 考虑该题期望等于总和除以总方案数，而总方案数又刚好是 $n^k$，不难想到转化。 然后问题转化到对于所有 $\prod\limits_{i=1}^n{(a_i-b_i)}$ 求和，考虑以下两个性质： + 最后答案只与每个数被减了多少次有关。 + 减少数字的序列是有序的。 第一个性质让我们想到分别生成函数再乘起来，第二个性质让我们想到指数型生成函数。 设第 $i$ 个数的生成函数为 $G(i)$，则有： $$ \begin{aligned} G(i)&=\sum_{j=0}^\infty\frac{(a_i-j)}{j!}x^j\\ &=\sum_{j=0}^\infty (\frac{a_i}{j!}-\frac{1}{(j-1)!})x^j\\ &=e^x(a_i-x)\\ \therefore Res&=[x^k]e^{nx}\prod_{i=1}^n(a_i-x)\\ &=[x^k]e^{nx}\sum_{i=0}^nc_{n-i}x^i(-1)^i\\ &=k!\sum_{i=0}^nn^{k-i}c_{n-i}(-1)^i\frac{1}{(k-i)!}\\ \therefore Ans&=k!\sum_{i=0}^n\frac{c_{n-i}(-1)^i}{n^i(k-i)!} \end{aligned} $$ 其中 $c_i$ 是恰选 $i$ 个数，这 $i$ 个数的乘积之和，可使用 $O(n^2)$ 的递推解决，当然可以 $O(n\log^2n)$ 解决，不过就大材小用了。 ## [十二重计数法](https://www.luogu.com.cn/problem/P5824) 今天做了一遍，发现还挺简单的。 ~~代码有点长，找时间补。~~ 补了。 ### 盒子与球互相区分 1. 明显 $m^n$，入门难度。 2. 简单东西，相当于选 $m$ 个盒子使之有球，再区分球，${m\choose n}*n!$。 3. 用 EGF + 二项式定理推一下，可以发现一个盒子的生成函数是 $e^x-1$：$[x^n](e^x-1)^n=\sum\limits_{i=0}^m{m\choose i}(-1)^{m-i}i^n$。 ### 只有球相互区分 1. 第二类斯特林数，枚举有几个盒子有球：$\sum\limits_{i=0}^m{n\brace i}$。 2. 答案为 $[n\leq m]$。 3. 第二类斯特林数：${n\brace m}$。 ### 只有盒子互相区分 1. 允许为空的插板法：${n+m-1\choose m-1}$。 2. 选 $n$ 个有球的盒子：${m\choose n}$。 3. 先强制填满下界，然后当允许为空做：${n-1\choose m-1}$。 ### 都不互相区分 ~~唯一有难度的地方属于是。~~ #### 无其他限制 发现只和有多少个盒子选了多少个球有关。 考虑一下经典 DP 或者 Ferrers 图，可以发现答案是 $[x^n]\prod\limits_{i=1}^m\frac{1}{1-x^i}$。 就是付公主的背包，先求出 $\sum\limits_{i=1}^m\ln(1-x^i)$，再 exp + 求逆就可以解决。 可以考虑泰勒展开，但是今天想换一种推法： $$ \begin{aligned} G(x)&=\ln(1-x^i)\\ G(x)^{'}&=\frac{-ix^{i-1}}{1-x^i}\\ &=\sum_{j=0}^\infty x^{ij}*(-ix^{i-1})\\ &=-\sum_{j=1}^\infty ix^{ij-1}\\ G(x)&=-\sum_{j=1}^\infty\frac{x^{ij}}{j}\\ \end{aligned} $$ 在模意义下暴力做，复杂度是调和级数。 ### 至多一个球 还是 $[n\leq m]$。 ### 至少一个球 先填满下界，相当于无限制的第 $n-m$ 项，顺手做了。 ```cpp #include<bits/stdc++.h> using namespace std; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; } }; const modint g=3; const modint invg=(mod+1)/3; const int N=1e6+10+2e5; int rev[N],n,m; modint inv[N],frac[N],Inv[N]; void GetRev(int lim){ int bit=log2(lim); for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<(bit-1))); } inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } void Init(int len){ inv[1]=frac[0]=Inv[0]=1; for(int i=2;i<=len;++i) inv[i]=inv[mod%i]*(-mod/i); for(int i=1;i<=len;++i){frac[i]=frac[i-1]*i;Inv[i]=Inv[i-1]*inv[i];} } inline modint C(const int n,const int m){return ((n<m)||(n<0))?0:frac[n]*Inv[m]*Inv[n-m];} struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void derivation(const int len){for(int i=0;i<len;++i) a[i]=a[i+1]*(i+1);a[len-1]=0;} void integral(const int len){for(int i=len-1;i;--i) a[i]=a[i-1]*inv[i];a[0]=0;} void NTT(int limit,int type){ for(int i=0;i<limit;++i) if(i>rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; } } } if(type==-1){ modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } } }s1,s2,p; void getinv(Poly &y,Poly &x,const int len){ static Poly s,tmp; int lim=2,bas=1;s[0]=expow(x[0],mod-2); while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; GetRev(lim);s.NTT(lim,1);tmp.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=(s[i]*2)-s[i]*s[i]*tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; bas<<=1;lim<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=tmp[i]=0; } void ln(Poly &y,Poly &x,const int len){ static Poly t1,t2;int lim=1; for(int i=0;i<len;++i) t1[i]=x[i]; getinv(t2,t1,len);t1.derivation(len); for(;lim<(len<<1);lim<<=1) ;GetRev(lim); t1.NTT(lim,1);t2.NTT(lim,1);for(int i=0;i<lim;++i) t1[i]*=t2[i]; t1.NTT(lim,-1);t1.integral(len);for(int i=0;i<len;++i) y[i]=t1[i]; for(int i=0;i<=(len<<2);++i) t1[i]=t2[i]=0; } void exp(Poly &y,Poly &x,const int len){ static Poly s,tmp,t; int lim=2,bas=1;s[0]=1; while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; ln(t,s,bas);GetRev(lim);tmp[0]+=1; for(int i=0;i<bas;++i) tmp[i]-=t[i]; tmp.NTT(lim,1);s.NTT(lim,1); for(int i=0;i<lim;++i) s[i]*=tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; lim<<=1;bas<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=t[i]=tmp[i]=0; } inline int CalRev(const int len){int lim=1;for(;lim<=len;lim<<=1);GetRev(lim);return lim;} inline void Solve1(){cout<<expow(m,n)<<endl;} inline void Solve2(){cout<<C(m,n)*frac[n]<<endl;} inline void Solve3(){ modint res=0; for(int i=0;i<=m;++i) res+=C(m,i)*expow(i,n)*((m-i)&1?-1:1); cout<<res<<endl; } inline void Solve4(){ modint res=0; for(int i=1;i<=m;++i) res+=s1[i]; cout<<res<<endl; } inline void Solve5(){cout<<(n<=m)<<endl;} inline void Solve6(){cout<<s1[m]<<endl;} inline void Solve7(){cout<<C(n+m-1,m-1)<<endl;} inline void Solve8(){cout<<C(m,n)<<endl;} inline void Solve9(){cout<<C(n-1,m-1)<<endl;} inline void Solve10(){cout<<p[n]<<endl;} inline void Solve11(){cout<<(n<=m)<<endl;} inline void Solve12(){cout<<(n>=m?p[n-m]:0)<<endl;} signed main(){ cin>>n>>m;Init(400000);int lim; //第二类斯特林数 for(int j=0;j<=m;++j) s1[j]=Inv[j]*expow(j,n); for(int j=0;j<=m;++j) s2[j]=Inv[j]*(j&1?-1:1); lim=CalRev(m<<1);s1.NTT(lim,1);s2.NTT(lim,1); for(int i=0;i<lim;++i) s1[i]*=s2[i]; s1.NTT(lim,-1); //拆分数（付公主的背包） for(int i=1;i<=m;++i) for(int j=i;j<=n+5;j+=i) p[j]-=inv[j/i]; exp(p,p,n+5);getinv(p,p,n+5); Solve1();Solve2();Solve3(); Solve4();Solve5();Solve6(); Solve7();Solve8();Solve9(); Solve10();Solve11();Solve12(); return 0; } ``` ## [如何优雅的求和](https://www.luogu.com.cn/problem/P6667) 插值很难过去，考虑套路地将 $f(k)$ 转为下降幂多项式。 令 $f(x)=\sum\limits_{i=0}^ms_i{x\choose i}$，因为组合数本质上是一个下降幂的多项式，所以 $s$ 数列一定存在。 发现有 $f(k)=\sum\limits_{i=0}^k{k\choose i}s_i$，所以 $s_k=\sum\limits_{i=0}^k{k\choose i}(-1)^{k-i}f(i)$，一次 NTT 带走。 下面代入： $$ \begin{aligned} Q(f,n,x)&=\sum_{k=0}^nf(k){n\choose k}x^k(1-x)^{n-k}\\ &=\sum_{k=0}^n\sum_{i=0}^ks_i{k\choose i}{n\choose k}x^k(1-x)^{n-k}\\ &=\sum_{k=0}^n\sum_{i=0}^ks_i{n\choose i}{n-i\choose k-i}x^k(1-x)^{n-k}\\ &=\sum_{i=0}^ns_i{n\choose i}\sum_{k=i}^n{n-i\choose k-i}x^k(1-x)^{n-k}\\ &=\sum_{i=0}^ms_i{n\choose i}x^i\sum_{k=0}^{n-i}{n-i\choose k}x^k(1-x)^{n-k-i}\\ &=\sum_{i=0}^ms_i{n\choose i}x^i(x+1-x)^{n-i}\\ &=\sum_{i=0}^ms_i{n\choose i}x^i\\ \end{aligned} $$ 最后遍历一遍带走。 ```cpp #include<bits/stdc++.h> using namespace std; inline int read(){ int x=0,f=0; char c=getchar(); while(c<'0'||c>'9'){f|=(c=='-');c=getchar();} while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();} return f?-x:x; } const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } static constexpr int get_mod() { return mod; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; } }; const modint g=3; const int N=6e5+10; const modint invg=(mod+1)/3; modint frac[N],inv[N],ans,C[N]; int rev[N],n,d,m; inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } void GetRev(const int lim){ int bit=log2(lim)-1; for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<bit)); } int Rev(const int c){int lim=1;for(;lim<=c;lim<<=1);GetRev(lim);return lim;} struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void NTT(const int limit,const int type){ for(int i=0;i<limit;++i) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y;a[i+j+k]=x-y; } } } if(~type) return ; modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } void clear(int len){for(int i=0;i<=len;++i) a[i]=0;} void out(int len){for(int i=0;i<=len;++i) printf("%d ",a[i].val);printf("\n");} }s,f,G; void init(const int len){ frac[0]=1; for(int i=1;i<=len;++i) frac[i]=frac[i-1]*i;inv[len]=frac[len].inv(); for(int i=len-1;~i;--i) inv[i]=inv[i+1]*(i+1); } int main(){ init(200000); n=read();m=read();d=read();int lim=Rev(m<<1); for(int i=0;i<=m;++i) f[i]=read(); for(int i=0;i<=m;++i) f[i]*=inv[i]; for(int i=0;i<=m;++i) G[i]=inv[i]*((i&1)?-1:1); f.NTT(lim,1);G.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=f[i]*G[i]; s.NTT(lim,-1); for(int i=0;i<=m;++i) s[i]*=frac[i]; C[0]=1; for(int i=1;i<=m;++i) C[i]=C[i-1]*expow(i,mod-2)*(n-i+1); // for(int i=0;i<=m;++i) cout<<C[i]<<' ';cout<<endl; // for(int i=0;i<=m;++i) cout<<s[i]<<' ';cout<<endl; for(int i=0;i<=m;++i) ans+=s[i]*C[i]*expow(d,i); printf("%d\n",ans.val); return 0; } ``` ## [玩游戏](https://www.luogu.com.cn/problem/P4705) 有一个前置题目：[P7431](https://www.luogu.com.cn/problem/P7431)。 大部分不牵扯无穷的期望是假期望。 不难发现 $k$ 次的答案是 $\frac{\sum\limits_{i=1}^n\sum\limits_{j=1}^m(a_i+b_j)^k}{nm}$。 把式子的分子用二项式定理展开来看： $$ \begin{aligned} \sum\limits_{i=1}^n\sum\limits_{j=1}^m(a_i+b_j)^k&=\sum_{r=0}^k\sum_{i=1}^n\sum_{j=1}^m{k\choose r}a_i^rb_j^{k-r}\\ &=k!\sum_{r=0}^k(\frac1{r!}\sum_{i=1}^na_i^r)(\frac{1}{(k-r)!}\sum_{j=0}^mb_j^{k-r})\\ \end{aligned} $$ 现在问题变成了求 $\sum\limits_{i=1}^na_i^r$，因为求后面部分是同理的。 接下来就是我没想到的了。 $[x^r]\sum\limits_{i=1}^n\frac{1}{1-a_ix}=\sum\limits_{i=1}^na_i^r$，因为 $\frac{1}{1-a_ix}=\sum\limits_{j=0}^\infty a_i^jx^j$，所以将其累加起来就是我们想要的值。 处理加法可以 $\ln$ 再求导什么的，有种常数更大的做法（NTT+分治）。 具体地：$\frac{A(x)}{B(x)}+\frac{C(x)}{D(x)}=\frac{A(x)D(x)+B(x)C(x)}{B(x)D(x)}$，暴力通分后加起来，时间复杂度 $O(n\log^2n)$。 可以想一想为什么这样加封闭形式复杂度就对了。 常数：$\infty$。要注意很多地方的清空。 ```cpp #include<bits/stdc++.h> using namespace std; namespace Fread { const int SIZE=1<<21;char buf[SIZE],*S,*T; inline char getchar() {if(S==T){T=(S=buf)+fread(buf,1,SIZE,stdin);if(S==T)return '\n';}return *S++;} } namespace Fwrite { const int SIZE=1<<21; char buf[SIZE],*S=buf,*T=buf+SIZE; inline void flush(){fwrite(buf,1,S-buf,stdout);S=buf;} inline void putchar(char c){*S++=c;if(S==T)flush();} struct POPOSSIBLE{~POPOSSIBLE(){flush();}}ztr; } #define getchar Fread :: getchar #define putchar Fwrite :: putchar namespace Fastio{ struct Reader{ template<typename T> Reader& operator >> (T& x) { char c=getchar();x=0; while(c<'0'||c>'9') c=getchar(); while(c>='0'&&c<='9'){x=x*10+(c-'0');c=getchar();} return *this; } Reader(){} }cin; struct Writer{ template<typename T> Writer& operator << (T x) { if(x==0){putchar('0');return *this;} static int sta[45];int top=0; while(x){sta[++top]=x%10;x/=10;} while(top){putchar(sta[top]+'0');--top;} return *this; } Writer& operator << (char c) {putchar(c);return *this;} Writer(){} }cout; } #define endl '\n' #define cin Fastio :: cin #define cout Fastio :: cout const int g=3; const int mod=998244353; const int invg=(mod+1)/3; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } inline int Mod(const int x){return x>=mod?x-mod:x;} modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint& operator+=(const modint& o) { return val = Mod(val + o.val), *this; } modint& operator-=(const modint& o) { return val = norm(val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend Fastio :: Writer& operator<<(Fastio :: Writer& os, modint a) { return os << a.val; } friend Fastio :: Reader& operator>>(Fastio :: Reader& is, modint a) { return is >> a.val; } }; const int N=8e5+10; int rev[N],n,m,c[N],t,a[N],b[N]; modint frac[N],inv[N]; inline void GetRev(const int lim){int bit=log2(lim)-1;for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<bit));} inline int CalRev(const int len){int lim=1;for(;lim<=len;lim<<=1);GetRev(lim);return lim;} inline modint expow(modint x,int y){modint res=1;for(;y;y>>=1){if(y&1)res*=x;x*=x;}return res;} struct Poly{ modint a[N]; modint& operator [](const int x){return a[x];} inline void Readin(const int len){for(int i=0;i<len;++i) cin>>a[i].val;} inline void out(const int len){for(int i=0;i<len;++i) cout<<a[i]<<' ';cout<<endl;} inline void NTT(const int limit,const int type){ for(int i=0;i<limit;++i) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=a[i+j+k]*G; a[j+k]=x+y;a[i+j+k]=x-y; } } } if(~type) return ; modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } }x,y; struct Frac{ Poly s[2]; inline void NTT(int lim,int op){s[0].NTT(lim,op);s[1].NTT(lim,op);} }s[20],res; inline void getinv(Poly &y,Poly &x,const int len){ static Poly s,tmp; int lim=2,bas=1;s[0]=expow(x[0],mod-2); while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; GetRev(lim);tmp.NTT(lim,1);s.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=s[i]*2-s[i]*s[i]*tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; bas<<=1;lim<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<bas;++i) s[i]=tmp[i]=0; } void Solve(Frac &t,int dep,int l,int r){ if(l==r){ t.s[0][0]=1;t.s[1][0]=1; t.s[1][1]=-c[l]; return ; } int mid=(l+r)>>1,len=r-l+1; for(int i=0;i<=(len<<2);++i) t.s[0][i]=t.s[1][i]=0; Solve(s[dep],dep+1,l,mid); Solve(t,dep+1,mid+1,r); int lim=CalRev(len); s[dep].NTT(lim,1);t.NTT(lim,1); for(int i=0;i<lim;++i){ modint tmp=t.s[1][i]; t.s[0][i]=s[dep].s[0][i]*tmp+s[dep].s[1][i]*t.s[0][i]; t.s[1][i]=s[dep].s[1][i]*tmp; s[dep].s[0][i]=s[dep].s[1][i]=0; } t.NTT(lim,-1); } int main(){ cin>>n>>m;frac[0]=inv[0]=1;modint Inv=expow(1ll*n*m%mod,mod-2); for(int i=1;i<=n;++i) cin>>a[i]; for(int i=1;i<=m;++i) cin>>b[i]; cin>>t;int lim;for(int i=1;i<=t;++i) frac[i]=frac[i-1]*i; inv[t]=frac[t].inv();for(int i=t-1;i;--i) inv[i]=inv[i+1]*(i+1); for(int i=1;i<=n;++i) c[i]=a[i]; Solve(res,0,1,n);getinv(res.s[1],res.s[1],t+1);lim=CalRev((t+1)<<1); for(int i=t+1;i<=n;++i) res.s[0][i]=res.s[1][i]=0; res.NTT(lim,1);for(int i=0;i<lim;++i) res.s[0][i]*=res.s[1][i]; res.s[0].NTT(lim,-1);for(int i=0;i<=t;++i) x[i]=res.s[0][i]*inv[i]; for(int i=0;i<lim;++i) res.s[0][i]=res.s[1][i]=0; for(int i=1;i<=m;++i) c[i]=b[i]; Solve(res,0,1,m);getinv(res.s[1],res.s[1],t+1);lim=CalRev((t+1)<<1); for(int i=t+1;i<=m;++i) res.s[0][i]=res.s[1][i]=0; res.NTT(lim,1);for(int i=0;i<lim;++i) res.s[0][i]*=res.s[1][i]; res.s[0].NTT(lim,-1);for(int i=0;i<=t;++i) y[i]=res.s[0][i]*inv[i]; x.NTT(lim,1);y.NTT(lim,1); for(int i=0;i<lim;++i) x[i]*=y[i]; x.NTT(lim,-1);for(int i=1;i<=t;++i) cout<<(x[i]*frac[i]*Inv).val<<endl; return 0; } ``` ## [Partitions ](https://www.luogu.com.cn/problem/CF961G) 发现直接算是明显不可做的，有一个显然的式子：$Ans=\sum\limits_{i=1}^nw_i\sum\limits_{j=1}^n{n-1\choose j-1}j{n-j\brace k-1}$，对于每个数枚举其所在集合的大小与对应的方案，然后发现方案与前面的 $i$ 不挂钩，所以可以只化简后面的 $\sum\limits_{j=1}^n{n-1\choose j-1}j{n-j\brace k-1}$，但是这样也太麻烦了。 后面补这个式子的化简，其实有一个更简单的组合意义。 考虑元素 $j$ 加入元素 $i$ 所在的集合时，会使 $i$ 的贡献增加，也就是在划分好集合后，集合内每个点都对当前点造成 $i$ 的贡献。 所以式子：$Ans=\sum\limits_{i=1}^nw_i({n\brace k}+(n-1){n-1\brace k})$。 解释一下：自己对自己的贡献显然是 ${n\brace k}$，然后后面部分是指在集合划分好后，将当前点加入其中一个集合造成的一个贡献。 然后第二类斯特林数是有通项公式的，于是问题可以线性解决。 ```cpp #include<bits/stdc++.h> using namespace std; bool st; namespace Fread { const int SIZE=1<<21;char buf[SIZE],*S,*T; inline char getchar() {if(S==T){T=(S=buf)+fread(buf,1,SIZE,stdin);if(S==T)return '\n';}return *S++;} } namespace Fwrite { const int SIZE=1<<21; char buf[SIZE],*S=buf,*T=buf+SIZE; inline void flush(){fwrite(buf,1,S-buf,stdout);S=buf;} inline void putchar(char c){*S++=c;if(S==T)flush();} struct POPOSSIBLE{~POPOSSIBLE(){flush();}}ztr; } #define getchar Fread :: getchar #define putchar Fwrite :: putchar namespace Fastio{ struct Reader{ template<typename T> Reader& operator >> (T& x) { char c=getchar();T f=1; while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}x=0; while(c>='0'&&c<='9'){x=x*10+(c-'0');c=getchar();}x*=f; return *this; } Reader(){} }cin; struct Writer{ template<typename T> Writer& operator << (T x) { if(x==0){putchar('0');return *this;} if(x<0){putchar('-');x=-x;} static int sta[45];int top=0; while(x){sta[++top]=x%10;x/=10;} while(top){putchar(sta[top]+'0');--top;} return *this; } Writer& operator << (char c) {putchar(c);return *this;} Writer(){} }cout; } #define endl '\n' #define cin Fastio :: cin #define cout Fastio :: cout const int N=2e5+10; const int mod=1e9+7; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint operator-() const { return modint(norm(-val)); } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend Fastio :: Writer& operator<<(Fastio :: Writer& os, modint a) { return os << a.val; } }frac[N],inv[N]; int n,k;modint x,sum; void init(const int len=200000){ frac[0]=inv[0]=1; for(int i=1;i<=len;++i) frac[i]=frac[i-1]*i; inv[len]=frac[len].inv();for(int i=len-1;i;--i) inv[i]=inv[i+1]*(i+1); } inline modint expow(modint x,int y){modint res=1;for(;y;y>>=1){if(y&1)res*=x;x*=x;}return res;} inline modint Get(int n,int i){ modint res=0; for(int j=0;j<=i;++j) res+=expow(j,n)*inv[j]*inv[i-j]*((i-j)&1?-1:1); return res; } int main(){ init();cin>>n>>k; for(int i=1;i<=n;++i){cin>>x;sum+=x;} cout<<(sum*(Get(n,k)+(Get(n-1,k)*(n-1))))<<endl; return 0; } ``` ## Permutation 鉴于原题解对我来说语焉不详，这里写一发记录一下赛时思路。 有 $1,2,3,\cdots,n$ 这 $n$ 个数排成一排，你需要解决下面两个问题： 1. 你每次可以交换任意两个数，求**恰好** $k$ 次交换后能得到多少种不同的排列。 2. 你每次可以交换相邻两个数，求**恰好** $k$ 次交换后能得到多少种不同的排列。 由于答案可能很大，你只需要输出其对 $998244353$ 取模的结果即可。 $n\leq 10^9,k\leq 10^5$。 ### 交换任意两个位置 有几个经典结论： 1. 一个排列变成 $1,2,3,\cdots,n$ 的最小步数为 $n$ 减去轮换数。 2. 交换是可逆的，从原排列 $k$ 步能到排列 $P$ 等价于 $P$ 用 $k$ 步能到原排列。 3. 一个排列能恰好 $k$ 步到达当且仅当到其的最小步数小于等于 $k$，且奇偶性与 $k$ 相同。 证明：结论二三都是显然的，容易发现 $1,2,3,\cdots,n$ 的轮换数为 $n$，每次交换最多多出 $1$ 个轮换，故得证。 综上，用 $f_i$ 表示 ${n\brack n-i}$，则到时候对 $f$ 求一下和就是答案了。 问题在于如何求第一类斯特林数行的后 $k$ 项。 发现 $G(x)=\sum\limits_{i\geq 0}{n\brack i}x^i=x^{\overline n}$，基础常识。 所以 $G(x)=\prod\limits_{i=0}^{n-1}(x+i)$，发现我们要求的是后 $k$ 项，于是将其翻转得 $\prod\limits_{i=0}^{n-1}(1+ix)$，发现是[P5850](https://www.luogu.com.cn/problem/P5850) 于是做出来了。 关于反转的正确性：发现后 $k$ 项代表只有 $k$ 个地方不选 $x$，所以组合意义正确。 ### 交换两个相邻位置 首先要知道两个结论： 1. 一个排列变成 $1,2,3,\cdots,n$ 的最小步数为其逆序对数。 2. 一个排列能恰好 $k$ 步到达当且仅当到其的最小步数小于等于 $k$，且奇偶性与 $k$ 相同。 然后只需要会求：长为 $n$，逆序对数为 $k$ 的排列的个数就行了。 默认你做过 [P2513](https://www.luogu.com.cn/problem/P2513)，一道简单黄题。 发现 DP 式子是 $f_{i,j}=\sum\limits_{k=j-i+1}^j f_{i-1,k}=f_{i,j-1}+f_{i-1,j}-f_{i,j-i}$。 记录 $G_i(x)=\sum\limits_{j\geq 0}f_{i,j}x^j$，有 $G_i(x)=xG_i(x)+G_{i-1}(x)-x^iG_{i-1}(x)$。 所以有 $G_n(x)=\frac{1-x^n}{1-x}G_{n-1}(x)=\frac{1}{(1-x)^n}\prod\limits_{i=1}^n(1-x^i)$。 不难发现后面是 [P4389](https://www.luogu.com.cn/problem/P4389) 状物，$\ln+\exp$ 一套带走。 ```cpp #include<bits/stdc++.h> using namespace std; namespace Fread { const int SIZE=1<<21;char buf[SIZE],*S,*T; inline char getchar() {if(S==T){T=(S=buf)+fread(buf,1,SIZE,stdin);if(S==T)return '\n';}return *S++;} } namespace Fwrite { const int SIZE=1<<21; char buf[SIZE],*S=buf,*T=buf+SIZE; inline void flush(){fwrite(buf,1,S-buf,stdout);S=buf;} inline void putchar(char c){*S++=c;if(S==T)flush();} struct POPOSSIBLE{~POPOSSIBLE(){flush();}}ztr; } #define getchar Fread :: getchar #define putchar Fwrite :: putchar namespace Fastio{ struct Reader{ template<typename T> Reader& operator >> (T& x) { char c=getchar();T f=1; while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}x=0; while(c>='0'&&c<='9'){x=x*10+(c-'0');c=getchar();}x*=f; return *this; } Reader(){} }cin; struct Writer{ template<typename T> Writer& operator << (T x) { if(x==0){putchar('0');return *this;} if(x<0){putchar('-');x=-x;} static int sta[45];int top=0; while(x){sta[++top]=x%10;x/=10;} while(top){putchar(sta[top]+'0');--top;} return *this; } Writer& operator << (char c) {putchar(c);return *this;} Writer(){} }cout; } #define endl '\n' #define cin Fastio :: cin #define cout Fastio :: cout const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint a) { return os << a.val; } friend Fastio :: Writer& operator<<(Fastio :: Writer& os, modint a) { return os << a.val; } friend Fastio :: Reader& operator>>(Fastio :: Reader& is, modint a) { return is >> a.val; } }; const modint g=3; const modint invg=(mod+1)/3; const int N=10+6e5; int rev[N],n,m,k; modint inv[N],frac[N],Inv[N]; void GetRev(int lim){ int bit=log2(lim); for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<(bit-1))); } inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } void Init(int len){ inv[1]=frac[0]=Inv[0]=1; for(int i=2;i<=len;++i) inv[i]=inv[mod%i]*(-mod/i); for(int i=1;i<=len;++i){frac[i]=frac[i-1]*i;Inv[i]=Inv[i-1]*inv[i];} } inline modint C(const int n,const int m){return ((n<m)||(n<0))?0:frac[n]*Inv[m]*Inv[n-m];} struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void derivation(const int len){for(int i=0;i<len;++i) a[i]=a[i+1]*(i+1);a[len-1]=0;} void integral(const int len){for(int i=len-1;i;--i) a[i]=a[i-1]*inv[i];a[0]=0;} inline void clear(const int len){for(int i=0;i<len;++i) a[i]=modint();} void NTT(int limit,int type){ for(int i=0;i<limit;++i) if(i>rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; } } } if(type==-1){ modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } } }s1,s2,p,x,y,G; void getinv(Poly &y,Poly &x,const int len){ static Poly s,tmp; int lim=2,bas=1;s[0]=expow(x[0],mod-2); while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; GetRev(lim);s.NTT(lim,1);tmp.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=(s[i]*2)-s[i]*s[i]*tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; bas<<=1;lim<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=tmp[i]=0; } void ln(Poly &y,Poly &x,const int len){ static Poly t1,t2;int lim=1; for(int i=0;i<len;++i) t1[i]=x[i]; getinv(t2,t1,len);t1.derivation(len); for(;lim<(len<<1);lim<<=1) ; GetRev(lim);t1.NTT(lim,1);t2.NTT(lim,1); for(int i=0;i<lim;++i) t1[i]*=t2[i]; t1.NTT(lim,-1);t1.integral(len);for(int i=0;i<len;++i) y[i]=t1[i]; for(int i=0;i<=(len<<2);++i) t1[i]=t2[i]=0; } void exp(Poly &y,Poly &x,const int len){ static Poly s,tmp,t; int lim=2,bas=1;s[0]=1; while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; ln(t,s,bas);GetRev(lim);tmp[0]+=1; for(int i=0;i<bas;++i) tmp[i]-=t[i]; tmp.NTT(lim,1);s.NTT(lim,1); for(int i=0;i<lim;++i) s[i]*=tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; lim<<=1;bas<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<=(len<<2);++i) s[i]=t[i]=tmp[i]=0; } inline int CalRev(const int len){int lim=1;for(;lim<=len;lim<<=1);GetRev(lim);return lim;} void init(int k,int m){ for(int i=0;i<=m;++i) x[i]=Inv[i+1]; for(int i=0;i<=m;++i) y[i]=expow(k+1,i+1)*Inv[i+1]; getinv(x,x,m+1);int lim=CalRev((m+1)<<1); x.NTT(lim,1);y.NTT(lim,1); for(int i=0;i<lim;++i) G[i]=x[i]*y[i]; G.NTT(lim,-1);for(int i=0;i<lim;++i) G[i]*=frac[i]; G[0]=0;x.clear(lim);y.clear(lim); for(int i=1;i<=m;++i) G[i]*=inv[i]*((i&1)?1:-1); exp(G,G,m+1);G[1]=1ll*k*(k+1)/2;G[0]=1; } signed main(){ cin>>n>>k;Init(400000);init(n-1,k); for(int i=2;i<=k;++i) G[i]+=G[i-2]; cout<<G[k]<<endl;int lim; for(int i=1;i<=min(n,k);++i) for(int j=i;j<=k;j+=i) p[j]-=inv[j/i]; exp(p,p,k+1);s1[0]=1;for(int i=1;i<=k;++i) s1[i]=s1[i-1]*(n+i-1)*inv[i]; lim=CalRev((k+1)<<1);s1.NTT(lim,1);p.NTT(lim,1); for(int i=0;i<lim;++i) p[i]*=s1[i]; p.NTT(lim,-1);for(int i=2;i<=k;++i) p[i]+=p[i-2]; cout<<p[k]<<endl; return 0; } ``` ## [Bags with Balls](https://www.luogu.com.cn/problem/CF1716F) 啥也不会就只好来水水数学了。 考虑直接列式子爆算。 记录 $x=\lceil \frac m2 \rceil$，也就是一个袋子中标号为奇数的数量，$y=m-x$。 考虑枚举有多少个标号是奇数的球 $Ans=\sum\limits_{i=1}^nx^iy^{n-i}{n\choose i}i^k$。 套路地将 $i^k$ 展开，并使用经典套路化简： $$ \begin{aligned} Ans&=\sum_{i=1}^nx^iy^{n-i}{n \choose i}\sum_{j=1}^k{i\choose j}j!{k\brace j}\\ &=\sum_{j=1}^k {k\brace j}j!\sum_{i=1}^n \frac{n!}{j!(i-j)!(n-i)!}x^iy^{n-i}\\ &=\sum_{j=1}^k{k\brace j}j!{n\choose j}\sum_{i=1}^n{n-j\choose i-j}x^iy^{n-i}\\ &=\sum_{j=1}^k{k\brace j}j!{n\choose j}x^j\sum_{i=0}^{n-j}{n-j\choose i}x^iy^{n-i-j}\\ &=\sum_{j=1}^k{k\brace j}n^{\underline j}x^jm^{n-j}\\ \end{aligned} $$ 化简完毕，$O(k^2)$ 预处理斯特林数，计算时线性递推组合数和幂次，$O(k^2+Tk)$，可以过了。 考虑进一步优化呢？使用通项公式再次将斯特林数展开………也就是用 Cards 加强版的 Trick。反正能过，不推。 ```cpp #include<bits/stdc++.h> using namespace std; namespace Fread { const int SIZE=1<<21;char buf[SIZE],*S,*T; inline char getchar() {if(S==T){T=(S=buf)+fread(buf,1,SIZE,stdin);if(S==T)return '\n';}return *S++;} } namespace Fwrite { const int SIZE=1<<21; char buf[SIZE],*S=buf,*T=buf+SIZE; inline void flush(){fwrite(buf,1,S-buf,stdout);S=buf;} inline void putchar(char c){*S++=c;if(S==T)flush();} struct POPOSSIBLE{~POPOSSIBLE(){flush();}}ztr; } #define getchar Fread :: getchar #define putchar Fwrite :: putchar namespace Fastio{ struct Reader{ template<typename T> Reader& operator >> (T& x) { char c=getchar();T f=1; while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}x=0; while(c>='0'&&c<='9'){x=x*10+(c-'0');c=getchar();}x*=f; return *this; } Reader(){} }cin; struct Writer{ template<typename T> Writer& operator << (T x) { if(x==0){putchar('0');return *this;} if(x<0){putchar('-');x=-x;} static int sta[45];int top=0; while(x){sta[++top]=x%10;x/=10;} while(top){putchar(sta[top]+'0');--top;} return *this; } Writer& operator << (char c) {putchar(c);return *this;} Writer(){} }cout; } #define endl '\n' #define cin Fastio :: cin #define cout Fastio :: cout #define int unsigned long long const int mod=998244353; struct modint { int val; modint() : val(0) {} modint(const int& m) : val(m) {} modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend Fastio :: Writer& operator<<(Fastio :: Writer& os, modint a) { return os << a.val; } }S[2010][2010],ans,inv; int n,m,T,k,x; inline modint expow(modint x,int y){modint res=1;for(;y;y>>=1){if(y&1)res*=x;x*=x;}return res;} signed main(){ cin>>T; S[0][0]=1;for(int i=1;i<=2000;++i) for(int j=1;j<=2000;++j) S[i][j]=S[i-1][j-1]+S[i-1][j]*j; while(T--){ cin>>n>>m>>k;x=(m+1)>>1; ans=0;inv=expow(m,mod-2); modint tmp=1,tmp2=1,tmp3=expow(m,n); for(int j=1;j<=k;++j) ans+=(tmp*=(n-j+1))*S[k][j]*(tmp2*=x)*(tmp3*=inv); cout<<ans<<endl; } return 0; } ``` ## [重返现世](https://www.luogu.com.cn/problem/P4707) 做的第一道 $k\min\max$ 容斥，但是最开始没看到 $m$ 的限制。 $E(\max\limits_k(S))=\sum\limits_{T\subset S}(-1)^{|T|-k}{|T|-1\choose k-1}E(\min(T))$。 考虑求 $E(\min(T))=\frac {m}{\sum\limits_{x\in T}p_x}$，理由显然，进行 DP。 状态为 $f_{i,j,k}$ DP 到第 $i$ 个数，$\sum p=j$，求第 $k$ 大的 $(-1)^{|T|-k}{|T|-1\choose k-1}$ 之和，先不考虑初值怎么设（有些麻烦）。 首先可以不选这个数，$f_{i,j,k}+=f_{i-1,j,k}$。 然后选，此时 $|T|$ 加 $1$，考虑对每一项计算影响。 对于 $(-1)^{|T|-k}$ 这一项，发现刚好多了一个 $-1$，乘上，然后发现 ${|T|-1\choose k-1}={|T|-2 \choose k-1}+{|T|-2\choose k-1}$。 此时我们存 $k$ 这一项的用处就显现出来了，$f_{i,j,k}+=f_{i-1,j-p_i,k-1}-f_{i-1,j-p_i,k}$，注意是减号。 发现是背包，滚动一下，最后大力计算就行了。 发现初值并不好设，对于不为 $0$ 的 $k$，令 $f_{0,0,k}=-1$ 就行了。 考场想不出来的话可以解个方程（因为你大概率只有初值不知道，逆推就行了）。 ```cpp #include<bits/stdc++.h> using namespace std; namespace Fread { const int SIZE=1<<21;char buf[SIZE],*S,*T; inline char getchar() {if(S==T){T=(S=buf)+fread(buf,1,SIZE,stdin);if(S==T)return '\n';}return *S++;} } namespace Fwrite { const int SIZE=1<<21; char buf[SIZE],*S=buf,*T=buf+SIZE; inline void flush(){fwrite(buf,1,S-buf,stdout);S=buf;} inline void putchar(char c){*S++=c;if(S==T)flush();} struct POPOSSIBLE{~POPOSSIBLE(){flush();}}ztr; } #define getchar Fread :: getchar #define putchar Fwrite :: putchar namespace Fastio{ struct Reader{ template<typename T> Reader& operator >> (T& x) { char c=getchar();T f=1; while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}x=0; while(c>='0'&&c<='9'){x=x*10+(c-'0');c=getchar();}x*=f; return *this; } Reader(){} }cin; struct Writer{ template<typename T> Writer& operator << (T x) { if(x==0){putchar('0');return *this;} if(x<0){putchar('-');x=-x;} static int sta[45];int top=0; while(x){sta[++top]=x%10;x/=10;} while(top){putchar(sta[top]+'0');--top;} return *this; } Writer& operator << (char c) {putchar(c);return *this;} Writer(){} }cout; } #define endl '\n' #define cin Fastio :: cin #define cout Fastio :: cout const int N=1e4+10; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend Fastio :: Writer& operator<<(Fastio :: Writer& os, modint a) { return os << a.val; } }f[N][12]; int n,k,m,p[N]; int main(){ cin>>n>>k>>m;modint ans=0;k=n-k+1; for(int i=1;i<=n;++i) cin>>p[i]; for(int i=1;i<=k;++i) f[0][i]=-1; for(int i=1;i<=n;++i) for(int j=m;j>=p[i];--j) for(int c=k;c;--c) f[j][c]+=f[j-p[i]][c-1]-f[j-p[i]][c]; for(int i=1;i<=m;++i) ans+=f[i][k]*m*modint(i).inv(); cout<<ans<<endl; return 0; } ``` ## [[HEOI2016/TJOI2016] 求和](https://www.luogu.com.cn/problem/P4091) 太简单了。 $$ \begin{aligned} f(n)&=\sum_{i=0}^n\sum_{j=0}^n{i\brace j}2^j j!\\ &=\sum_{i=0}^n\sum_{j=0}^n2^j\sum_{k=0}^j{j\choose k}(-1)^{j-k}k^i\\ &=\sum_{j=0}^n2^j\sum_{k=0}^j{j\choose k}(-1)^{j-k}[n+1]_k\\ \end{aligned} $$ 可以卷积。 ```cpp #include<bits/stdc++.h> using namespace std; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint& a) { return os << a.val; } }; const modint g=3; const int N=6e5+10; const modint invg=(mod+1)/3; modint frac[N],inv[N],ans,C[N]; int rev[N],n,d,m; inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } void GetRev(const int lim){ int bit=log2(lim)-1; for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<bit)); } int Rev(const int c){int lim=1;for(;lim<=c;lim<<=1);GetRev(lim);return lim;} struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} void NTT(const int limit,const int type){ for(int i=0;i<limit;++i) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y;a[i+j+k]=x-y; } } } if(~type) return ; modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } void clear(int len){for(int i=0;i<=len;++i) a[i]=0;} void out(int len){for(int i=0;i<=len;++i) printf("%d ",a[i].val);printf("\n");} }x,y; void init(const int len){ frac[0]=1; for(int i=1;i<=len;++i) frac[i]=frac[i-1]*i; inv[len]=frac[len].inv(); for(int i=len-1;~i;--i) inv[i]=inv[i+1]*(i+1); } int main(){ ios::sync_with_stdio(0); cin.tie(0);cout.tie(0);init(200000); cin>>n;x[0]=1;x[1]=n+1;int lim=1;for(;lim<=(n<<1);lim<<=1); for(int i=2;i<=n;++i) x[i]=(expow(i,n+1)-1)*(modint(i-1).inv())*inv[i]; for(int i=0;i<=n;++i) y[i]=inv[i]*((i&1)?-1:1); GetRev(lim);x.NTT(lim,1);y.NTT(lim,1); for(int i=0;i<lim;++i) x[i]*=y[i]; x.NTT(lim,-1);modint ans=0; for(int i=0;i<=n;++i) ans+=expow(2,i)*x[i]*frac[i]; cout<<ans<<endl; return 0; } ``` ## [猎人杀](https://www.luogu.com.cn/problem/P5644) [交了题解的](https://www.luogu.com.cn/article/az6d62ok) 稍微有点经验就能看出来概率一直变化是假的。 这题与百鸽笼的差距就是白鸽笼求得是期望次数，所以要实时处理，这题不用，所以可以转化一下。 每一个猎人死亡后仍然可能被击中，问 $1$ 号最后死的概率。 然后考虑容斥原理，不失套路地设 $W(S)=\sum\limits_{x\in S}w_x$，$P(S)$ 为在 $1$ 之后死的至少包括 $S$ 的概率。 $Ans=\sum\limits_{S}(-1)^{|S|}P(S)$，也就是常规的偶数减奇数。 然后考虑 $P(S)$ 怎么求，发现 $S$ 肯定一直不能死，$1$ 肯定是最后死，我们枚举 $1$ 什么时候死。 $P(S)=\sum\limits_{i=0}^\infty (\frac{tot-a_1-W(S)}{tot})^i\frac{a_1}{tot}=\frac {a_1}{tot}\frac{tot}{a_1+W(S)}=\frac{a_1}{a_1+W(S)}$。 所以 $Ans=\sum\limits_{S}(-1)^{|S|}\frac{a_1}{a_1+W(S)}$。 发现现在只跟元素个数和权值和有关了，而权值和很小，不难想到可以计算 $W(S)$ 为某个权值的 $(-1)^{|S|}$ 之和。 即使用生成函数计算 $\prod\limits_{i=2}^n(1-x^{a_i})$ 可以做到 $O(n\log n)$，不过本题数据范围下放过了 $O(n\log^2 n)$。 ```cpp #include<bits/stdc++.h> using namespace std; namespace Fread { const int SIZE=1<<21;char buf[SIZE],*S,*T; inline char getchar() {if(S==T){T=(S=buf)+fread(buf,1,SIZE,stdin);if(S==T)return '\n';}return *S++;} } namespace Fwrite { const int SIZE=1<<21; char buf[SIZE],*S=buf,*T=buf+SIZE; inline void flush(){fwrite(buf,1,S-buf,stdout);S=buf;} inline void putchar(char c){*S++=c;if(S==T)flush();} struct POPOSSIBLE{~POPOSSIBLE(){flush();}}ztr; } #define getchar Fread :: getchar #define putchar Fwrite :: putchar namespace Fastio{ struct Reader{ template<typename T> Reader& operator >> (T& x) { char c=getchar();T f=1; while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}x=0; while(c>='0'&&c<='9'){x=x*10+(c-'0');c=getchar();}x*=f; return *this; } Reader(){} }cin; struct Writer{ template<typename T> Writer& operator << (T x) { if(x==0){putchar('0');return *this;} if(x<0){putchar('-');x=-x;} static int sta[45];int top=0; while(x){sta[++top]=x%10;x/=10;} while(top){putchar(sta[top]+'0');--top;} return *this; } Writer& operator << (char c) {putchar(c);return *this;} Writer& operator << (const char* str){int cur=0;while(str[cur])putchar(str[cur++]);return *this;} Writer(){} }cout; } #define endl '\n' #define cin Fastio :: cin #define cout Fastio :: cout const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } }; const modint g=3; const modint invg=(mod+1)/3; const int N=1e6+10; int rev[N],n,A,B,sum[N],w[N]; modint frac[N],inv[N]; void GetRev(int lim){ int bit=log2(lim); for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<(bit-1))); } inline modint expow(modint x,int y){ modint res=1; for(;y;y>>=1){if(y&1) res*=x;x*=x;} return res; } struct Poly{ modint a[N]; modint& operator [](const int &x){return a[x];} inline void Readin(const int len){for(int i=0;i<len;++i) cin>>a[i].val;} inline void out(const int len){for(int i=0;i<len;++i) cout<<a[i].val<<' ';putchar('\n');} inline void clear(const int len){for(int i=0;i<len;++i) a[i]=modint();} void NTT(int limit,int type){ for(int i=0;i<limit;++i) if(i>rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=G*a[i+j+k]; a[j+k]=x+y;a[i+j+k]=x-y; } } } if(type==-1){ modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } } }f[18]; void mul(int dep,int l,int r){ int len=sum[r]-sum[l-1]; int mid=(l+r)>>1;f[dep].clear(len+5); if(l==r){f[dep][0]=1;f[dep][w[l]]=-1;return ;} mul(dep+1,l,mid);int len1=sum[mid]-sum[l-1]; for(int i=0;i<=len1;++i) f[dep][i]=f[dep+1][i]; f[dep+1].clear(len1+5);mul(dep+1,mid+1,r); int lim=1;while(lim<=len) lim<<=1; GetRev(lim);f[dep+1].NTT(lim,1);f[dep].NTT(lim,1); for(int i=0;i<lim;++i) f[dep][i]*=f[dep+1][i]; f[dep].NTT(lim,-1);f[dep+1].clear(lim); } signed main(){ cin>>n; for(int i=1;i<=n;++i) cin>>w[i]; for(int i=1;i<=n;++i) sum[i]=sum[i-1]+w[i]; mul(0,2,n);modint ans=0; // for(int i=0;i<=sum[n];++i) cout<<f[0][i].val<<' '; // cout<<endl; for(int i=0;i<=sum[n];++i) ans+=((modint(w[1]+i)).inv())*f[0][i]*w[1]; cout<<ans.val<<endl; return 0; } ``` ## 【BZOJ5093】图的价值 简单无向图是指无重边、无自环的无向图（不一定连通）。 一个带标号的图的价值定义为每个点度数的k次方的和。 给定 $n$ 和 $k$，请计算所有 $n$ 个点的带标号的简单无向图的价值之和。 因为答案很大，请对 $998244353$ 取模输出。 不难想到对于每个点算贡献，最后再乘上 $n$。 每个点贡献为： $$ \begin{aligned} n\sum_{i=0}^{n-1}{n-1\choose i}2^{{n-1\choose 2}}i^k&=n\times 2^{{n-1\choose 2}}\sum_{i=0}^{n-1}{n-1\choose i}i^k \end{aligned} $$ 只看求和号里的部分，将 $i^k$ 使用第二类斯特林数展开，有： $$ \begin{aligned} \sum_{i=0}^{n-1}{n-1\choose i}i^k&=\sum_{i=0}^{n-1}{n-1\choose i}\sum_{j=0}^i{k\brace j}{i \choose j}j!\\ &=\sum_{j=0}^kj!{k\brace j}\sum_{i=0}^{n-1}{n-1\choose i}{i\choose j}\\ &=\sum_{j=0}^kj!{k\brace j}{n-1\choose j}\sum_{i=0}^{n-1}{n-j-1\choose i-j}\\ &=\sum_{j=0}^kj!{k\brace j}{n-1\choose j}2^{n-j-1}\\ &=\sum_{j=0}^k{k\brace j}(n-1)^{\underline{j}}2^{n-j-1} \end{aligned} $$ 考虑第二类斯特林数的通项公式： $$ \begin{aligned} i^k&=\sum_{j=0}^i {k\brace j}{i\choose j}j!\\ {k\brace i}&=\sum_{j=0}^i\frac{(-1)^{i-j}j^k}{j!(i-j)!} \end{aligned} $$ ```cpp #include<bits/stdc++.h> using namespace std; const int g=3; const int N=6e5+10; const int mod=998244353; const int invg=(mod+1)/3; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint a) { return os << a.val; } }frac[N],inv[N],f[N],ans; int rev[N<<2],n,k,m; inline void GetRev(int lim){int bit=__lg(lim)-1;for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<bit));} inline int CalRev(int len){int lim=1;while(lim<=(len))lim<<=1;GetRev(lim);return lim;} inline modint expow(modint x,int y){modint res=1;for(;y;y>>=1){if(y&1)res*=x;x*=x;}return res;} struct Poly{ modint a[N<<2]; modint& operator [](const int x){return a[x];} inline void NTT(int limit,int typ){ for(int i=0;i<limit;++i) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow((typ==1)?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=a[i+j+k]*G; a[j+k]=x+y;a[i+j+k]=x-y; } } } if(~typ) return ; modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } }x,y; inline void init(const int len=200000){//f_i:(n-1) 的 i 次下降幂 frac[0]=inv[0]=1;for(int i=1;i<=len;++i) frac[i]=frac[i-1]*i; inv[len]=frac[len].inv();for(int i=len-1;i;--i) inv[i]=inv[i+1]*(i+1); f[0]=1;for(int i=1;i<=k;++i) f[i]=f[i-1]*(n-i); } inline void Stirling(const int k){//k:行号 for(int i=0;i<=k;++i) x[i]=inv[i]*expow(i,k); for(int i=0;i<=k;++i) y[i]=inv[i]*((i&1)?-1:1); int lim=CalRev((k<<1));x.NTT(lim,1);y.NTT(lim,1); for(int i=0;i<lim;++i) x[i]*=y[i]; x.NTT(lim,-1); } inline modint C(int n,int m){return (m<0||n<m)?0:frac[n]*inv[m]*inv[n-m];} int main(){ ios::sync_with_stdio(0);cin.tie(0);cout.tie(0); cin>>n>>k;init();Stirling(k); for(int i=0;i<=k;++i) ans+=x[i]*f[i]*expow(2,n-i-1); ans*=expow(2,(1ll*(n-2)*(n-1)/2)%(mod-1))*n; cout<<ans<<endl; return 0; } ``` ## [ARC154F](https://www.luogu.com.cn/problem/AT_arc154_f) 也是奇怪概论题。 考虑使用比较 Bot 的做法，记 $f_{i,j}$ 表示抛了 $j$ 次抛出了 $i$ 种数的概率，有 $f_{i,j}=\frac inf_{i,j-1}+\frac{n-i+1}nf_{i-1,j-1}$，记 $F_i(x)=\sum\limits_{j=0}^\infty f_{i,j}x^j$。 有 $F_i(x)=x(\frac inF_i(x)+\frac{n-i+1}nF_{i-1}(x))$，得到 $F_i(x)=\frac{x(n-i+1)F_{i-1}(x)}{n-ix}=\frac{x^i n^{\underline{i}}}{\prod\limits_{j=1}^{i}(n-jx)}$。 有 $F_n(x)=\frac xnF_{n-1}(x)=\frac{(n-1)!x^n}{\prod\limits_{j=1}^{n-1}(n-ix)}$。 这个可以分治的过程中 NTT 一下，即可求出 $F_n(x)$，$F_n(e^x)$ 就是答案的 EGF。 思考为什么 $F_n(e^x)$ 是答案的 EGF，将之展开： $$ \begin{aligned} F_n(e^x)&=\sum_{i=0}^\infty f_{n,i}e^{ix}\\ &=\sum_{i=0}^\infty f_{n,i}\sum_{j=0}^\infty \frac{(ix)^j}{j!}\\ &=\sum_{j=0}^\infty \frac {x^j}{j!}\sum_{i=0}^\infty i^jf_{n,i}\\ \end{aligned} $$ 明显 $[\frac{x^k}{k!}]F_n(e^x)$ 是次数的 $k$ 方和的期望值。 然后考虑求 $F_n(e^x)$，是套路的，再来一次 NTT 分治，需要会 EGF 转 OGF。 具体地，记 $G(x)=\prod\limits_{j=1}^{n-1}(n-ix)$，那么 $F(e^x)=\frac{(n-1)!e^{nx}}{G(e^x)}$。 考虑在知晓 $G(x)$ 的情况下如何求 $G(e^x)$。 可以使用 EGF 转 OGF 的 trick。 再 NTT 分治一下就行了。 NTT 分治合并分子分母有一个[经典题](https://www.luogu.com.cn/problem/P7431)。 代码咕咕咕，有时间再写！ ## [ARC156D](https://www.luogu.com.cn/problem/AT_arc156_d) 感觉简单于 C。 记 $c_i$ 为 $\sum\limits_{j=1}^K [p_j=i]$，那么相当于求 $\sum\limits_{i=1}^na_ic_i$ 的异或和。 对于指定的一组 $c$，产生其的方案有 ${K\choose c_1,c_2,\dots,c_n}$ 种，也就是多重集排列。 发现如果其方案为偶数，则这组 $c$ 没有贡献，考虑找到方案为奇数的充要条件。 ${K\choose c_1,c_2,\dots,c_n}=\frac{K!}{c_1!c_2!\dots c_n!}$，也就是说 $c_i!$ 中 $2$ 因子的个数和要刚好等于 $K!$ 的 $2$ 因子个数。 用 $f(x)$ 表示 $x!$ 中 $2$ 因子的个数，不难发现 $f(x)=\lfloor \frac x2\rfloor+f(\lfloor \frac x2\rfloor)$，手玩一下可以看出 $f(x)=x-\operatorname{bit(x)}$，$\operatorname{bit}(x)$ 表示 $x$ 二进制下 $1$ 的个数。 所以 ${K\choose c_1,c_2,\dots,c_n}$ 为奇数的一个必要条件是 $\operatorname{bit}(K)=\sum\limits_{i=1}^n \operatorname{bit}(c_i)$，又因为 $K=\sum\limits_{i=1}^n c_i$。 所以有贡献的 $c_i$ 一定是对于 $K$ 的二进制分解，即将 $K$ 二进制下每一位都分配给一个 $c_i$（注意两个二进制位分给了同一个 $c_i$ 也是合法的）。 现在考虑 DP。 使用 $dp(i,j)$ 表示正准备使用第 $i$ 位，当前和为除以 $2^{i-1}$ 后为 $j$ 的方案的奇偶性，稍微 DP 一下就行了。 ```cpp #include<bits/stdc++.h> #define int long long #define endl '\n' using namespace std; const int B=61; const int N=1<<13; const int Max=1<<11; int f[B][N],a[N],n,k,ans; signed main(){ ios::sync_with_stdio(0);cin.tie(0);cout.tie(0); cin>>n>>k;for(int i=1;i<=n;++i) cin>>a[i]; if(k&1){for(int i=1;i<=n;++i) f[0][a[i]]^=1;} else f[0][0]=1; for(int i=1;i<B;++i){ if((k>>i)&1){for(int j=0;j<=Max;++j) for(int k=1;k<=n;++k) f[i][(j>>1)+a[k]]^=f[i-1][j];} else {for(int j=0;j<=Max;++j) f[i][j>>1]^=f[i-1][j];} } for(int i=0;i<B;++i){ if((n&1)||!(k>>i)){ int tmp=0; for(int j=0;j<=Max;++j) if(j&1) tmp^=f[i][j]; if(tmp) ans+=1ll<<i; } } cout<<ans<<endl; return 0; } ``` ## [ARC156E](https://www.luogu.com.cn/problem/AT_arc156_e) 考虑允许 $(i,i+1)$ 的边时如何判定是否合法。 首先显然地 $S=\sum\limits_{i=1}^n X_i$ 得是偶数，因为一条无向边可以使两个点的度数加一，如果度数为奇数显然不合法。 然后发现 $\max\limits_{i=1}^n(X_i)\le\frac S2$，否则只能有自环了。 那么不允许 $(i,i+1)$ 时怎么做呢？可以类似地发现 $\max(X_i,X_{i+1})\le\frac S2$，否则要么有自环，要么 $(i,i+1)$ 有边。 发现直接计数好像不是很方便，考虑容斥。 记总方案数为 $A$，存在一个 $i$ 使得 $X_i+X_{i-1}> \frac S2$ 的方案数为 $B$，存在一个 $i$ 使得 $X_i+X_{i-1}>\frac S2,X_{i+1}+X_i>\frac S2$ 的方案数为 $C$。 则答案为 $A-B+C$。 #### 考虑求总方案数 $A$： 每个位置的生成函数是 $F(x)=\sum\limits_{i=0}^m x^i=\frac{1-x^{m+1}}{1-x}$，所以序列的生成函数是 $G(x)=F^n(x)=\frac{(1-x^{m+1})^n}{(1-x)^n}$，分成上下两部分。 1. $(1-x^{m+1})^n=\sum\limits_{i=0}^n(-1)^{n-i}x^{i(m+1)}$。 2. $\frac 1{(1-x)^n}=\sum\limits_{i=0}^\infty {n+i-1 \choose i}x^i$。 要求生成函数偶次项系数之和，可以处理下面部分的前缀和，枚举上面部分。容易的。 #### 求至少有一个地方不满足条件的方案数 $B$： 以下记 $\operatorname{Calc}(x,y)$ 表示表示将 $x$ 个无标号小球放进 $y$ 个有标号盒子中的方案数。 类似计算总方案数 $A$ 的方法计算，不难发现是 $O(\frac xm)$ 的，因为生成函数上面部分每 $m$ 位有一个值。 假设这两个位置是 $X_1,X_2$，到时候把方案乘上 $n$ 即可。 枚举 $S$ 和 $T=X_1+X_2$，不难发现此时 $S\le 4\times m$，所以枚举是 $O(m^2)$ 的。 然后求 $X_{3\dots n}$ 的分配方案是容易的，就是 $\operatorname{Calc}(S-T,n-2)$，此时 $\operatorname{Calc}$ 是常数时间的。 #### 求两个地方不满足条件的方案数 $C$： 枚举 $T=X_1+X_2+X_3$ 和 $S$。 发现随着 $X_2$ 的增加，方案数是等差数列，直接做了。 代码也是咕咕咕！ ## [ARC160D](https://www.luogu.com.cn/problem/AT_arc160_d) 有点招笑但很经典的套路题，每一步都没有跳出套路。 套路地考虑将操作序列和权值序列构建双射。 思考为什么会不构成双射，一定是某一些权值序列有多种方法被构造出来。 我们考虑尽可能使用一操作去构造序列，换句话说，如果一个序列有多种方式可以构造出来，我们选择一操作最多的那种，明显就构成双射了。 具体来讲呢？发现每一个元素被操作二所选定为左端点的次数不会超过 $k-1$ 次，否则一定可以使用 $k$ 次操作一来代替 $k$ 次操作二。 所以解法就出来了，令 $s=\frac mk$，表示操作的总次数，求出有多少种合法的操作方案即可，容斥一下好像可以完成，先胡一下。 人傻智力低，生成函数启动！！！ 操作一的生成函数是 $(\frac 1{1-x})^n$，操作二是 $(\frac{1-x^k}{1-x})^{n-k+1}$ 将这两个玩意暴力卷起来： $$ \begin{aligned} (\frac 1{1-x})^n(\frac{1-x^k}{1-x})^{n-k+1}&=(\frac{1}{1-x})^{2n-k+1}(1-x^k)^{n-k+1}\\ &=\Big( \sum_{i=0}^\infty{2n-k+i\choose 2n-k}x^i\Big)(\sum_{j=0}^{n-k+1} {n-k+1\choose j}(-1)^{j}x^{jk})\\ \end{aligned} $$ 要求第 $s$ 项，但是发现后面只有 $n-k+1$ 项所以直接枚举后边，组合数都可以 $O(n)$ 所以总时间是 $O(n^2)$。 ```cpp #include<bits/stdc++.h> #define int long long using namespace std; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint a) { return os << a.val; } }ans; int n,m,k; inline modint C(int n,int m){//O(m) if(m<0||n<m) return 0; modint res=1,inv=1; for(int i=1;i<=m;++i) res*=(n-i+1),inv*=i; return res*inv.inv(); } inline modint Get(int i){return C(2*n-k+i,2*n-k);} inline modint Get2(int j){return C(n-k+1,j)*((j)&1?-1:1);} signed main(){ cin>>n>>m>>k; if(m%k){cout<<0<<endl;return 0;} m/=k; for(int j=0;j<=n-k+1;++j){ ans+=Get2(j)*Get(m-j*k); } cout<<ans<<endl; return 0; } ``` ## [ARC162D](https://www.luogu.com.cn/problem/AT_arc162_d) 为了这个题才会的 Prufer，没救了…… 先考虑求出好树的个数，如果你会 Prufer 就不难发现为 ${n-2\choose d_1-1,d_2,d_3\dots d_n}=\frac{d_1(n-2)!}{\prod\limits_{i=1}^n{d_i!}}$。 根和叶子始终是好的，直接加上，现在对于每一个节点求其贡献。 假设我们要求节点 $u$ 的贡献，先想想给定其子树的点集 $T$ 时怎么做。 方案自然是 $\frac {d_u(|T|-2)!}{\prod\limits_{x\in T}d_x!}$，然后还要算子树外的节点，此时我们直接把 $d_u$ 子树当成一个点算 Prufer，可以得到是 $\frac {d_1(n-|T|-1)!}{\prod\limits_{x\not\in T}d_x!}$。 把两部分乘起来，变成了 $\frac{d_1d_u(|T|-2)!(n-|T|-1)!}{\prod\limits_{i=1}^nd_i!}$，只与 $|T|$ 和 $u$ 有关。 直接 DP 就行了，记得要从后往前 DP，因为子树内只能放编号大于当前点的。 注意 $|T|-1=\sum\limits_{x\in T} d_x$ 的条件，不然你无法构成一棵树。 ```cpp #include<bits/stdc++.h> using namespace std; const int N=510; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint a) { return os << a.val; } }dp[N][N],w,frac[N],tmp,res,T; int n,d[N]; int main(){ ios::sync_with_stdio(0);cin.tie(0);cout.tie(0); dp[0][0]=1;cin>>n;for(int i=n;i;--i,w=d[n]) cin>>d[i]; frac[0]=1;for(int i=1;i<=n;++i) frac[i]=frac[i-1]*i; tmp=w*frac[n-2];T=1;for(int i=1;i<=n;++i) T*=frac[d[i]].inv(); for(int i=1;i<=n;++i) if(!d[i]||i==n) res+=tmp*T; for(int i=1;i<n;++i){ for(int j=d[i];j<i;++j){ int x=j+1,y=j-d[i]; res+=frac[x-2]*frac[n-x-1]*T*d[i]*w*dp[j][y]; } for(int j=i-1;~j;--j){ for(int k=0;k<=n;++k){ if(!dp[j][k].val) continue; dp[j+1][k+d[i]]+=dp[j][k]; } } } cout<<res<<endl; return 0; } ``` ## [ARC163F](https://www.luogu.com.cn/problem/AT_arc163_f) 真是神仙题呀，看完题面的第一反应居然觉得是缝合怪，还是太菜了。 Increasing Problem 的标准做法是 slope trick，考虑前 $i$ 个数之后得到的分段函数是 $G_i(x)$，则令 $F_{i}(x)=G_{i-1}(x)+|a_i-x|$，然后让 $F_i(x)$ 做一下前缀取 $\min$ 就可以得到 $G_i(x)$ 了。 就是： ```cpp #include<bits/stdc++.h> #define endl '\n' using namespace std; typedef long long ll; const int N=1e6+10; int a[N],n,b[N];ll ans; priority_queue < int > S; int main(){ ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin>>n; for(int i=1;i<=n;++i){ cin>>a[i];S.push(a[i]); if(S.top()>a[i]){int x=S.top();S.pop();ans+=x-a[i];S.push(a[i]);} b[i]=S.top(); } cout<<ans<<endl; return 0; } ``` 思考 Many Increasing Problems 能不能类似做？ $n,m\le 10^5$ 有点吓人。 发现本质上算的是 $\sum a_i-\sum\limits_{x\in S}x$，$\sum a_i$ 是容易算的，考虑计算 $\sum\limits_{x\in S} x$。 形式化地考虑 $S$ 是怎么算的： 1. 朝 $S$ 中加入两个 $a_i$。 2. 弹出 $S$ 的最大值。 发现还是不好做，如果序列只有 $01$ 就好了。 考虑差分计算，将数列中 $\le x$ 的部分变成 $0$，$> x$ 的部分变成 $1$，那么每次有 $x$ 种方法令 $S$ 减少一个 $1$，有 $m-x$ 种方法使 $S$ 多一个 $1$。 考虑转到网格上，坐标 $(a,b)$ 表示考虑了前 $a$ 个点，目前有 $b$ 个 $1$，那么现在有两种走法： 1. 走到 $(a+1,\max(b-1,0))$，有 $x$ 种方案。 2. 走到 $(a+1,b+1)$，有 $m-x$ 种方案。 发现 $\max(b-1,0)$ 是一个很难刻画的条件，思考如何替代，这里有一个比较牛的 trick。 先不考虑对 $0$ 取 $\max$，假如一条路径走过的最低点的在 $0$ 下第 $z$ 个格子，那么就把路径整体抬高 $z$ 个格子，需要注意的是 $z$ 明显始终是非负数。 这样做为什么是对的呢？举个例子吧（红线是原路径，绿线是抬高后的路径）： ![](https://cdn.luogu.com.cn/upload/image_hosting/09vzk3rk.png) 所以终点其实是没变的，问题转化到了求从 $(0,a)$ 开始到 $(n,b)$ 结束，至少经过一次 $y=0$ 但不经过 $y=-1$ 的路径的方案数了。 发现确定起点和终点，可以直接做两次反射容斥求出方案。 发现我们在知道起点和终点的情况下减少 $1$ 和增加 $1$ 的次数都清楚，就只用看方案数就行了。方案数是不经过 $y=-1$ 的方案数减去不经过 $y=0$ 的方案数。 不经过 $y=-1$ 的方案数： $$ {n\choose \frac{n+b-a}{2}}-{n\choose\frac{n-b-a}{2}-1} $$ 不经过 $y=0$ 的方案数： $$ {n\choose \frac{n+b-a}{2}}-{n\choose\frac{n-b-a}{2}} $$ 做差得 ${n\choose\frac{n-b-a}{2}}-{n\choose\frac{n-b-a}{2}-1}$。 令 $i=\frac{n-b-a}2,j=\frac{n+b-a}2$，则有 $b=j-i$。 注意到 $i$ 的上界为 $\lfloor\frac{n}{2}\rfloor$，下界为 $0$，$j$ 的上界为 $n-i$，下界为 $i$。 写出式子，准备开心消消乐，注意 $0^0=1$： $$ \begin{aligned} Ans&=\sum_{x=0}^m\sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}\sum_{j=i}^{n-i}\left({n\choose i}-{n\choose i-1}\right)x^{n-j}(m-x)^j(j-i)\\ &=\sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}\sum_{j=i}^{n-i}(j-i)\left({n\choose i}-{n\choose i-1}\right)\sum_{x=0}^mx^{n-j}(m-x)^j\\ \end{aligned} $$ 令 $F(z)=\sum\limits_{i=0}^\infty z^i\sum\limits_{x=0}^mx^{n-i}(m-x)^i$，处理一下： $$ \begin{aligned} F(z)&=\sum_{i=0}^\infty z^i\sum_{x=0}^mx^{n-i}(m-x)^i\\ &=\sum_{i=0}^\infty \sum_{x=0}^mx^nx^{-i}(m-x)^iz^i\\ &=\sum_{x=0}^mx^n\sum_{i=0}^\infty(\frac{m-x}{x}z)^i\\ &=\sum_{x=0}^m\frac{x^n}{1-\frac{m-x}xz} \end{aligned} $$ 常规地 NTT 分治暴力合并分子分母就行了，这一部分是 $O(n\log^2n)$。 约定 $f_i=[z^i]F(z)$，继续往下看： $$ \begin{aligned} Ans &=\sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}\left({n\choose i}-{n\choose i-1}\right)\sum_{j=i}^{n-i}(j-i)f_j\\ \end{aligned} $$ 发现 $\sum\limits_{j=i}^{n-j}(j-i)f_j$ 是好预处理的，线性做了。 最后再线性遍历一次，总时间是 $O(n\log^2 n)$，瓶颈在于 NTT 分治。 ```cpp #include<bits/stdc++.h> using namespace std; namespace Fread { const int SIZE=1<<21;char buf[SIZE],*S,*T; inline char getchar() {if(S==T){T=(S=buf)+fread(buf,1,SIZE,stdin);if(S==T)return '\n';}return *S++;} } namespace Fwrite { const int SIZE=1<<21; char buf[SIZE],*S=buf,*T=buf+SIZE; inline void flush(){fwrite(buf,1,S-buf,stdout);S=buf;} inline void putchar(char c){*S++=c;if(S==T)flush();} struct POPOSSIBLE{~POPOSSIBLE(){flush();}}ztr; } #define getchar Fread :: getchar #define putchar Fwrite :: putchar namespace Fastio{ struct Reader{ template<typename T> Reader& operator >> (T& x) { char c=getchar();x=0; while(c<'0'||c>'9') c=getchar(); while(c>='0'&&c<='9'){x=x*10+(c-'0');c=getchar();} return *this; } Reader(){} }cin; struct Writer{ template<typename T> Writer& operator << (T x) { if(x==0){putchar('0');return *this;} static int sta[45];int top=0; while(x){sta[++top]=x%10;x/=10;} while(top){putchar(sta[top]+'0');--top;} return *this; } Writer& operator << (char c) {putchar(c);return *this;} Writer(){} }cout; } #define endl '\n' #define cin Fastio :: cin #define cout Fastio :: cout const int g=3; const int mod=998244353; const int invg=(mod+1)/3; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } inline int Mod(const int x){return x>=mod?x-mod:x;} modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint operator-() const { return modint(norm(-val)); } modint& operator+=(const modint& o) { return val = Mod(val + o.val), *this; } modint& operator-=(const modint& o) { return val = norm(val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend Fastio :: Writer& operator<<(Fastio :: Writer& os, modint a) { return os << a.val; } friend Fastio :: Reader& operator>>(Fastio :: Reader& is, modint& a) { return is >> a.val; } }; const int N=1e6+10; int rev[N],n,m; template < int Max > struct Choose{ modint frac[Max+10],inv[Max+10]; Choose(){ frac[0]=inv[0]=1;for(int i=1;i<=Max;++i) frac[i]=frac[i-1]*i; inv[Max]=frac[Max].inv();for(int i=Max-1;i;--i) inv[i]=inv[i+1]*(i+1); } modint operator ()(const int x,const int y){return x<y||y<0?0:frac[x]*inv[y]*inv[x-y];} }; Choose < N > C; inline void GetRev(const int lim){int bit=log2(lim)-1;for(int i=0;i<lim;++i) rev[i]=((rev[i>>1]>>1)|((i&1)<<bit));} inline int CalRev(const int len){int lim=1;for(;lim<=len;lim<<=1);GetRev(lim);return lim;} inline modint expow(modint x,int y){modint res=1;for(;y;y>>=1){if(y&1)res*=x;x*=x;}return res;} struct Poly{ modint a[N]; modint& operator [](const int x){return a[x];} inline void Readin(const int len){for(int i=0;i<len;++i) cin>>a[i].val;} inline void out(const int len){for(int i=0;i<len;++i) cout<<a[i]<<' ';cout<<endl;} inline void NTT(const int limit,const int type){ for(int i=0;i<limit;++i) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int i=1;i<limit;i<<=1){ modint Gn=expow(type==1?g:invg,(mod-1)/(i<<1)); for(int j=0;j<limit;j+=(i<<1)){ modint G=1; for(int k=0;k<i;++k,G*=Gn){ modint x=a[j+k],y=a[i+j+k]*G; a[j+k]=x+y;a[i+j+k]=x-y; } } } if(~type) return ; modint inv=expow(limit,mod-2); for(int i=0;i<limit;++i) a[i]*=inv; } }x,y; struct Frac{ Poly s[2]; inline void NTT(int lim,int op){s[0].NTT(lim,op);s[1].NTT(lim,op);} }s[20],res; inline void getinv(Poly &y,Poly &x,const int len){ static Poly s,tmp; int lim=2,bas=1;s[0]=expow(x[0],mod-2); while(bas<(len<<1)){ for(int i=0;i<bas;++i) tmp[i]=x[i]; GetRev(lim);tmp.NTT(lim,1);s.NTT(lim,1); for(int i=0;i<lim;++i) s[i]=s[i]*2-s[i]*s[i]*tmp[i]; s.NTT(lim,-1);for(int i=bas;i<lim;++i) s[i]=0; bas<<=1;lim<<=1; } for(int i=0;i<len;++i) y[i]=s[i]; for(int i=0;i<bas;++i) s[i]=tmp[i]=0; } void Solve(Frac &t,int dep,int l,int r){ if(l==r){ t.s[0][0]=expow(l,n); t.s[1][0]=1; t.s[1][1]=-((modint(l).inv())*(m-l)); return ; } int mid=(l+r)>>1,len=r-l+1; for(int i=0;i<=(len<<2);++i) t.s[0][i]=t.s[1][i]=0; Solve(s[dep],dep+1,l,mid); Solve(t,dep+1,mid+1,r); int lim=CalRev(len); s[dep].NTT(lim,1);t.NTT(lim,1); for(int i=0;i<lim;++i){ modint tmp=t.s[1][i]; t.s[0][i]=s[dep].s[0][i]*tmp+s[dep].s[1][i]*t.s[0][i]; t.s[1][i]=s[dep].s[1][i]*tmp; s[dep].s[0][i]=s[dep].s[1][i]=0; } t.NTT(lim,-1); } modint F[N],G[N]; inline modint Get(modint *v,int l,int r){return v[r]-(l?v[l-1]:0);} inline modint Calc(){ Solve(res,0,1,m);modint ans=0; for(int i=n+1;i<N;++i) res.s[0][i]=res.s[1][i]=0; getinv(res.s[1],res.s[1],n+1);int lim=CalRev((n+2)<<1); res.NTT(lim,1);for(int i=0;i<lim;++i) res.s[0][i]*=res.s[1][i]; res.s[0].NTT(lim,-1); for(int i=0;i<=n;++i){ F[i]=res.s[0][i]; G[i]=F[i]*i; if(i){ F[i]+=F[i-1]; G[i]+=G[i-1]; } } F[n]+=expow(m,n);G[n]+=expow(m,n)*n; // for(int i=0;i<=n/2;++i){ // for(int j=i;j<=n-i;++j){ // ans+=(C(n,i)-C(n,i-1))*(j-i)*F[j]; // } // } for(int i=0;i<=n/2;++i){ ans+=(C(n,i)-C(n,i-1))*(Get(G,i,n-i)-Get(F,i,n-i)*i); } // cout<<ans<<endl; return ans; } int main(){ modint ans=0;cin>>n>>m; for(int i=1;i<=m;++i) ans+=expow(m,n-1)*i*n; cout<<ans-Calc(); return 0; } ``` ## [ARC167C](https://www.luogu.com.cn/problem/AT_arc167_c) 提供一种瓶颈在排序的 $O(n\log n)$ 做法。 先排序明显对答案没有影响。 常规的，我们肯定要拆贡献来做，只要求出每条边在最小生成树中出现了几次，再加起来。 考虑这样一个问题： >两个点 $i,j$ 有边当且仅当 $|P_i-P_j|\le K$，且 $\max(P_i,P_j)\le x$，问最多能选几条边使得不出现环。 对于所有的排列，将该问题的答案求和，称为 $f_x$。 那么答案是 $\sum\limits_{i=1}^n (f_i-f_{i-1})A_i$。 考虑对于一个固定的排列，如何贪心求出 $f_x$： 对于所有 $p_j\le x$，将 $x$ 加入序列 $g$，然后将 $g$ 从小到大排序一下，答案即为 $\sum\limits_{i=2}^x[|g_i-g_{i-1}|\le K]$。 接下来考验基本功： 再次拆贡献，对于 $g_i$ 和 $g_{i-1}$ 求出它们贡献多少次，再乘上 $x-1$ 就是 $f_x$。理由很简单，发现对于所有 $i$，其贡献次数是一样的。 令 $k=g_i-g_{i-1}$，那么有 ${n-k\choose x-1}$ 个排列能使 $i$ 产生贡献。 所以 $f_x=x!(n-x)!(x-1)\sum\limits_{k=1}^K{n-k\choose x-1}$，$x!(n-x)!$ 是在分配 $p$。 $O(nK)$ 算出 $f$ 这个题就做完了，但是这个题可以优化到 $O(n)$？ 不对算上排序是 $O(n\log n)$。 具体来讲就是注意到 $\sum\limits_{k=1}^K{n-k\choose x-1}=\left({n\choose x}-{n-K\choose x}\right)$，然后优化了。 ```cpp #include<bits/stdc++.h> #define int long long using namespace std; const int N=5010; const int mod=998244353; const int inv2=(mod+1)/2; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint a) { return os << a.val; } }f[N],ans; template < int Max > struct Choose{ modint frac[Max+10],inv[Max+10]; Choose(){ frac[0]=inv[0]=1;for(int i=1;i<=Max;++i) frac[i]=frac[i-1]*i; inv[Max]=frac[Max].inv();for(int i=Max-1;i;--i) inv[i]=inv[i+1]*(i+1); } modint* operator [](const int x){return x?frac:inv;} modint operator ()(const int x,const int y){return x<y||y<0?0:frac[x]*inv[y]*inv[x-y];} }; int n,a[N],K; Choose < N > C; signed main(){ ios::sync_with_stdio(0);cin.tie(0);cout.tie(0); cin>>n>>K;for(int i=1;i<=n;++i) cin>>a[i]; sort(a+1,a+1+n); for(int x=1;x<=n;++x){ f[x]=C.frac[x]*C.frac[n-x]*(x-1)*(C(n,x)-C(n-K,x)); } for(int i=1;i<=n;++i) ans+=(f[i]-f[i-1])*a[i]; cout<<ans<<endl; return 0; } ``` ## [ARC176D](https://www.luogu.com.cn/problem/AT_arc176_d) 感觉很棒的 trick 啊！给出了这一类问题的通解。 单纯对于这个题有更简洁的做法。 对于任意一个位置 $x$ 来讲，其他值都是等价的。也就是说除了 $P_x$ 之外的值最后留在 $x$ 的概率是相等的。 注意到这一点之后，因为答案牵扯到两个位置，所以我们枚举两个位置 $x,y$ 讨论它们的权值关系，不妨假设 $x<y$，不难发现对于任意 $x<y$ 来讲，求出来的东西亦等价。 需要注意的是，我们并没有强调 $x,y$ 一定相邻，不过不难发现 $x,y$ 是否相邻最后结果都是一样的。 记 $a=P_x,b=P_y$，$c$ 表示不为 $a,b$ 的值，那么最后位置 $x,y$ 上的值有 $7$ 种情况：$(a,b),(b,a),(c,b),(b,c),(a,c),(c,a),(c,c)$。 考虑构造矩阵 $M$，其中 $M_{i,j}$ 表示原本是情况 $i$ 通过一次交换变成情况 $j$ 的方案数，构造如下： $$ \begin{bmatrix} {n-2\choose 2}&1&n-2&0&n-2&0&0\\ 1&{n-2\choose 2}&0&n-2&0&n-2&0\\ 1&0&{n-2\choose 2}+n-3&1&0&1&n-3\\ 0&1&1&{n-2\choose 2}+n-3&1&0&n-3\\ 1&0&0&1&{n-2\choose 2}+n-3&1&n-3\\ 0&1&1&0&1&{n-2\choose 2}+n-3&n-3\\ 0&0&1&1&1&1&{n-2\choose 2}+2(n-4)+1\\ \end{bmatrix} $$ 定义 $1\times 7$ 的行向量 $e$：$\begin{bmatrix} 1&0&0&0&0&0&0\end{bmatrix}$，也就是只有开头为一。 再定义行向量 $f=e\times M^{k}$，那么 $f_{1,i}$ 就是经过 $k$ 次操作后，$x,y$ 上值是 $i$ 种情况得方案数了，当然为了方便起见，在没有歧义的情况下我们用 $f_i$ 表示 $f_{1,i}$。 回归到本问题。 假设现在正在处理位置 $i,i=1$，考虑 $7$ 种情况分别对答案的贡献。 $(a,b),(b,a)$ 的贡献为：$(f_1+f_2)\times|P_i-P_{i+1}|$。 设 $s_x=\sum\limits_{i=1}^n |i-x|=\frac{x(x-1)}{2}+\frac{(n-x+1)(n-x)}{2}$，继续考虑接下来的贡献。 $(c,b),(b,c)$ 的贡献：$\frac{(f_3+f_4)\times(s_{P_{i+1}}-|P_{i}-P_{i+1}|)}{n-2}$。 除以 $n-2$ 是因为选出 $c$ 的方案已经在 $M$ 中算过了，但是在 $s$ 中又算了一遍所以要去重。 同理 $(a,c),(c,a)$ 的贡献：$\frac{(f_5+f_6)\times(s_{P_i}-|P_{i}-P_{i+1}|)}{n-2}$。 记 $sum=\sum\limits_{i=1}^n s_i$，继续计算 $(c,c)$ 的贡献为： $$ \frac{f_7\times(sum-2(s_{P_i}+s_{P_{i+1}}-|P_i-P_{i+1}|))}{(n-2)(n-3)} $$ 把所有贡献加起来即可。 ### 更简洁的做法 直接转差分，枚举阈值 $c$，将小于 $c$ 的当作 $0$，大于等于 $c$ 的当作 $1$。 那么一个 $01$ 串对答案产生的贡献就是 $01$ 交界处的个数。 然后发现情况数被压到了三种，直接沿用上面的做法即可！ 这里提供通解代码： ```cpp #include<bits/stdc++.h> using namespace std; const int N=2e5+10; const int mod=998244353; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint a) { return os << a.val; } }s[N],ans,sum,f[N]; int n,m,a[N]; const int Max=7; struct Matrix{ modint a[Max+1][Max+1]; modint* operator [](const int x){return a[x];} inline void clear(){memset(a,0,sizeof(a));} inline void Init(){for(int i=1;i<=Max;++i) for(int j=1;j<=Max;++j) a[i][j]=(i==j);} Matrix operator *(const Matrix &t){ Matrix z;z.clear(); for(int i=1;i<=Max;++i) for(int k=1;k<=Max;++k) for(int j=1;j<=Max;++j) z[i][j]+=a[i][k]*t.a[k][j]; return z; } }bas,res; inline void init(){ modint tmp=1ll*(n-2)*(n-3)/2;res[1][1]=1; bas=(Matrix){{ {0,0,0,0,0,0,0,0}, {0,tmp,1,n-2,0,n-2,0,0}, {0,1,tmp,0,n-2,0,n-2,0}, {0,1,0,tmp+n-3,1,0,1,n-3}, {0,0,1,1,tmp+n-3,1,0,n-3}, {0,1,0,0,1,tmp+n-3,1,n-3}, {0,0,1,1,0,1,tmp+n-3,n-3}, {0,0,0,1,1,1,1,tmp+2*(n-4)+1} }}; } int main(){ ios::sync_with_stdio(0);cin.tie(0);cout.tie(0); cin>>n>>m;init();for(int i=1;i<=n;++i) cin>>a[i]; for(int i=1;i<=n;++i){ s[i]=1ll*i*(i-1)/2+1ll*(n-i+1)*(n-i)/2; sum+=s[i]; } for(int i=m;i;i>>=1,bas=bas*bas) if(i&1) res=res*bas; for(int i=1;i<=7;++i) f[i]=res[1][i]; modint d=modint(n-2).inv(),div=modint(1ll*(n-2)*(n-3)).inv(); for(int i=1;i<n;++i){ modint v=abs(a[i]-a[i+1]),x=s[a[i]],y=s[a[i+1]]; ans+=(f[1]+f[2])*v; ans+=(f[3]+f[4])*(y-v)*d; ans+=(f[5]+f[6])*(x-v)*d; ans+=f[7]*(sum-(x+y-v)*2)*div; } cout<<ans<<endl; return 0; } ``` ## [ARC184E](https://www.luogu.com.cn/problem/AT_arc184_e) 怎么 B 题杨表，C 题奇怪注意题。还是做做 E 吧。 第一步就是想不到的转化：考虑将前缀和变成差分，称得到的新函数为 $g(i,j)$，那么由于前缀和与差分是逆运算，有关系 $f(i,j)=g(j,i)$。 现在考虑求 $\sum\limits_{i=1}^n\sum\limits_{j=1}^i g(i,j)$。 我们考虑进行一次差分相当于乘上一个 $1-x$，由于是模二意义下，所以也相当于乘 $1+x$。 考虑序列 $F$ 其被操作了 $k$ 次，那么它变成了 $(1+x)^kF$，将 $k$ 二进制拆分为 $2^{k_1}+2^{k_2}+\dots$，那么也等价于 $F\prod\limits_{i=1}(1+x)^{2^{k_i}}=F\prod\limits_{i=1}(1+x^{2^{k_i}})$，理由：${n\choose m}$ 是偶数当且仅当 $m\subset n$，使用 Lucas 即可证明。 发现模意义下的差分操作必然存在周期性，进一步地，每一个序列必然存在且仅存在一条入边和出边，证明： 1. 差分得到的结果唯一，所以有恰一条入边； 2. 前缀和得到结果唯一，所以恰一条出边。 所以如果把序列建成图，必然是由若干个简单环形成的。 如果我们对于每一个环确定一个代表元，并求出每个序列 $x$ 位置 $p_x$，那么对于两个序列 $a,b$，$g(a,b)$ 就是 $((p_a-p_b)\bmod l)$，$l$ 是环长。 问题落在了三件事上面： 1. 如何求 $p$。 2. 如何求 $l$。 3. 如何求在哪个环里。 先来确定什么是代表元，我们钦定一个环的代表元是环上字典序最小的一个序列。 那么对于一个序列 $a$，我们从小到大枚举 $(1+x^{2^{k_i}})$ 并找到最小的 $j$ 使得 $a_j=1$，那么每次变化的第一个位置一定是 $a_{j+2^{k_i}}$，我们可以根据这个来判断是否选上这个 $2^{k_i}$。 由上述过程可以得到 $p$，$l$ 明显是最小的 $2^k$ 使得 $i+2^k>m$，两个序列在一个环里当且仅当它们操作完后相同。 直接处理出 $p$ 之后使用树状数组即可。 代码被我吃了。 ## [ARC180F](https://www.luogu.com.cn/problem/AT_arc180_f) 传奇题啊！ 我们记 $a_n(y)$ 表示以下问题的答案： > 在 $[0,y]$ 上均匀随机 $n$ 个实数 $x_{1,\dots,n}$，若满足 $x_1<x_2<\dots<x_n$，则答案为 $\prod\limits_{i=1}^n(y^A+\sum\limits_{j=i+1}^n x_j^A)$，否则为 $0$。 注：此处如果保证 $x_1<x_2<\dots<x_n$ 的话貌似也可以做，只不过要把下面的 EGF 换成 OGF。 答案明显是 $n!\times a_n(1)$。 考虑记录 $a_n(y)$ 的生成函数 $f(y,x)=\sum\limits_{i\ge 0} x^ia_i(y)$。 考虑 $\prod\limits_{i=1}^n(y^A+\sum\limits_{j=i+1}^n x_j^A)$ 的组合意义是什么，相当于对于每一个 $i$ 选出一个 $j>i,x_j^A$ 或者 $y^A$，考虑刻画一个选择方案的值，如果 $i$ 选择了 $x_j^A$，则令 $i$ 向 $j$ 连边，否则向超级源点（姑且称作 $0$ 号节点），明显构成了一棵树，将节点 $j$ 的权值看作 $x_j^A$，那么值为每个子树的节点的权值之积再乘起来。 考虑拆贡献，来看子树 $i$ 的贡献是啥，设 $x_i=a$，那么其权值为 $\exp(a^Ax\int_0^af(t,x)\textrm dt)$。 又因为 $0$ 的点权值为 $y$，所以答案是 $f(y,x)=\exp(y^Ax\int_0^yf(t,x)\textrm d t)$。 所以不难发现 $f(y,x)$ 是关于 $x,y$ 的二元生成函数，更进一步地，$f(y,x)$ 可以表示成 $g(z)$，$z=y^{A+1}x$，以上两步 Atcoder 的题解似乎并没有给出形式证明，可以感性理解一下，更好的方式是喂给打 MO 的同学让他们讲。 接下来设 $g(z)=f(y,x)=\sum\limits_{i\ge 0}g_iz^i$，我们再设 $h(z)$ 等于括号里的那一坨东西，则 $h(z)=\sum\limits_{i\ge 0}\frac{g_i}{(A+1)i+1}z^{i+1}$，直接把后面的积分式算出来就行了。 所以又有 $g=\exp h$，直接大力递推即可。 ```cpp #include<bits/stdc++.h> using namespace std; const int N=1e5+10; const int mod=1e9+7; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } static constexpr int get_mod() { return mod; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint expow(int y){ modint res=1,x=val; for(;y;y>>=1,x*=x) if(y&1) res*=x; return res; } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint a) { return os << a.val; } }f[N],g[N]; int n,A; template < int Max > struct Choose{ modint frac[Max+10],inv[Max+10]; Choose(){ frac[0]=inv[0]=1;for(int i=1;i<=Max;++i) frac[i]=frac[i-1]*i; inv[Max]=frac[Max].inv();for(int i=Max-1;i;--i) inv[i]=inv[i+1]*(i+1); } modint* operator [](const int x){return x?frac:inv;} modint operator ()(const int x,const int y){return x<y||y<0?0:frac[x]*inv[y]*inv[x-y];} }; Choose < N > C; void Solve() { cin>>n>>A;g[0]=1; for(int i=1;i<=n;i++) { f[i]=g[i-1]*(modint(1ll*(A+1)*(i-1)+1).inv()); for(int j=1;j<=i;j++) g[i]+=f[j]*g[i-j]*j; g[i]*=C[0][i]*C[1][i-1]; } cout<<g[n]*C[1][n]<<'\n'; } int main(){ Solve(); return 0; } ``` ## [树的数量](https://www.luogu.com.cn/problem/P2767) 设 $[x^i]F(x)$ 表示有 $i$ 个点时答案为多少，特别地 $[x^0]F(x)=0$。 则有 $F(x)=x(1+F(x))^m$，故有 $x=\frac {F(x)}{(1+F(x))^m}$。 考虑拉格朗日反演，定义 $F(x)$ 的复合逆 $G(x)$ 满足 $G(F(x))=x=\frac{F(x)}{(1+F(x))^m}$。 有 $G(x)=\frac {x}{(1+x)^m}$，有 $[x^n]F(x)=\frac 1n[x^{n-1}]\left(\frac x{G(x)}\right)^n$。 有答案为 $\frac{nm\choose n-1}{n}$。 ```cpp #include<bits/stdc++.h> #define endl '\n' using namespace std; const int N=2e4+10; const int mod=10007; struct modint { int val; static int norm(const int& x) { return x < 0 ? x + mod : x; } modint inv() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) t = a / b, swap(a -= t * b, b), swap(u -= t * v, v); return modint(u); } modint() : val(0) {} modint(const int& m) : val(norm(m)) {} modint(const long long& m) : val(norm(m % mod)) {} modint operator-() const { return modint(norm(-val)); } bool operator==(const modint& o) { return val == o.val; } bool operator<(const modint& o) { return val < o.val; } modint& operator+=(const modint& o) { return val = (1ll * val + o.val) % mod, *this; } modint& operator-=(const modint& o) { return val = norm(1ll * val - o.val), *this; } modint& operator*=(const modint& o) { return val = static_cast<int>(1ll * val * o.val % mod), *this; } modint operator-(const modint& o) const { return modint(*this) -= o; } modint operator+(const modint& o) const { return modint(*this) += o; } modint operator*(const modint& o) const { return modint(*this) *= o; } friend std::ostream& operator<<(std::ostream& os, const modint a) { return os << a.val; } }; template < int Max > struct Choose{ modint frac[Max+10],inv[Max+10]; Choose(){ frac[0]=inv[0]=1;for(int i=1;i<=Max;++i) frac[i]=frac[i-1]*i; inv[Max]=frac[Max].inv();for(int i=Max-1;i;--i) inv[i]=inv[i+1]*(i+1); } modint operator ()(const int x,const int y){return x<y||y<0?0:frac[x]*inv[y]*inv[x-y];} }; int n,m; Choose < mod-1 > C; inline modint Cs(int x,int y){ if(x<mod) return C(x,y); else return Cs(x/mod,y/mod)*C(x%mod,y%mod); } int main(){ ios::sync_with_stdio(0);cin.tie(0);cout.tie(0); cin>>n>>m;cout<<(Cs(n*m,n-1)*modint(n).inv())<<endl; return 0; } ``` ## [生成树问题](https://www.luogu.com.cn/problem/P4002) $$ \begin{aligned} Ans&=\sum_{T}\left(\prod_{i=1}^n d_i^m\right)\left(\sum_{i=1}^nd_i^m\right)\left(\prod_{i=1}^n a_i^{c_i}\right)\\ \end{aligned} $$ 转而枚举 Prufer 序列而非枚举树。 用 $c_i$ 表示第 $i$ 个连通块在 Prufer 序列中出现的次数： $$ \begin{aligned} Ans&=\sum_{(\sum c)=n-2}\frac {(n-2)!}{\prod\limits_{i=1}^n(c_i!)}\sum_{i=1}^n(c_i+1)^{2m}a_i^{c_i}\left(\prod_{i\neq j} (c_j+1)^ma_j^{c_j+1}\right) \end{aligned} $$ 发现 $(\sum c)=n-2$ 这个限制非常特殊，能够快速处理这种限制的工具无疑是生成函数。 但是列出生成函数仍然是不容易的，考虑转换视角，考虑每一个 $i$ 带来的贡献并对这个贡献生成函数，进一步地，我们列 $2n$ 个生成函数形如： $$ \begin{aligned} A_i(x)&=\sum_{n\ge 0}\frac{x^n}{n!}(n+1)^{2m}a_i^{n+1}\\ B_i(x)&=\sum_{n\ge 0}\frac{x^n}{n!}(n+1)^{m}a_i^{n+1} \end{aligned} $$ 直接来处理 $A_i(x)$，发现 $(n+1)^{2m}$ 和 $a_i^{n+1}$ 都牵扯到 $n+1$。 套路地，此时可以走两条路：可以积分也可以选择直接展开。 考虑积分： $$ \begin{aligned} \int A_i(x)dx&=\sum_{n\ge 1}\frac{x^n}{n!}n^{2m}a_i^n\\ &=\sum_{n\ge 1}\frac{a_i^nx^n}{n!}\sum_{j=0}^n{n\choose j}{2m\brace j}j!\\ &=\sum_{n\ge 1}a_i^nx^n\sum_{j=0}^n{2m\brace j}\frac{1}{(n-j)!}\\ &=\sum_{j=0}^{2m}{2m\brace j}a_i^jx^j\sum_{n\ge j}\frac{(a_ix)^{n-j}}{(n-j)!}\\ &=\sum_{j=0}^{2m}{2m\brace j}a_i^jx^je^{a_ix} \end{aligned} $$ 再通过链式法则求导回去： $$ \begin{aligned} A_i(x)&=\sum_{j=1}^{2m}{2m\brace j}a_i^jjx^{j-1}e^{a_ix}+\sum_{j=0}^{2m}{2m\brace j}a_i^{j+1}x^je^{a_ix}\\ &=e^{a_ix}\sum_{j=0}^{2m}a_i^{j+1}x^j\left({2m\brace j}+(j+1){2m\brace j+1}\right)\\ &=e^{a_ix}\sum_{j=0}^{2m}a_i^{j+1}x^j{2m+1\brace j+1} \end{aligned} $$ 同理有： $$ B_i(x)=e^{a_ix}\sum_{j=0}^{m}a_i^{j+1}x^j{m+1\brace j+1} $$ 定义： $$ \begin{aligned} F_i(x)&=\sum_{j=0}^{2m}a_i^{j+1}x^j{2m+1\brace j+1}\\ G_i(x)&=\sum_{j=0}^{m}a_i^{j+1}x^j{m+1\brace j+1} \end{aligned} $$ $$ Ans=(n-2)![x^{n-2}](e^{\sum\limits_{i=1}^na_ix})\sum_{i=1}^nF_i(x)\prod_{i\neq j}G_j(x) $$ 求这个是容易的。