「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;
}