#line 1 "verify/fps/LC_inv_of_formal_power_series_2d.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/inv_of_formal_power_series_2d"
#line 2 "template/template.hpp"
#include <bits/stdc++.h>
using namespace std ;
#line 2 "template/macro.hpp"
#define rep(i, a, b) for (int i = (a); i < (int)(b); i++)
#define rrep(i, a, b) for (int i = (int)(b) - 1; i >= (a); i--)
#define ALL(v) (v).begin(), (v).end()
#define UNIQUE(v) sort(ALL(v)), (v).erase(unique(ALL(v)), (v).end())
#define SZ(v) (int)v.size()
#define MIN(v) *min_element(ALL(v))
#define MAX(v) *max_element(ALL(v))
#define LB(v, x) int(lower_bound(ALL(v), (x)) - (v).begin())
#define UB(v, x) int(upper_bound(ALL(v), (x)) - (v).begin())
#define YN(b) cout << ((b) ? "YES" : "NO") << "\n";
#define Yn(b) cout << ((b) ? "Yes" : "No") << "\n";
#define yn(b) cout << ((b) ? "yes" : "no") << "\n";
#line 6 "template/template.hpp"
#line 2 "template/util.hpp"
using uint = unsigned int ;
using ll = long long int ;
using ull = unsigned long long ;
using i128 = __int128_t ;
using u128 = __uint128_t ;
template < class T >
using priority_queue_asc = priority_queue < T , vector < T > , greater < T >> ;
template < class T , class S = T >
S SUM ( const vector < T >& a ) {
return accumulate ( ALL ( a ), S ( 0 ));
}
template < class T1 , class T2 >
inline bool chmin ( T1 & a , T2 b ) {
if ( a > b ) {
a = b ;
return true ;
}
return false ;
}
template < class T1 , class T2 >
inline bool chmax ( T1 & a , T2 b ) {
if ( a < b ) {
a = b ;
return true ;
}
return false ;
}
template < class T1 , class T2 >
inline bool chmin_opt ( optional < T1 >& a , T2 b ) {
if ( ! a || a > b ) {
a = b ;
return true ;
}
return false ;
}
template < class T1 , class T2 >
inline bool chmax_opt ( optional < T1 >& a , T2 b ) {
if ( ! a || a < b ) {
a = b ;
return true ;
}
return false ;
}
template < class T >
int popcnt ( T x ) {
return __builtin_popcountll ( x );
}
template < class T >
int topbit ( T x ) {
return ( x == 0 ? - 1 : 63 - __builtin_clzll ( x ));
}
template < class T >
int lowbit ( T x ) {
return ( x == 0 ? - 1 : __builtin_ctzll ( x ));
}
#line 8 "template/template.hpp"
#line 2 "template/inout.hpp"
struct Fast {
Fast () {
cin . tie ( nullptr );
ios_base :: sync_with_stdio ( false );
cout << fixed << setprecision ( 15 );
}
} fast ;
ostream & operator << ( ostream & os , __uint128_t x ) {
char buf [ 40 ];
size_t k = 0 ;
while ( x > 0 ) buf [ k ++ ] = ( char )( x % 10 + '0' ), x /= 10 ;
if ( k == 0 ) buf [ k ++ ] = '0' ;
while ( k ) os << buf [ -- k ];
return os ;
}
ostream & operator << ( ostream & os , __int128_t x ) {
return x < 0 ? ( os << '-' << ( __uint128_t )( - x )) : ( os << ( __uint128_t ) x );
}
template < class T , size_t N >
ostream & operator << ( ostream & os , const array < T , N >& a );
template < class T1 , class T2 >
istream & operator >> ( istream & is , pair < T1 , T2 >& p ) {
return is >> p . first >> p . second ;
}
template < class T1 , class T2 >
ostream & operator << ( ostream & os , const pair < T1 , T2 >& p ) {
return os << p . first << " " << p . second ;
}
template < class T >
istream & operator >> ( istream & is , vector < T >& a ) {
for ( auto & v : a ) is >> v ;
return is ;
}
template < class T >
ostream & operator << ( ostream & os , const vector < T >& a ) {
for ( auto it = a . begin (); it != a . end ();) {
os << * it ;
if ( ++ it != a . end ()) os << " " ;
}
return os ;
}
template < class T , size_t N >
ostream & operator << ( ostream & os , const array < T , N >& a ) {
for ( auto it = a . begin (); it != a . end ();) {
os << * it ;
if ( ++ it != a . end ()) os << " " ;
}
return os ;
}
template < class T >
ostream & operator << ( ostream & os , const set < T >& st ) {
os << "{" ;
for ( auto it = st . begin (); it != st . end ();) {
os << * it ;
if ( ++ it != st . end ()) os << "," ;
}
os << "}" ;
return os ;
}
template < class T1 , class T2 >
ostream & operator << ( ostream & os , const map < T1 , T2 >& mp ) {
os << "{" ;
for ( auto it = mp . begin (); it != mp . end ();) {
os << it -> first << ":" << it -> second ;
if ( ++ it != mp . end ()) os << "," ;
}
os << "}" ;
return os ;
}
void in () {}
template < typename T , class ... U >
void in ( T & t , U & ... u ) {
cin >> t ;
in ( u ...);
}
template < class ... T >
void in_zip ( int n , T & ... t ) {
assert ( n >= 0 && (( size ( t ) >= static_cast < size_t > ( n )) && ...));
for ( int i = 0 ; i < n ; i ++ ) in ( t [ i ]...);
}
void out () { cout << " \n " ; }
template < typename T , class ... U , char sep = ' ' >
void out ( const T & t , const U & ... u ) {
cout << t ;
if ( sizeof ...( u )) cout << sep ;
out ( u ...);
}
template < class T , class U >
void out_opt ( const optional < T >& opt , const U & fallback , ostream & os = cout ) {
if ( opt . has_value ())
os << opt . value ();
else
os << fallback ;
os << " \n " ;
}
template < class T , class U >
void out_opt ( const vector < optional < T >>& vec , const U & fallback , ostream & os = cout ) {
for ( auto it = vec . begin (); it != vec . end ();) {
if (( * it ). has_value ())
os << ( * it ). value ();
else
os << fallback ;
if ( ++ it != vec . end ()) os << " " ;
}
os << " \n " ;
}
namespace IO {
template < class T , class ... U >
T read ( U && ... u ) {
T t = T ( forward < U > ( u )...);
in ( t );
return t ;
}
namespace Graph {
vector < vector < int >> unweighted ( int n , int m , bool directed = false , int offset = 1 ) {
vector < vector < int >> g ( n );
for ( int i = 0 ; i < m ; i ++ ) {
int u , v ;
cin >> u >> v ;
u -= offset , v -= offset ;
g [ u ]. push_back ( v );
if ( ! directed ) g [ v ]. push_back ( u );
}
return g ;
}
template < class T >
vector < vector < pair < int , T >>> weighted ( int n , int m , bool directed = false , int offset = 1 ) {
vector < vector < pair < int , T >>> g ( n );
for ( int i = 0 ; i < m ; i ++ ) {
int u , v ;
T w ;
cin >> u >> v >> w ;
u -= offset , v -= offset ;
g [ u ]. push_back ({ v , w });
if ( ! directed ) g [ v ]. push_back ({ u , w });
}
return g ;
}
} // namespace Graph
namespace Tree {
vector < vector < int >> unweighted ( int n , bool directed = false , int offset = 1 ) {
return Graph :: unweighted ( n , n - 1 , directed , offset );
}
template < class T >
vector < vector < pair < int , T >>> weighted ( int n , bool directed = false , int offset = 1 ) {
return Graph :: weighted < T > ( n , n - 1 , directed , offset );
}
vector < vector < int >> rooted ( int n , bool to_root = true , bool to_leaf = true , int offset = 1 ) {
vector < vector < int >> g ( n );
for ( int i = 1 ; i < n ; i ++ ) {
int p ;
cin >> p ;
p -= offset ;
if ( to_root ) g [ i ]. push_back ( p );
if ( to_leaf ) g [ p ]. push_back ( i );
}
return g ;
}
} // namespace Tree
} // namespace IO
#line 10 "template/template.hpp"
#line 2 "template/debug.hpp"
#ifdef LOCAL
#define debug 1
#define show(...) _show(0, #__VA_ARGS__, __VA_ARGS__)
#else
#define debug 0
#define show(...) true
#endif
template < class T >
void _show ( int , T ) {
cerr << '\n' ;
}
template < class T1 , class T2 , class ... T3 >
void _show ( int i , const T1 & a , const T2 & b , const T3 & ... c ) {
for (; a [ i ] != ',' && a [ i ] != '\0' ; i ++ ) cerr << a [ i ];
cerr << ":" << b << " " ;
_show ( i + 1 , a , c ...);
}
#line 2 "math/util.hpp"
namespace Math {
template < class T >
T safe_mod ( T a , T b ) {
assert ( b != 0 );
if ( b < 0 ) a = - a , b = - b ;
a %= b ;
return a >= 0 ? a : a + b ;
}
template < class T >
T floor ( T a , T b ) {
assert ( b != 0 );
if ( b < 0 ) a = - a , b = - b ;
return a >= 0 ? a / b : ( a + 1 ) / b - 1 ;
}
template < class T >
T ceil ( T a , T b ) {
assert ( b != 0 );
if ( b < 0 ) a = - a , b = - b ;
return a > 0 ? ( a - 1 ) / b + 1 : a / b ;
}
long long isqrt ( long long n ) {
if ( n <= 0 ) return 0 ;
long long x = sqrt ( n );
while (( __int128 )( x + 1 ) * ( x + 1 ) <= n ) x ++ ;
while (( __int128 ) x * x > n ) x -- ;
return x ;
}
long long floor_root ( long long n , int k ) {
assert ( n >= 0 );
if ( n == 0 ) return 0 ;
assert ( k >= 1 );
if ( k == 1 ) return n ;
if ( k > 64 ) return 1 ;
long long x = round ( pow (( long double ) n , 1.0 L / k ));
auto check = [ & ]( long long a ) {
if ( a <= 0 ) return true ;
__int128_t p = 1 ;
for ( int i = 0 ; i < k ; ++ i )
if (( p *= a ) > n ) return false ;
return true ;
};
while ( check ( x + 1 )) x ++ ;
while ( ! check ( x )) x -- ;
return x ;
}
unsigned long long floor_root_unsigned ( unsigned long long n , int k ) {
assert ( k >= 1 );
if ( n <= 1 || k == 1 ) return n ;
if ( k >= 64 ) return 1 ;
int bits = ( 64 + k - 1 ) / k ;
unsigned long long ok = 1 , ng = min ( n , 1ULL << bits );
auto check = [ & ]( unsigned long long a ) {
__uint128_t p = 1 ;
for ( int i = 0 ; i < k ; i ++ ) {
p *= a ;
if ( p > n ) return false ;
}
return true ;
};
while ( ok + 1 < ng ) {
unsigned long long mid = ok + ( ng - ok ) / 2 ;
( check ( mid ) ? ok : ng ) = mid ;
}
return ok ;
}
// return g=gcd(a,b)
// a*x+b*y=g
// - b!=0 -> 0<=x<|b|/g
// - b=0 -> ax=g
template < class T >
T ext_gcd ( T a , T b , T & x , T & y ) {
T a0 = a , b0 = b ;
bool sgn_a = a < 0 , sgn_b = b < 0 ;
if ( sgn_a ) a = - a ;
if ( sgn_b ) b = - b ;
if ( b == 0 ) {
x = sgn_a ? - 1 : 1 ;
y = 0 ;
return a ;
}
T x00 = 1 , x01 = 0 , x10 = 0 , x11 = 1 ;
while ( b != 0 ) {
T q = a / b , r = a - b * q ;
x00 -= q * x01 ;
x10 -= q * x11 ;
swap ( x00 , x01 );
swap ( x10 , x11 );
a = b , b = r ;
}
x = x00 , y = x10 ;
if ( sgn_a ) x = - x ;
if ( sgn_b ) y = - y ;
if ( b0 != 0 ) {
a0 /= a , b0 /= a ;
if ( b0 < 0 ) a0 = - a0 , b0 = - b0 ;
T q = x >= 0 ? x / b0 : ( x + 1 ) / b0 - 1 ;
x -= b0 * q ;
y += a0 * q ;
}
return a ;
}
constexpr long long inv_mod ( long long x , long long m ) {
x %= m ;
if ( x < 0 ) x += m ;
long long a = m , b = x ;
long long y0 = 0 , y1 = 1 ;
while ( b > 0 ) {
long long q = a / b ;
swap ( a -= q * b , b );
swap ( y0 -= q * y1 , y1 );
}
if ( y0 < 0 ) y0 += m / a ;
return y0 ;
}
long long pow_mod ( long long x , long long n , long long m ) {
if ( m == 1 ) return 0 ;
x = ( x % m + m ) % m ;
long long y = 1 ;
while ( n ) {
if ( n & 1 ) y = y * x % m ;
x = x * x % m ;
n >>= 1 ;
}
return y ;
}
constexpr long long pow_mod_constexpr ( long long x , long long n , int m ) {
if ( m == 1 ) return 0 ;
unsigned int _m = ( unsigned int )( m );
unsigned long long r = 1 ;
unsigned long long y = x % m ;
if ( y >= m ) y += m ;
while ( n ) {
if ( n & 1 ) r = ( r * y ) % _m ;
y = ( y * y ) % _m ;
n >>= 1 ;
}
return r ;
}
constexpr bool is_prime_constexpr ( int n ) {
if ( n <= 1 ) return false ;
if ( n == 2 || n == 7 || n == 61 ) return true ;
if ( n % 2 == 0 ) return false ;
long long d = n - 1 ;
while ( d % 2 == 0 ) d /= 2 ;
constexpr long long bases [ 3 ] = { 2 , 7 , 61 };
for ( long long a : bases ) {
long long t = d ;
long long y = pow_mod_constexpr ( a , t , n );
while ( t != n - 1 && y != 1 && y != n - 1 ) {
y = y * y % n ;
t <<= 1 ;
}
if ( y != n - 1 && t % 2 == 0 ) {
return false ;
}
}
return true ;
}
template < int n >
constexpr bool is_prime = is_prime_constexpr ( n );
}; // namespace Math
#line 3 "modint/modint.hpp"
template < unsigned int m = 998244353 >
struct ModInt {
using mint = ModInt ;
static constexpr unsigned int get_mod () { return m ; }
static mint raw ( int v ) {
mint x ;
x . _v = v ;
return x ;
}
ModInt () : _v ( 0 ) {}
ModInt ( int64_t v ) {
long long x = ( long long )( v % ( long long )( umod ()));
if ( x < 0 ) x += umod ();
_v = ( unsigned int )( x );
}
unsigned int val () const { return _v ; }
mint & operator ++ () {
_v ++ ;
if ( _v == umod ()) _v = 0 ;
return * this ;
}
mint & operator -- () {
if ( _v == 0 ) _v = umod ();
_v -- ;
return * this ;
}
mint operator ++ ( int ) {
mint result = * this ;
++* this ;
return result ;
}
mint operator -- ( int ) {
mint result = * this ;
--* this ;
return result ;
}
mint & operator += ( const mint & rhs ) {
_v += rhs . _v ;
if ( _v >= umod ()) _v -= umod ();
return * this ;
}
mint & operator -= ( const mint & rhs ) {
_v -= rhs . _v ;
if ( _v >= umod ()) _v += umod ();
return * this ;
}
mint & operator *= ( const mint & rhs ) {
unsigned long long z = _v ;
z *= rhs . _v ;
_v = ( unsigned int )( z % umod ());
return * this ;
}
mint & operator /= ( const mint & rhs ) { return * this *= rhs . inv (); }
mint operator + () const { return * this ; }
mint operator - () const { return mint () - * this ; }
mint pow ( long long n ) const {
assert ( 0 <= n );
mint x = * this , r = 1 ;
while ( n ) {
if ( n & 1 ) r *= x ;
x *= x ;
n >>= 1 ;
}
return r ;
}
mint inv () const {
if ( is_prime ) {
assert ( _v );
return pow ( umod () - 2 );
} else {
auto inv = Math :: inv_mod ( _v , umod ());
return raw ( inv );
}
}
friend mint operator + ( const mint & lhs , const mint & rhs ) { return mint ( lhs ) += rhs ; }
friend mint operator - ( const mint & lhs , const mint & rhs ) { return mint ( lhs ) -= rhs ; }
friend mint operator * ( const mint & lhs , const mint & rhs ) { return mint ( lhs ) *= rhs ; }
friend mint operator / ( const mint & lhs , const mint & rhs ) { return mint ( lhs ) /= rhs ; }
friend bool operator == ( const mint & lhs , const mint & rhs ) { return lhs . _v == rhs . _v ; }
friend bool operator != ( const mint & lhs , const mint & rhs ) { return lhs . _v != rhs . _v ; }
friend istream & operator >> ( istream & is , mint & x ) {
int64_t v ;
is >> v ;
x = mint ( v );
return is ;
}
friend ostream & operator << ( ostream & os , const mint & x ) { return os << x . val (); }
private:
unsigned int _v ;
static constexpr unsigned int umod () { return m ; }
static constexpr bool is_prime = Math :: is_prime < m > ;
};
using ModInt998244353 = ModInt < 998244353 > ;
using ModInt1000000007 = ModInt < 1000000007 > ;
#line 5 "verify/fps/LC_inv_of_formal_power_series_2d.test.cpp"
using mint = ModInt < 998244353 > ;
#line 2 "fps/fps-2d-ntt-friendly.hpp"
#line 2 "fft/ntt.hpp"
template < class mint >
struct NTT {
static constexpr unsigned int mod = mint :: get_mod ();
static constexpr unsigned long long pow_constexpr ( unsigned long long x , unsigned long long n , unsigned long long m ) {
unsigned long long y = 1 ;
while ( n ) {
if ( n & 1 ) y = y * x % m ;
x = x * x % m ;
n >>= 1 ;
}
return y ;
}
static constexpr unsigned int get_g () {
unsigned long long x = 2 ;
while ( pow_constexpr ( x , ( mod - 1 ) >> 1 , mod ) == 1 ) x += 1 ;
return x ;
}
static constexpr unsigned int g = get_g ();
static constexpr int rank2 = __builtin_ctzll ( mod - 1 );
array < mint , rank2 + 1 > root ;
array < mint , rank2 + 1 > iroot ;
array < mint , max ( 0 , rank2 - 2 + 1 ) > rate2 ;
array < mint , max ( 0 , rank2 - 2 + 1 ) > irate2 ;
array < mint , max ( 0 , rank2 - 3 + 1 ) > rate3 ;
array < mint , max ( 0 , rank2 - 3 + 1 ) > irate3 ;
NTT () {
root [ rank2 ] = mint ( g ). pow (( mod - 1 ) >> rank2 );
iroot [ rank2 ] = root [ rank2 ]. inv ();
for ( int i = rank2 - 1 ; i >= 0 ; i -- ) {
root [ i ] = root [ i + 1 ] * root [ i + 1 ];
iroot [ i ] = iroot [ i + 1 ] * iroot [ i + 1 ];
}
{
mint prod = 1 , iprod = 1 ;
for ( int i = 0 ; i <= rank2 - 2 ; i ++ ) {
rate2 [ i ] = root [ i + 2 ] * prod ;
irate2 [ i ] = iroot [ i + 2 ] * iprod ;
prod *= iroot [ i + 2 ];
iprod *= root [ i + 2 ];
}
}
{
mint prod = 1 , iprod = 1 ;
for ( int i = 0 ; i <= rank2 - 3 ; i ++ ) {
rate3 [ i ] = root [ i + 3 ] * prod ;
irate3 [ i ] = iroot [ i + 3 ] * iprod ;
prod *= iroot [ i + 3 ];
iprod *= root [ i + 3 ];
}
}
}
void ntt ( vector < mint >& a ) {
int n = int ( a . size ());
int h = __builtin_ctzll (( unsigned int ) n );
assert ( h <= rank2 );
a . resize ( 1 << h );
int len = 0 ; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed
while ( len < h ) {
if ( h - len == 1 ) {
int p = 1 << ( h - len - 1 );
mint rot = 1 ;
for ( int s = 0 ; s < ( 1 << len ); s ++ ) {
int offset = s << ( h - len );
for ( int i = 0 ; i < p ; i ++ ) {
auto l = a [ i + offset ];
auto r = a [ i + offset + p ] * rot ;
a [ i + offset ] = l + r ;
a [ i + offset + p ] = l - r ;
}
if ( s + 1 != ( 1 << len )) rot *= rate2 [ __builtin_ctzll ( ~ ( unsigned int )( s ))];
}
len ++ ;
} else {
// 4-base
int p = 1 << ( h - len - 2 );
mint rot = 1 , imag = root [ 2 ];
for ( int s = 0 ; s < ( 1 << len ); s ++ ) {
mint rot2 = rot * rot ;
mint rot3 = rot2 * rot ;
int offset = s << ( h - len );
for ( int i = 0 ; i < p ; i ++ ) {
auto mod2 = 1ULL * mint :: get_mod () * mint :: get_mod ();
auto a0 = 1ULL * a [ i + offset ]. val ();
auto a1 = 1ULL * a [ i + offset + p ]. val () * rot . val ();
auto a2 = 1ULL * a [ i + offset + 2 * p ]. val () * rot2 . val ();
auto a3 = 1ULL * a [ i + offset + 3 * p ]. val () * rot3 . val ();
auto a1na3imag = 1ULL * mint ( a1 + mod2 - a3 ). val () * imag . val ();
auto na2 = mod2 - a2 ;
a [ i + offset ] = a0 + a2 + a1 + a3 ;
a [ i + offset + 1 * p ] = a0 + a2 + ( 2 * mod2 - ( a1 + a3 ));
a [ i + offset + 2 * p ] = a0 + na2 + a1na3imag ;
a [ i + offset + 3 * p ] = a0 + na2 + ( mod2 - a1na3imag );
}
if ( s + 1 != ( 1 << len )) rot *= rate3 [ __builtin_ctzll ( ~ ( unsigned int )( s ))];
}
len += 2 ;
}
}
}
void intt ( vector < mint >& a ) {
int n = int ( a . size ());
int h = __builtin_ctzll (( unsigned int ) n );
assert ( h <= rank2 );
a . resize ( 1 << h );
int len = h ; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed
while ( len ) {
if ( len == 1 ) {
int p = 1 << ( h - len );
mint irot = 1 ;
for ( int s = 0 ; s < ( 1 << ( len - 1 )); s ++ ) {
int offset = s << ( h - len + 1 );
for ( int i = 0 ; i < p ; i ++ ) {
auto l = a [ i + offset ];
auto r = a [ i + offset + p ];
a [ i + offset ] = l + r ;
a [ i + offset + p ] = ( unsigned long long )( mint :: get_mod () + l . val () - r . val ()) * irot . val ();
}
if ( s + 1 != ( 1 << ( len - 1 ))) irot *= irate2 [ __builtin_ctzll ( ~ ( unsigned int )( s ))];
}
len -- ;
} else {
// 4-base
int p = 1 << ( h - len );
mint irot = 1 , iimag = iroot [ 2 ];
for ( int s = 0 ; s < ( 1 << ( len - 2 )); s ++ ) {
mint irot2 = irot * irot ;
mint irot3 = irot2 * irot ;
int offset = s << ( h - len + 2 );
for ( int i = 0 ; i < p ; i ++ ) {
auto a0 = 1ULL * a [ i + offset + 0 * p ]. val ();
auto a1 = 1ULL * a [ i + offset + 1 * p ]. val ();
auto a2 = 1ULL * a [ i + offset + 2 * p ]. val ();
auto a3 = 1ULL * a [ i + offset + 3 * p ]. val ();
auto a2na3iimag = 1ULL * mint (( mint :: get_mod () + a2 - a3 ) * iimag . val ()). val ();
a [ i + offset ] = a0 + a1 + a2 + a3 ;
a [ i + offset + 1 * p ] = ( a0 + ( mint :: get_mod () - a1 ) + a2na3iimag ) * irot . val ();
a [ i + offset + 2 * p ] = ( a0 + a1 + ( mint :: get_mod () - a2 ) + ( mint :: get_mod () - a3 )) * irot2 . val ();
a [ i + offset + 3 * p ] = ( a0 + ( mint :: get_mod () - a1 ) + ( mint :: get_mod () - a2na3iimag )) * irot3 . val ();
}
if ( s + 1 != ( 1 << ( len - 2 ))) irot *= irate3 [ __builtin_ctzll ( ~ ( unsigned int )( s ))];
}
len -= 2 ;
}
}
mint e = mint ( n ). inv ();
for ( auto & x : a ) x *= e ;
}
vector < mint > multiply ( const vector < mint >& a , const vector < mint >& b ) {
if ( a . empty () || b . empty ()) return vector < mint > ();
int n = a . size (), m = b . size ();
int sz = n + m - 1 ;
if ( n <= 30 || m <= 30 ) {
if ( n > 30 ) return multiply ( b , a );
vector < mint > res ( sz );
for ( int i = 0 ; i < n ; i ++ )
for ( int j = 0 ; j < m ; j ++ ) res [ i + j ] += a [ i ] * b [ j ];
return res ;
}
int sz1 = 1 ;
while ( sz1 < sz ) sz1 <<= 1 ;
vector < mint > res ( sz1 );
for ( int i = 0 ; i < n ; i ++ ) res [ i ] = a [ i ];
ntt ( res );
if ( a == b )
for ( int i = 0 ; i < sz1 ; i ++ ) res [ i ] *= res [ i ];
else {
vector < mint > c ( sz1 );
for ( int i = 0 ; i < m ; i ++ ) c [ i ] = b [ i ];
ntt ( c );
for ( int i = 0 ; i < sz1 ; i ++ ) res [ i ] *= c [ i ];
}
intt ( res );
res . resize ( sz );
return res ;
}
// c[i]=sum[j]a[j]b[i+j]
vector < mint > middle_product ( const vector < mint >& a , const vector < mint >& b ) {
if ( b . empty () || a . size () > b . size ()) return {};
int n = a . size (), m = b . size ();
int sz = m - n + 1 ;
if ( n <= 30 || sz <= 30 ) {
vector < mint > res ( sz );
for ( int i = 0 ; i < sz ; i ++ )
for ( int j = 0 ; j < n ; j ++ ) res [ i ] += a [ j ] * b [ i + j ];
return res ;
}
int sz1 = 1 ;
while ( sz1 < m ) sz1 <<= 1 ;
vector < mint > res ( sz1 ), b2 ( sz1 );
reverse_copy ( a . begin (), a . end (), res . begin ());
copy ( b . begin (), b . end (), b2 . begin ());
ntt ( res );
ntt ( b2 );
for ( int i = 0 ; i < res . size (); i ++ ) res [ i ] *= b2 [ i ];
intt ( res );
res . resize ( m );
res . erase ( res . begin (), res . begin () + n - 1 );
return res ;
}
void ntt_doubling ( vector < mint >& a ) {
int n = ( int ) a . size ();
auto b = a ;
intt ( b );
mint r = 1 , zeta = mint ( g ). pow (( mint :: get_mod () - 1 ) / ( n << 1 ));
for ( int i = 0 ; i < n ; i ++ ) b [ i ] *= r , r *= zeta ;
ntt ( b );
copy ( b . begin (), b . end (), back_inserter ( a ));
}
};
/**
* @brief NTT (数論変換)
* @docs docs/fft/ntt.md
*/
#line 2 "fps/fps-2d.hpp"
/**
* @brief 二変数形式的冪級数
* @docs docs/fps/fps-2d.md
*/
template < class mint >
struct FormalPowerSeries2D : vector < vector < mint >> {
using Base = vector < vector < mint >> ;
using FPS2D = FormalPowerSeries2D ;
using Base :: Base ;
FormalPowerSeries2D ( int n , int m ) : Base ( n , vector < mint > ( m )) {}
FormalPowerSeries2D ( const Base & r ) : Base ( r ) {}
FormalPowerSeries2D ( Base && r ) : Base ( std :: move ( r )) {}
int height () const { return ( int ) this -> size (); }
int width () const { return this -> empty () ? 0 : ( int )( * this )[ 0 ]. size (); }
void resize ( int n , int m ) {
assert ( n >= 0 && m >= 0 );
if ( n == 0 || m == 0 ) {
this -> clear ();
return ;
}
Base :: resize ( n );
for ( auto & row : * this ) row . resize ( m );
}
FPS2D & operator = ( const Base & r ) {
Base :: operator = ( r );
return * this ;
}
FPS2D & operator = ( Base && r ) {
Base :: operator = ( std :: move ( r ));
return * this ;
}
FPS2D & operator += ( const FPS2D & r ) {
int n = max ( height (), r . height ()), m = max ( width (), r . width ());
resize ( n , m );
for ( int i = 0 ; i < r . height (); i ++ ) {
assert (( int ) r [ i ]. size () == r . width ());
for ( int j = 0 ; j < r . width (); j ++ ) ( * this )[ i ][ j ] += r [ i ][ j ];
}
return * this ;
}
FPS2D & operator += ( const mint & r ) {
if ( this -> empty ()) resize ( 1 , 1 );
( * this )[ 0 ][ 0 ] += r ;
return * this ;
}
FPS2D & operator -= ( const FPS2D & r ) {
int n = max ( height (), r . height ()), m = max ( width (), r . width ());
resize ( n , m );
for ( int i = 0 ; i < r . height (); i ++ ) {
assert (( int ) r [ i ]. size () == r . width ());
for ( int j = 0 ; j < r . width (); j ++ ) ( * this )[ i ][ j ] -= r [ i ][ j ];
}
return * this ;
}
FPS2D & operator -= ( const mint & r ) {
if ( this -> empty ()) resize ( 1 , 1 );
( * this )[ 0 ][ 0 ] -= r ;
return * this ;
}
FPS2D & operator *= ( const mint & r ) {
for ( auto & row : * this )
for ( auto & x : row ) x *= r ;
return * this ;
}
FPS2D & operator /= ( const mint & r ) { return * this *= r . inv (); }
FPS2D & operator /= ( const FPS2D & r ) {
int n = height (), m = width ();
assert ( n > 0 && m > 0 );
return * this = (( * this ) * r . inv ( n , m )). pre ( n , m );
}
FPS2D operator + ( const FPS2D & r ) const { return FPS2D ( * this ) += r ; }
FPS2D operator + ( const mint & r ) const { return FPS2D ( * this ) += r ; }
FPS2D operator - ( const FPS2D & r ) const { return FPS2D ( * this ) -= r ; }
FPS2D operator - ( const mint & r ) const { return FPS2D ( * this ) -= r ; }
FPS2D operator * ( const FPS2D & r ) const { return FPS2D ( * this ) *= r ; }
FPS2D operator * ( const mint & r ) const { return FPS2D ( * this ) *= r ; }
FPS2D operator / ( const FPS2D & r ) const { return FPS2D ( * this ) /= r ; }
FPS2D operator / ( const mint & r ) const { return FPS2D ( * this ) /= r ; }
FPS2D operator - () const {
FPS2D ret ( * this );
for ( auto & row : ret )
for ( auto & x : row ) x = - x ;
return ret ;
}
friend FPS2D operator + ( const mint & l , const FPS2D & r ) { return r + l ; }
friend FPS2D operator - ( const mint & l , const FPS2D & r ) { return - r + l ; }
friend FPS2D operator * ( const mint & l , const FPS2D & r ) { return r * l ; }
FPS2D shift ( int di , int dj ) const {
assert ( di >= 0 && dj >= 0 );
if ( this -> empty ()) return {};
int m = width ();
FPS2D ret ( * this );
for ( auto & row : ret ) {
assert (( int ) row . size () == m );
row . insert ( row . begin (), dj , mint ( 0 ));
}
ret . insert ( ret . begin (), di , vector < mint > ( m + dj ));
return ret ;
}
FPS2D pre ( int n , int m ) const {
assert ( n >= 0 && m >= 0 );
n = min ( n , height ()), m = min ( m , width ());
if ( n == 0 || m == 0 ) return {};
FPS2D ret ( n , m );
for ( int i = 0 ; i < n ; i ++ ) {
assert (( int )( * this )[ i ]. size () == width ());
copy_n (( * this )[ i ]. begin (), m , ret [ i ]. begin ());
}
return ret ;
}
FPS2D log ( int n = - 1 , int m = - 1 ) const {
assert ( ! this -> empty () && width () > 0 && ( * this )[ 0 ][ 0 ] == mint ( 1 ));
if ( n == - 1 ) n = height ();
if ( m == - 1 ) m = width ();
assert ( n >= 0 && m >= 0 );
if ( n == 0 || m == 0 ) return {};
FPS2D d = this -> pre ( n , m );
d . resize ( n , m );
for ( int i = 0 ; i < n ; i ++ )
for ( int j = 0 ; j < m ; j ++ ) d [ i ][ j ] *= mint ( i + j );
FPS2D ret = ( d * this -> inv ( n , m )). pre ( n , m );
ret . resize ( n , m );
ret [ 0 ][ 0 ] = mint ( 0 );
for ( int i = 0 ; i < n ; i ++ )
for ( int j = 0 ; j < m ; j ++ )
if ( i + j > 0 ) ret [ i ][ j ] /= mint ( i + j );
return ret ;
}
static void * convolution_ptr ;
static void set_convolution ();
static vector < mint > convolution ( const vector < mint >& a , const vector < mint >& b );
FPS2D & operator *= ( const FPS2D & r );
FPS2D inv ( int n = - 1 , int m = - 1 ) const ;
FPS2D exp ( int n = - 1 , int m = - 1 ) const ;
};
template < class mint >
void * FormalPowerSeries2D < mint >:: convolution_ptr = nullptr ;
#line 5 "fps/fps-2d-ntt-friendly.hpp"
template < class mint >
void FormalPowerSeries2D < mint >:: set_convolution () {
if ( ! convolution_ptr ) convolution_ptr = new NTT < mint > ;
}
template < class mint >
vector < mint > FormalPowerSeries2D < mint >:: convolution ( const vector < mint >& a , const vector < mint >& b ) {
set_convolution ();
return static_cast < NTT < mint >*> ( convolution_ptr ) -> multiply ( a , b );
}
template < class mint >
FormalPowerSeries2D < mint >& FormalPowerSeries2D < mint >:: operator *= ( const FPS2D & r ) {
if ( this -> empty () || r . empty ()) {
this -> clear ();
return * this ;
}
int n = height (), m = width (), rn = r . height (), rm = r . width ();
for ( const auto & row : * this ) assert (( int ) row . size () == m );
for ( const auto & row : r ) assert (( int ) row . size () == rm );
int on = n + rn - 1 , om = m + rm - 1 ;
vector < mint > a (( n - 1 ) * om + m ), b (( rn - 1 ) * om + rm );
for ( int i = 0 ; i < n ; i ++ ) copy (( * this )[ i ]. begin (), ( * this )[ i ]. end (), a . begin () + i * om );
for ( int i = 0 ; i < rn ; i ++ ) copy ( r [ i ]. begin (), r [ i ]. end (), b . begin () + i * om );
vector < mint > c = convolution ( a , b );
FPS2D ret ( on , om );
for ( int i = 0 ; i < on ; i ++ ) copy_n ( c . begin () + i * om , om , ret [ i ]. begin ());
return * this = std :: move ( ret );
}
template < class mint >
FormalPowerSeries2D < mint > FormalPowerSeries2D < mint >:: inv ( int n , int m ) const {
assert ( ! this -> empty () && width () > 0 && ( * this )[ 0 ][ 0 ] != mint ( 0 ));
if ( n == - 1 ) n = height ();
if ( m == - 1 ) m = width ();
assert ( n >= 0 && m >= 0 );
if ( n == 0 || m == 0 ) return {};
FPS2D ret {{ mint ( 1 ) / ( * this )[ 0 ][ 0 ]}};
int deg = 1 ;
while ( deg < n + m - 1 ) {
deg <<= 1 ;
int nn = min ( n , deg ), mm = min ( m , deg );
ret . resize ( nn , mm );
FPS2D c = ( this -> pre ( nn , mm ) * ret ). pre ( nn , mm );
c = - c ;
c [ 0 ][ 0 ] += mint ( 2 );
ret = ( ret * c ). pre ( nn , mm );
}
ret . resize ( n , m );
return ret ;
}
template < class mint >
FormalPowerSeries2D < mint > FormalPowerSeries2D < mint >:: exp ( int n , int m ) const {
assert ( ! this -> empty () && width () > 0 && ( * this )[ 0 ][ 0 ] == mint ( 0 ));
if ( n == - 1 ) n = height ();
if ( m == - 1 ) m = width ();
assert ( n >= 0 && m >= 0 );
if ( n == 0 || m == 0 ) return {};
FPS2D ret {{ mint ( 1 )}};
int deg = 1 ;
while ( deg < n + m - 1 ) {
deg <<= 1 ;
int nn = min ( n , deg ), mm = min ( m , deg );
ret . resize ( nn , mm );
FPS2D c = this -> pre ( nn , mm ) - ret . log ( nn , mm ) + mint ( 1 );
ret = ( ret * c ). pre ( nn , mm );
}
ret . resize ( n , m );
return ret ;
}
#line 7 "verify/fps/LC_inv_of_formal_power_series_2d.test.cpp"
using fps2d = FormalPowerSeries2D < mint > ;
int main () {
int n , m ;
in ( n , m );
fps2d a ( n , m );
in ( a );
out ( a . inv ());
} 
verify/fps/LC_inv_of_formal_power_series_2d.test.cpp
math/util.hpp