• 题解 仙人掌计数


    题目传送门

    题目大意

    给出(q)个查询,每次查询(n)个点的无根有标号仙人掌有多少个。

    (qle 5 imes 10^4,n< 131072)

    思路

    因为这道题太难码了,所以先把题解写了再写代码(好奇怪啊)终于码出来了,果然还是( ext {vector})好用(雾

    为了方便,我们下面的答案其实求的是有根有标号的答案,最后除以(n)就好了。我们设(c_n)(n)个点时的答案,我们对其构造指数型生成函数(C(x)):

    [C(x)=sum_{i=0} c_ifrac{x^i}{i!} ]

    我们考虑钦定一个根,然后连边,显然要么连出去一条边,要么连出去一个环。需要注意的是,连边相当于把根与一个仙人掌相连,连环的环其实是仙人掌连成一个环。

    于是,我们可以得到转移式:

    [C(x)=xexp (C(x)+frac{1}{2}sum_{i=2}^{infty} C^i(x)) ]

    除以(2)是因为一个环无论正反都是一样,但是我们计算时却视作两种情况。

    将上式用等比数列求和公式可以得到:

    [C(x)=xexp(frac{2C(x)-C^2(x)}{2-2C(x)}) ]

    我们发现我们如果构造

    [F(C(x))=xexp(frac{2C(x)-C^2(x)}{2-2C(x)})-C(x) ]

    我们即是要求导(F(C(x)))的根,对此,我们可以使用多项式牛顿迭代法。我们就需要求到(F(C(x)))的导:

    [F^{'}(C(x))=(xexp(frac{2C(x)-C^2(x)}{2-2C(x)})-C(x))^{'} ]

    [=xexp (frac{2C(x)-C^2(x)}{2-2C(x)})(frac{2C(x)-C^2(x)}{2-2C(x)})^{'}-1 ]

    [=xexp (frac{2C(x)-C^2(x)}{2-2C(x)})frac{C^2(x)-2C(x)+2}{2C^2(x)-4C(x)+2}-1 ]

    于是,我们可以得到牛顿迭代法的式子:

    [C(x)=C_0(x)-frac{F(C_0(x))}{F^{'}(C_0(x))} ]

    [=C_0(x)-frac{xexp(frac{2C_0(x)-C_0^2(x)}{2-2C_0(x)})-C(x)}{xexp (frac{2C_0(x)-C_0^2(x)}{2-2C_0(x)})frac{C_0^2(x)-2C_0(x)+2}{2C_0^2(x)-4C_0(x)+2}-1} ]

    [=C_0(x)-frac{2xexp (frac{2C_0(x)-C_0^2(x)}{2-2C_0(x)})-2C_0(x)}{(1+frac{1}{(C_0(x)-1)^2})xexp (frac{2C_0(x)-C_0^2(x)}{2-2C_0(x)})-2} ]

    我们发现如果我们设(G=xexp (frac{2C_0(x)-C_0^2(x)}{2-2C_0(x)})),那么式子就可以化为:

    [C(x)=C_0(x)-frac{2G-2C_0(x)}{(1+frac{1}{(C_0(x)-1)^2})G-2} ]

    ( ext {Code})

    #pragma GCC optimize("Ofast")
    #pragma GCC optimize("inline", "no-stack-protector", "unroll-loops")
    #pragma GCC diagnostic error "-fwhole-program"
    #pragma GCC diagnostic error "-fcse-skip-blocks"
    #pragma GCC diagnostic error "-funsafe-loop-optimizations"
    
    #include <bits/stdc++.h>
    using namespace std;
    
    #define SZ(x) ((int)x.size())
    #define Int register int
    #define mod 998244353
    #define MAXN 1000005
    
    int mul (int a,int b){return 1ll * a * b % mod;}
    int dec (int a,int b){return a >= b ? a - b : a + mod - b;}
    int add (int a,int b){return a + b >= mod ? a + b - mod : a + b;}
    int qkpow (int a,int k){
    	int res = 1;for (;k;k >>= 1,a = 1ll * a * a % mod) if (k & 1) res = 1ll * res * a % mod;
    	return res;
    }
    int inv (int x){return qkpow (x,mod - 2);}
    
    typedef vector <int> poly;
    
    int w[MAXN],rev[MAXN];
    
    void init_ntt (){
    	int lim = 1 << 18;
    	for (Int i = 0;i < lim;++ i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << 17);
    	int Wn = qkpow (3,(mod - 1) / lim);w[lim >> 1] = 1;
    	for (Int i = lim / 2 + 1;i < lim;++ i) w[i] = mul (w[i - 1],Wn);
    	for (Int i = lim / 2 - 1;i;-- i) w[i] = w[i << 1];
    }
    
    void ntt (poly &a,int lim,int type){
    #define G 3
    #define Gi 332748118
    	static unsigned long long d[MAXN];
    	for (Int i = 0,z = 18 - __builtin_ctz(lim);i < lim;++ i) d[rev[i] >> z] = a[i];
    	for (Int i = 1;i < lim;i <<= 1)
    		for (Int j = 0;j < lim;j += i << 1)
    			for (Int k = 0;k < i;++ k){
    				int x = d[i + j + k] * w[i + k] % mod;
    				d[i + j + k] = d[j + k] + mod - x,d[j + k] += x;
    			}
    	for (Int i = 0;i < lim;++ i) a[i] = d[i] % mod;
    	if (type == -1){
    		reverse (a.begin() + 1,a.begin() + lim);
    		for (Int i = 0,Inv = inv (lim);i < lim;++ i) a[i] = mul (a[i],Inv);
    	}
    #undef G
    #undef Gi 
    }
    
    poly operator + (poly a,poly b){
    	a.resize (max (SZ (a),SZ (b)));
    	for (Int i = 0;i < SZ (b);++ i) a[i] = add (a[i],b[i]);
    	return a;
    }
    
    poly operator - (poly a,poly b){
    	a.resize (max (SZ (a),SZ (b)));
    	for (Int i = 0;i < SZ (b);++ i) a[i] = dec (a[i],b[i]);
    	return a;
    }
    
    poly operator * (poly a,int b){
    	for (Int i = 0;i < SZ (a);++ i) a[i] = mul (a[i],b);
    	return a;
    }
    
    poly operator * (poly a,poly b){
    	int d = SZ (a) + SZ (b) - 1,lim = 1;while (lim < d) lim <<= 1;
    	a.resize (lim),b.resize (lim);
    	ntt (a,lim,1),ntt (b,lim,1);
    	for (Int i = 0;i < lim;++ i) a[i] = mul (a[i],b[i]);
    	ntt (a,lim,-1),a.resize (d);
    	return a;
    }
    
    poly operator << (poly a,int n){
    	a.resize (SZ (a) + n);
    	for (Int i = SZ (a) - 1;~i;-- i) a[i] = (i >= n ? a[i - n] : 0);
    	return a;
    }
    
    poly inv (poly a,int n){
    	poly b(1,inv (a[0])),c;
    	for (Int l = 4;(l >> 2) < n;l <<= 1){
    		c.resize (l >> 1);
    		for (Int i = 0;i < (l >> 1);++ i) c[i] = i < n ? a[i] : 0;
    		c.resize (l),b.resize (l);
    		ntt (c,l,1),ntt (b,l,1);
    		for (Int i = 0;i < l;++ i) b[i] = mul (b[i],dec (2,mul (b[i],c[i])));
    		ntt (b,l,-1),b.resize (l >> 1);
    	}
    	b.resize (n);
    	return b;
    }
    
    poly inv (poly a){return inv (a,SZ (a));}
    poly der (poly a){
    	for (Int i = 0;i < SZ (a) - 1;++ i) a[i] = mul (a[i + 1],i + 1);
    	a.pop_back ();return a;
    }
    poly ine (poly a){
    	a.push_back (0);
    	for (Int i = SZ (a) - 1;i;-- i) a[i] = mul (a[i - 1],inv (i));
    	a[0] = 0;return a;
    }
    
    poly ln (poly a,int n){
    	a = ine (der (a) * inv (a));
    	a.resize (n);
    	return a;
    }
    poly ln (poly a){return ln (a,SZ (a));}
    
    poly exp (poly a,int n){
    	poly b (1,1),c;
    	for (Int l = 2;(l >> 1) < n;l <<= 1){
    		b.resize (l),c = ln (b);
    		for (Int i = 0;i < l;++ i) c[i] = dec (i < n ? a[i] : 0,c[i]);
    		c[0] = add (c[0],1);
    		b = b * c,b.resize (l);
    	}
    	b.resize (n);
    	return b;
    }
    poly exp (poly a){return exp (a,SZ (a));}
    
    int n,m,x;poly A;
    
    poly work (int n){
    	poly a({0,1}),b,c;
    	for (Int l = 4;(l >> 1) < n;l <<= 1){
    		a.resize (l),c = inv (poly {1} - a),b = c * a,b.resize (l),b = (b + a) * inv (2);
    		b = exp (b),b = b << 1,c = c * c,c[0] = add (c[0],1),c.resize (l),c = c * inv (2) * b,c[0] = dec (c[0],1);
    		b = b - a;while (!b[0] && !c[0]) b.erase (b.begin()),c.erase (c.begin());
    		b.resize (l),c.resize (l);
    		a = a - b * inv (c),a.resize (l);
    	}
    	a.resize (n);
    	return a;
    }
    
    template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c == '-') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
    template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
    template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
    
    int q,len,fac[MAXN],que[MAXN];
    
    signed main(){
    	init_ntt(),read (q);for (Int i = 1;i <= q;++ i) read (que[i]),len = max (len,que[i]);
    	A = work (len);fac[0] = 1;for (Int i = 1;i <= len;++ i) fac[i] = mul (fac[i - 1],i);
    	for (Int i = 1;i <= q;++ i) write (1ll * A[que[i]] * fac[que[i] - 1] % mod),putchar ('
    ');
    	return 0;
    }
    

    话说预处理原根真的好快啊!!!

  • 相关阅读:
    雅虎公司C#笔试题及参考答案
    设计模式学习草稿:简单工厂,工厂,抽象工厂
    个人理解这三种工厂模式
    使用C#的BitmapData
    C#设计模式视频教程(不知道讲的好不好,刚刚看到)
    【HANA系列】SAP HANA XS的JavaScript API详解
    【HANA系列】SAP HANA XS使用服务器JavaScript Libraries详解
    【HANA系列】SAP HANA XS的JavaScript安全事项
    【HANA系列】SAP HANA XS使用JavaScript编程详解
    【HANA系列】SAP HANA XS使用JavaScript数据交互详解
  • 原文地址:https://www.cnblogs.com/Dark-Romance/p/13292609.html
Copyright © 2020-2023  润新知