题目
点这里看题目。
分析
由于本题明显涉及到了两个需要控制的变量,也就是序列长度和序列权值,我们引入二元生成函数,用 \(x\) 的指数描述长度,用 \(y\) 的指数描述权值。
容易发现,序列上可以被划分成两种极长连续段的交错组合:黑白段和灰色段。考虑到两个黑白段之间必然存在非空灰色段,我们把灰色段和黑白段打包并写出生成函数:
\[F(x,y)=\frac{x}{1-x}\cdot\frac{2x}{1-xy}
\]
这样,即可给出整个序列的生成函数:
\[\frac{1}{1-F}\cdot \frac{1}{1-x}=\frac{1-xy}{1-(y+1)x+(y-2)x^2}
\]
我们要求的实际上就是:
\[[x^{n+1}]\frac{1-xy}{1-(y+1)x+(y-2)x^2}\pmod{y^{k+1}}
\]
注意到一点:如果把 \(x\) 看作是主元,则我们实际上算的就是一个常系数齐次线性递推数列的远处的值,因此理论上来说可以直接套用那个算法,复杂度不清楚,具体还需移步 luogu 题解区。
但是,作为学过高中数学的同学,我们应当可以熟练地对于分母进行因式分解并直接得到通项,就像处理 Fibonacci 数列的生成函数那样。
不难解得:
\[\begin{aligned}
\Delta&=y^2-2y+9\\
\phi&=\frac{y+1+\sqrt \Delta}{2}\\
\hat\phi&=\frac{y+1-\sqrt \Delta}{2}\\
1-(y+1)x+(y-2)x^2&=(1-\phi x)(1-\hat\phi x)
\end{aligned}
\]
然后只需要再待定系数解一下即可得到:
\[\begin{aligned}
\frac{1-xy}{1-(y+1)x+(y-2)x^2}=\frac{1}{2\sqrt \Delta}\left(\frac{\sqrt\Delta+1-y}{1-\phi x}+\frac{\sqrt\Delta+y-1}{1-\hat\phi x}\right)
\end{aligned}
\]
因此可以很轻易地得到最终的答案为:
\[\frac{1}{2\sqrt\Delta}\left((\sqrt\Delta +1-y)\left(\frac{y+1+\sqrt\Delta}{2}\right)^{n+1}+(\sqrt\Delta + y-1)\left(\frac{y+1-\sqrt\Delta}{2}\right)^{n+1}\right)\pmod{y^{k+1}}
\]
直接硬算即可,需要用到多项式开根和多项式快速幂,然而最终复杂度仍然是 \(O(k\log k)\) 的。
小结:
-
适时引入新的元来表达必须表达的信息;
-
推生成函数,这种不太困难的主要就是看对于组合情景理解是否到位。基本功要扎实。
-
这样的选主元的想法比较基础;关键在于选取哪一个作主元。
不过既然这里有明确的题目背景:要求关于 \(y\) 的多项式,且这个多项式还是 \(x\) 的远处系数,也就不难想到让 \(x\) 成为主元。
想到这一点,之后的数学运算其实也就比较基础了。
代码
#include <cstdio>
#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 -- )
typedef long long LL;
const int mod = 998244353, inv2 = ( mod + 1 ) >> 1;
const int MAXN = 1 << 18;
template<typename _T>
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>
void write( _T x ) {
if( x < 0 ) putchar( '-' ), x = -x;
if( 9 < x ) write( x / 10 );
putchar( x % 10 + '0' );
}
int B[MAXN], C[MAXN], ans[MAXN];
int D[MAXN];
LL N; int K;
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 = 18, g = 3, phi = mod - 1;
int w[MAXN];
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] );
}
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 ) );
}
}
}
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 );
}
void PolyDer( int *ret, const int *coe, const int &n ) {
rep( i, 0, n - 2 ) ret[i] = Mul( coe[i + 1], i + 1 );
ret[n - 1] = 0;
}
void PolyInt( int *ret, const int *coe, const int &n ) {
per( i, n - 1, 0 ) ret[i + 1] = Mul( coe[i], Inv( i + 1 ) );
ret[0] = 0;
}
}
namespace PolyInv {
int P[MAXN], U[MAXN];
int A[MAXN], B[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 ) A[i] = B[i] = 0;
rep( i, 0, n - 1 ) A[i] = P[i];
rep( i, 0, h - 1 ) B[i] = U[i];
Basics :: DIF( A, L );
Basics :: DIF( B, L );
rep( i, 0, L - 1 ) A[i] = Mul( B[i], Sub( 2, Mul( B[i], A[i] ) ) );
Basics :: DIT( A, L );
rep( i, h, n - 1 ) U[i] = A[i];
}
void PolyInv( int *ret, const int *coe, const int &n ) {
rep( i, 0, n - 1 ) P[i] = coe[i];
Newton( n );
rep( i, 0, n - 1 ) ret[i] = U[i];
}
}
namespace PolySqrt {
int P[MAXN], U[MAXN];
int A[MAXN], B[MAXN];
void Newton( const int &n ) {
if( n == 1 ) {
U[0] = 3;
return ;
}
int h = ( n + 1 ) >> 1, L; Newton( h );
for( L = 1 ; L <= 2 * n - h - 1 ; L <<= 1 );
rep( i, 0, L - 1 ) A[i] = B[i] = 0;
rep( i, 0, n - 1 ) A[i] = P[i];
rep( i, 0, h - 1 ) B[i] = U[i];
PolyInv :: PolyInv( B, B, n );
Basics :: DIF( A, L );
Basics :: DIF( B, L );
rep( i, 0, L - 1 ) MulEq( A[i], B[i] );
Basics :: DIT( A, L );
rep( i, h, n - 1 ) U[i] = Mul( inv2, A[i] );
}
void PolySqrt( int *ret, const int *coe, const int &n ) {
rep( i, 0, n - 1 ) P[i] = coe[i];
Newton( n );
rep( i, 0, n - 1 ) ret[i] = U[i];
}
}
namespace PolyLn {
int P[MAXN];
int A[MAXN], B[MAXN];
void PolyLn( int *ret, const int *coe, const int &n ) {
rep( i, 0, n - 1 ) P[i] = coe[i];
int L; for( L = 1 ; L <= 2 * n - 2 ; L <<= 1 );
rep( i, 0, L - 1 ) A[i] = B[i] = 0;
Basics :: PolyDer( A, P, n );
PolyInv :: PolyInv( B, P, n );
Basics :: DIF( A, L );
Basics :: DIF( B, L );
rep( i, 0, L - 1 ) MulEq( A[i], B[i] );
Basics :: DIT( A, L );
Basics :: PolyInt( A, A, n );
rep( i, 0, n - 1 ) ret[i] = A[i];
}
}
namespace PolyExp {
int P[MAXN], U[MAXN];
int A[MAXN], B[MAXN];
void Newton( const int &n ) {
if( n == 1 ) {
U[0] = 1;
return ;
}
int h = ( n + 1 ) >> 1, L; Newton( h );
for( L = 1 ; L <= 2 * n - h - 1 ; L <<= 1 );
rep( i, 0, L - 1 ) A[i] = B[i] = 0;
rep( i, 0, h - 1 ) B[i] = U[i];
PolyLn :: PolyLn( A, B, n );
rep( i, 0, n - 1 ) A[i] = Sub( P[i], A[i] );
AddEq( A[0], 1 );
Basics :: DIF( A, L );
Basics :: DIF( B, L );
rep( i, 0, L - 1 ) MulEq( A[i], B[i] );
Basics :: DIT( A, L );
rep( i, h, n - 1 ) U[i] = A[i];
}
void PolyExp( int *ret, const int *coe, const int &n ) {
rep( i, 0, n - 1 ) P[i] = coe[i];
Newton( n );
rep( i, 0, n - 1 ) ret[i] = U[i];
}
}
namespace PolyPow {
int P[MAXN];
void PolyPow( int *ret, const int *coe, const int &n, const LL &k ) {
rep( i, 0, n - 1 ) P[i] = coe[i];
int c = P[0], inv = Inv( c );
rep( i, 0, n - 1 ) MulEq( P[i], inv );
PolyLn :: PolyLn( P, P, n );
rep( i, 0, n - 1 ) MulEq( P[i], k % mod );
PolyExp :: PolyExp( P, P, n );
c = Qkpow( c, k % ( mod - 1 ) );
rep( i, 0, n - 1 ) ret[i] = Mul( P[i], c );
}
}
int main() {
read( N ), read( K );
K ++, N ++, Basics :: NTTInit();
D[0] = 9, D[1] = mod - 2, D[2] = 1;
PolySqrt :: PolySqrt( D, D, K );
int L;
for( L = 1 ; L <= 2 * K - 2 ; L <<= 1 );
rep( i, 0, K - 1 ) B[i] = C[i] = D[i];
AddEq( B[0], 1 ), AddEq( B[1], 1 );
rep( i, 0, K - 1 ) MulEq( B[i], inv2 );
PolyPow :: PolyPow( B, B, K, N );
AddEq( C[0], 1 ), SubEq( C[1], 1 );
Basics :: DIF( B, L );
Basics :: DIF( C, L );
rep( i, 0, L - 1 ) MulEq( B[i], C[i] );
Basics :: DIT( B, L );
rep( i, 0, K - 1 ) AddEq( ans[i], B[i] );
rep( i, 0, L - 1 ) B[i] = C[i] = 0;
rep( i, 0, K - 1 ) B[i] = Mul( mod - 1, D[i] ), C[i] = D[i];
AddEq( B[0], 1 ), AddEq( B[1], 1 );
rep( i, 0, K - 1 ) MulEq( B[i], inv2 );
PolyPow :: PolyPow( B, B, K, N );
SubEq( C[0], 1 ), AddEq( C[1], 1 );
Basics :: DIF( B, L );
Basics :: DIF( C, L );
rep( i, 0, L - 1 ) MulEq( B[i], C[i] );
Basics :: DIT( B, L );
rep( i, 0, K - 1 ) AddEq( ans[i], B[i] );
rep( i, 0, K - 1 ) MulEq( D[i], 2 );
PolyInv :: PolyInv( D, D, K );
Basics :: DIF( ans, L );
Basics :: DIF( D, L );
rep( i, 0, L - 1 ) MulEq( ans[i], D[i] );
Basics :: DIT( ans, L );
rep( i, 0, K - 1 ) write( ans[i] ), putchar( " \n"[i == K - 1] );
return 0;
}