「HDU6593」Coefficient
题目
点这里看题目。
在后 open-hdu 的时代，我只能挂一个 vjudge 的链接充数了。
分析
首先观察到，和 \(x\) 具体关联的部分是 \(e^{ax+d}\)；而且我们还需要在 \(x_0=-\frac d a\) 处展开，这不是明示换元吗？
设 \(t=ax+d,g(t)=\frac{b}{e^t+c}\)，则有 \(f(x)=g(t)\)，且我们还有 \(\frac{f^{(n)}(x_0)}{n!}=a^n\frac{g^{(n)}(0)}{n!}\)。
事实上，由于我们这里不考虑收敛性，我们可以直接把 \(g(t)\) 当作形式幂级数处理，来求出它的 \([t^n]\)。这个结果便是所需的 \(\frac{g^{(n)}(0)}{n!}\)。
Note.
通常而言，封闭形式对应的完整的形式幂级数，对应的是它的麦克劳林级数。
观察一下特殊情况：如果 \(c=-1\)，则 \(g(0)\) 未定义，我们没有办法在 \(t_0=0\) 处展开它。因此，这种情况应当输出 \(-1\)。
在剩下的情况中，我们用一点小技巧，把整个问题转化成复合 \(e^x-1\) 的形式，也就是求 \([t^n]\frac{b}{(e^t-1)+(c+1)}\)。进一步的变形得到 \([t^n]\frac{\frac{b}{c+1}}{1-(-\frac{e^t-1}{c+1})}\)，然后展开：
\[\begin{aligned}
[t^n]\frac{1}{1-(-\frac{e^t-1}{c+1})}
&=[t^n]\sum_{k\ge 0}\left(-\frac{1}{c+1}\right)^k(e^t-1)^k\\
&=[t^n]\sum_{k=0}^{\bold{n}}\left(-\frac{1}{c+1}\right)^k\sum_{j=0}^k\binom{k}{j}e^{jt}(-1)^{k-j}\\
&=\sum_{k=0}^{n}\left(-\frac{1}{c+1}\right)^k\sum_{j=0}^k\binom{k}{j}\frac{j^n}{n!}(-1)^{k-j}\\
\end{aligned}
\]
这样的技巧的优秀之处在加粗的 \(n\) 处完全体现出来了。
Note.
事实是，构造良好的复合是一种很好的处理方法。
比如，你可能想问，为什么我们不选择直接展开成 \(\frac{\frac{b}{c}}{1-\left(-\frac{e^t}{c}\right)}\) 呢？
因为，这个时候复合的函数 \(-\frac{e^t}{c}\) 的级数形式不仅没有求和上界，还有常数项。这样的复合，实际上不够“良定义”的，是违背历史规律的行为。
那为什么 tly 直接复合 \(e^x\) 并解决本题呢？我想，这必然是因为祂是神，可以不受任何规律限制随心所欲地使用任何方式解决任何问题。
其它的例子包括斯特林数处理幂和问题。斯特林数的生成函数是：
\[\exp(y(e^x-1)) \]当我们具体处理幂和的时候，我们只需要强行将 \(e^x\) 拆成 \((e^x-1)+1\)，就可以得到斯特林数，并且还顺便缩减了求和的界。
接下来，我们可以注意到，第二个 \(\sum\) 里面的内容和 \(c\) 完全没有关系，则这实际上是关于 \(-\frac 1{c+1}\) 的多项式函数。算多项式的系数可以卷积处理，而算多次点值直接多点求值就好啦！
时间复杂度是 \(O(q\log^2 q)\)。
代码
#include <cstdio>
#include <vector>
#include <algorithm>
#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 mod = 998244353;
const int MAXN = ( 1 << 17 ) + 5;
template<typename _T>
void read( _T &x ) {
x = 0; char s = getchar(); int f = 1;
while( ! ( '0' <= s && s <= '9' ) ) { f = 1; if( s == '-' ) f = -1; s = getchar(); }
while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
x *= f;
}
template<typename _T>
void write( _T x ) {
if( x < 0 ) putchar( '-' ), x = -x;
if( 9 < x ) write( x / 10 );
putchar( x % 10 + '0' );
}
template<typename _T>
inline _T Max( const _T &a, const _T &b ) {
return a > b ? a : b;
}
typedef std :: vector<int> Poly;
int fac[MAXN], ifac[MAXN];
Poly prod[MAXN << 2];
int ans[MAXN], xVal[MAXN], id[MAXN];
Poly U, V, F;
int N, Q, M;
inline int Qkpow( int, int );
inline int Inv( const int &a ) { return Qkpow( a, mod - 2 ); }
inline int Mul( int x, const int &v ) { return 1ll * x * v % mod; }
inline int Sub( int x, const int &v ) { return ( x -= v ) < 0 ? x + mod : x; }
inline int Add( int x, const int &v ) { return ( x += v ) >= mod ? x - mod : x; }
inline int& MulEq( int &x, const int &v ) { return x = 1ll * x * v % mod; }
inline int& SubEq( int &x, const int &v ) { return ( x -= v ) < 0 ? ( x += mod ) : x; }
inline int& AddEq( int &x, const int &v ) { return ( x += v ) >= mod ? ( x -= mod ) : x; }
inline int Qkpow( int base, int indx ) {
int ret = 1;
while( indx ) {
if( indx & 1 ) MulEq( ret, base );
MulEq( base, base ), indx >>= 1;
}
return ret;
}
namespace Basics {
const int L = 17, g = 3, phi = mod - 1;
int w[MAXN];
inline void NTTInit( const int &n = 1 << L ) {
w[0] = 1, w[1] = Qkpow( g, phi >> L );
rep( i, 2, n - 1 ) w[i] = Mul( w[i - 1], w[1] );
}
inline void DIF( int *coe, const int &n ) {
int *wp, p, e, o;
for( int s = n >> 1 ; s ; s >>= 1 )
for( int i = 0 ; i < n ; i += s << 1 ) {
p = ( 1 << L ) / ( s << 1 ), wp = w;
for( int j = 0 ; j < s ; j ++, wp += p ) {
e = coe[i + j], o = coe[i + j + s];
coe[i + j] = Add( e, o );
coe[i + j + s] = Mul( *wp, Sub( e, o ) );
}
}
}
inline void DIT( int *coe, const int &n ) {
int *wp, p, k;
for( int s = 1 ; s < n ; s <<= 1 )
for( int i = 0 ; i < n ; i += s << 1 ) {
p = ( 1 << L ) / ( s << 1 ), wp = w;
for( int j = 0 ; j < s ; j ++, wp += p )
k = Mul( *wp, coe[i + j + s] ),
coe[i + j + s] = Sub( coe[i + j], k ),
coe[i + j] = Add( coe[i + j], k );
}
std :: reverse( coe + 1, coe + n );
int inv = Inv( n ); rep( i, 0, n - 1 ) MulEq( coe[i], inv );
}
}
namespace Operation {
inline Poly PlusMul( const Poly &A, const Poly &B, const int &rem ) {
static int P[MAXN] = {}, Q[MAXN] = {};
int n = A.size(), m = B.size(), L;
for( L = 1 ; L <= n + m - 2 ; L <<= 1 );
rep( i, 0, L - 1 ) P[i] = Q[i] = 0;
rep( i, 0, n - 1 ) P[i] = A[i];
rep( i, 0, m - 1 ) Q[i] = B[i];
Basics :: DIF( P, L );
Basics :: DIF( Q, L );
rep( i, 0, L - 1 ) MulEq( P[i], Q[i] );
Basics :: DIT( P, L ); Poly ret;
rep( i, 0, rem - 1 ) ret.push_back( P[i] );
return ret;
}
inline Poly SubMul( const Poly &A, const Poly &B, const int &rem ) {
static int P[MAXN] = {}, Q[MAXN] = {};
int n = A.size(), m = B.size(), L;
for( L = 1 ; L <= n + m - 2 ; L <<= 1 );
rep( i, 0, L - 1 ) P[i] = Q[i] = 0;
rep( i, 0, n - 1 ) P[i] = A[i];
rep( i, 0, m - 1 ) Q[i] = B[i];
std :: reverse( Q, Q + m );
Basics :: DIF( P, L );
Basics :: DIF( Q, L );
rep( i, 0, L - 1 ) MulEq( P[i], Q[i] );
Basics :: DIT( P, L ); Poly ret;
rep( i, 0, rem - 1 ) ret.push_back( P[i + m - 1] );
return ret;
}
}
namespace PolyInv {
int P[MAXN], Q[MAXN], U[MAXN], V[MAXN];
void Newton( const int &n ) {
if( n == 1 ) {
U[0] = Inv( P[0] );
return ;
}
int h = ( n + 1 ) >> 1, L; Newton( h );
for( L = 1 ; L <= n + h - 2 ; L <<= 1 );
rep( i, 0, L - 1 ) Q[i] = V[i] = 0;
rep( i, 0, n - 1 ) Q[i] = P[i];
rep( i, 0, h - 1 ) V[i] = U[i];
Basics :: DIF( Q, L );
Basics :: DIF( V, L );
rep( i, 0, L - 1 ) Q[i] = Mul( V[i], Sub( 2, Mul( V[i], Q[i] ) ) );
Basics :: DIT( Q, L );
rep( i, h, n - 1 ) U[i] = Q[i];
}
inline Poly PolyInv( const Poly &A, const int &rem ) {
int n = A.size();
rep( i, 0, rem - 1 ) P[i] = 0;
rep( i, 0, n - 1 ) P[i] = A[i];
Newton( rem ); Poly ret;
rep( i, 0, rem - 1 ) ret.push_back( U[i] );
return ret;
}
}
inline void Init( const int &n = 1 << 17 ) {
fac[0] = 1; rep( i, 1, n ) fac[i] = Mul( fac[i - 1], i );
ifac[n] = Inv( fac[n] ); per( i, n - 1, 0 ) ifac[i] = Mul( ifac[i + 1], i + 1 );
}
void Divide( const int &x, const int &l, const int &r ) {
using namespace Operation;
if( l > r ) return ;
if( l == r ) {
prod[x].clear(), prod[x].resize( 2 );
prod[x][0] = 1, prod[x][1] = Mul( mod - 1, xVal[l] );
return ;
}
int mid = ( l + r ) >> 1;
Divide( x << 1, l, mid );
Divide( x << 1 | 1, mid + 1, r );
prod[x] = PlusMul( prod[x << 1], prod[x << 1 | 1], r - l + 2 );
}
void Redivide( const int &x, const int &l, const int &r, const Poly &cur ) {
using namespace Operation;
if( l > r ) return ;
if( l == r ) { MulEq( ans[id[l]], cur[0] ); return ; }
int mid = ( l + r ) >> 1;
Redivide( x << 1, l, mid, SubMul( cur, prod[x << 1 | 1], mid - l + 2 ) );
Redivide( x << 1 | 1, mid + 1, r, SubMul( cur, prod[x << 1], r - mid + 1 ) );
}
int main() {
freopen( "coe.in", "r", stdin );
freopen( "coe.out", "w", stdout );
Init();
Basics :: NTTInit();
while( ~ scanf( "%d %d", &N, &Q ) ) {
U.clear(), V.clear();
U.resize( N + 1 ), V.resize( N + 1 );
rep( i, 0, N ) {
U[i] = Mul( Qkpow( i, N ), ifac[i] );
V[i] = i & 1 ? Mul( mod - 1, ifac[i] ) : ifac[i];
}
F = Operation :: PlusMul( U, V, N + 1 );
rep( i, 0, N ) MulEq( F[i], fac[i] );
M = 0;
rep( i, 1, Q ) {
int a, b, c, d;
read( a ), read( b );
read( c ), read( d );
a = ( a % mod + mod ) % mod;
b = ( b % mod + mod ) % mod;
c = ( c % mod + mod ) % mod;
d = ( d % mod + mod ) % mod;
if( c == mod - 1 ) ans[i] = -1;
else {
ans[i] = Mul( Mul( Qkpow( a, N ), Mul( b, Inv( Add( c, 1 ) ) ) ), ifac[N] );
M ++, xVal[M] = Mul( mod - 1, Inv( Add( c, 1 ) ) ), id[M] = i;
}
}
Divide( 1, 1, M );
U = PolyInv :: PolyInv( prod[1], Max( M + 1, N + 1 ) );
F = Operation :: SubMul( F, U, M + 1 );
Redivide( 1, 1, M, F );
rep( i, 1, Q ) write( ans[i] ), putchar( '\n' );
}
return 0;
}
