「LOJ6737」指数公式
题目
点这里看题目。
给定一个集合 \(S\subseteq [n]\)，其中 \([n]=\{1,2,3,\dots,n-1,n\}\)。设 \(F(x)=\sum_{k\in S}\frac{x^k}{k!}\)。
对于 \(1\le k\le n\)，计算 \((k![x^k]\exp F(x)) \bmod 2\)。
所有数据满足 \(1\le n\le 2\times 10^6\)。
分析
这题唯一难点就是如何在 \(\mathbb F_2\) 下算系数。
显然，我们不可能把阶乘逆元这种东西带着计算，所以考虑将一个指数型生成函数表示成 \(\sum_{k\ge 0}f_k\frac{x^k}{k!}\) 的形式。因为最后提取的就是 \(f\)，所以这个 \(f\) 可以直接对 \(2\) 取模。
在这种意义下，生成函数常见运算的含义是：
-
\(F(x)+G(x)\) 对应于系数的对应项求和。
-
\(F(x)G(x)\) 对应于系数的二项卷积：
\[F(x)G(x)=\sum_{n\ge 0}\frac{x^n}{n!}\sum_{k=0}^n\binom{n}{k}f_kg_{n-k} \] -
形式积分和形式微分对应于系数的平移。
有了基础运算后，我们姑且不顾复杂度，直接对于 \(\exp F(x)\) 施加牛顿迭代法。
考察 \(\exp F(x)\) 迭代需要的东西，除了上面的基础运算外还附加一个多项式 \(\ln\)。不过 \(\ln\) 计算时只需要积分、微分和求逆，求逆也可以用牛顿迭代算。
更进一步地，对于形式幂级数 \(F(x)\in \mathbb F_2[[x]]\)，如果 \(F(x)\) 本身是可逆的，则 \(F(x)^2=0\)，也即 \(F(x)\) 的逆即为其自身。这一点不难证明。
现在考虑二项卷积如何计算。注意到 \(\dbinom{n}{k}\bmod 2=1\) 当且仅当二进制集合意义下 \(k\subseteq n\)，所以形式幂级数的二项卷积与其对应集合幂级数的子集卷积相对应。于是，长度为 \(2^m\) 的两个形式幂级数的二项卷积可以在 \(O(2^mm^2)\) 的时间内计算。
此时，整个的复杂度为 \(O(n\log^2n)\)，还不够快。我们注意到子集卷积有一种基于占位多项式的计算方法，此处占位多项式系数为 \(01\)，并且度数不高，于是可以想到将系数压成整数，并用其模拟占位多项式的运算，从而做到 \(O(\frac{n\log^2n}{\omega})\)。
代码
#include <cstdio>
#include <cstring>
#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )
const int MAXN = ( 1 << 21 ) + 5;
template<typename _T>
inline void Read( _T &x ) {
x = 0; char s = getchar(); bool f = false;
while( ! ( '0' <= s && s <= '9' ) ) { f = s == '-', s = getchar(); }
while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
if( f ) x = -x;
}
template<typename _T>
inline void Write( _T x ) {
if( x < 0 ) putchar( '-' ), x = -x;
if( 9 < x ) Write( x / 10 );
putchar( x % 10 + '0' );
}
char buf[MAXN];
unsigned F[MAXN];
int N;
inline unsigned Mul( const unsigned &u, const unsigned &v ) {
if( ! v ) return 0;
int ret = 0;
for( unsigned x = v ; x ; x &= x - 1 )
ret ^= u << __builtin_ctz( x );
return ret;
}
inline void Transform( unsigned *coe, const int &n ) {
for( int s = 1 ; s < n ; s <<= 1 )
for( int i = 0 ; i < n ; i += s << 1 )
for( int j = 0 ; j < s ; j ++ )
coe[i | j | s] ^= coe[i | j];
}
namespace PolyInv {
unsigned P[MAXN], U[MAXN], Q[MAXN], V[MAXN];
inline void PolyInv( unsigned *ret, const unsigned *A, const int &n ) {
rep( i, 0, n - 1 ) ret[i] = A[i];
}
}
namespace PolyLn {
unsigned P[MAXN], PInv[MAXN];
inline void PolyLn( unsigned *ret, const unsigned *A, const int &n ) {
int L = 1;
for( ; L < n ; L <<= 1 );
rep( i, 0, L - 1 ) P[i] = PInv[i] = 0;
rep( i, 1, n - 1 ) P[i - 1] = A[i] << __builtin_popcount( i - 1 );
PolyInv :: PolyInv( PInv, A, n );
rep( i, 0, n - 1 ) PInv[i] <<= __builtin_popcount( i );
Transform( P, L );
Transform( PInv, L );
rep( i, 0, L - 1 ) P[i] = Mul( P[i], PInv[i] );
Transform( P, L );
rep( i, 0, n - 2 ) ret[i + 1] = P[i] >> __builtin_popcount( i ) & 1;
ret[0] = 0;
}
}
namespace PolyExp {
unsigned P[MAXN], U[MAXN], Q[MAXN], V[MAXN];
inline void PolyExp( unsigned *ret, const unsigned *A, const int &n ) {
rep( i, 0, n - 1 ) P[i] = A[i];
U[0] = 1;
for( int L = 1 ; L < n ; L <<= 1 ) {
PolyLn :: PolyLn( Q, U, L << 1 );
rep( i, 0, L - 1 )
V[i] = U[i] << __builtin_popcount( i ),
Q[i] = ( Q[i | L] ^ P[i | L] ) << __builtin_popcount( i );
Transform( V, L );
Transform( Q, L );
rep( i, 0, L - 1 ) V[i] = Mul( V[i], Q[i] );
Transform( V, L );
rep( i, 0, L - 1 )
U[i | L] = V[i] >> __builtin_popcount( i ) & 1;
}
rep( i, 0, n - 1 ) ret[i] = U[i];
}
}
int main() {
scanf( "%s", buf + 1 );
N = strlen( buf + 1 );
int L = 1;
for( ; L <= N ; L <<= 1 );
rep( i, 1, N ) F[i] = buf[i] == '1';
PolyExp :: PolyExp( F, F, L );
rep( i, 1, N ) putchar( F[i] + '0' );
puts( "" );
return 0;
}
