結果
| 問題 | No.118 門松列(2) |
| ユーザー |
 Wiiiiam104
|
| 提出日時 | 2024-06-19 15:07:13 |
| 言語 | C++23 (gcc 12.3.0 + boost 1.83.0) |
| 結果 |
AC
|
| 実行時間 | 9 ms / 5,000 ms |
| コード長 | 39,457 bytes |
| コンパイル時間 | 7,511 ms |
| コンパイル使用メモリ | 328,828 KB |
| 実行使用メモリ | 6,944 KB |
| 最終ジャッジ日時 | 2024-06-19 15:07:22 |
| 合計ジャッジ時間 | 7,471 ms |
|
ジャッジサーバーID （参考情報） |
judge5 / judge2 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 3 ms
6,816 KB |
| testcase_01 | AC | 9 ms
6,812 KB |
| testcase_02 | AC | 9 ms
6,944 KB |
| testcase_03 | AC | 9 ms
6,940 KB |
| testcase_04 | AC | 9 ms
6,940 KB |
| testcase_05 | AC | 9 ms
6,944 KB |
| testcase_06 | AC | 3 ms
6,940 KB |
| testcase_07 | AC | 3 ms
6,940 KB |
| testcase_08 | AC | 3 ms
6,940 KB |
| testcase_09 | AC | 9 ms
6,940 KB |
| testcase_10 | AC | 4 ms
6,940 KB |
| testcase_11 | AC | 5 ms
6,944 KB |
| testcase_12 | AC | 4 ms
6,940 KB |
| testcase_13 | AC | 3 ms
6,944 KB |
| testcase_14 | AC | 8 ms
6,944 KB |
| testcase_15 | AC | 3 ms
6,940 KB |
| testcase_16 | AC | 6 ms
6,940 KB |
| testcase_17 | AC | 3 ms
6,944 KB |
| testcase_18 | AC | 4 ms
6,944 KB |
| testcase_19 | AC | 6 ms
6,940 KB |
| testcase_20 | AC | 6 ms
6,944 KB |
| testcase_21 | AC | 8 ms
6,940 KB |
| testcase_22 | AC | 4 ms
6,940 KB |
| testcase_23 | AC | 7 ms
6,940 KB |
| testcase_24 | AC | 5 ms
6,944 KB |
| testcase_25 | AC | 6 ms
6,940 KB |
ソースコード
#if __INCLUDE_LEVEL__
#include <bits/stdc++.h>
#include <atcoder/all>
//#define int ll
using mint=atcoder::static_modint<1000000007>;
//using mint=modint; atcoder::modint::set_mod(int m);
#endif
#if !__INCLUDE_LEVEL__
#include __FILE__
signed main(){
//cout<<fixed<<setprecision(10);
cin.tie(0)->sync_with_stdio(0);
//int start_clock=clock();
//srand((unsigned)time(NULL));
//your code here!
int n;cin>>n;
vi a(n);cin>>a;
vi t(109);
rep(i,n) t[a[i]]++;
FPSFFT<mint> f={1};
rep(i,108) f*=FPSFFT<mint>{1,t[i]};
ans(f[3].val());
}
#else //template
using namespace std;
using ull = unsigned long long;
using ll = long long;
using vi = vector<int>;
using vl = vector<long>;
using vll = vector<long long>;
using vb = vector<bool>;
using vf = vector<float>;
using vd = vector<double>;
using vvi = vector<vi>;
using vvl = vector<vl>;
using vvll = vector<vll>;
using vvf = vector<vf>;
using vvd = vector<vd>;
using vvvi = vector<vvi>;
using vvvd = vector<vvd>;
using vs = vector<string>;
using pii = pair<int,int>;
using vpii = vector<pii>;
using vm=vector<mint>;
using vvm=vector<vm>;
using vvvm=vector<vvm>;
#define ERROR(...) {cerr<<"ERROR @ LINE "<<__LINE__<<endl; exit(0);}
#define REP_OVERLOAD(_1,_2,_3,NAME,...) NAME
#define REP2(i,n) for(int i=0;(i)<(int)(n);(i)++)
#define REP3(i,a,n) for(int i=(int)(a);(i)<(int)(n);(i)++)
#define rep(...) REP_OVERLOAD(__VA_ARGS__, REP3, REP2, ERROR)(__VA_ARGS__)
#define REPD_OVERLOAD(_1,_2,_3,NAME,...) NAME
#define REPD2(i,n) for(int i=(int)(n)-1;(i)>=0;(i)--)
#define REPD3(i,a,n) for(int i=(int)(n)-1;(i)>=(int)(a);(i)--)
#define repd(...) REPD_OVERLOAD(__VA_ARGS__, REPD3, REPD2, ERROR)(__VA_ARGS__)
template<typename T> inline T chmax(T &a,const T &b){if(a<b)a=b;return a;}
template<typename T> inline T chmin(T &a,const T &b){if(b<a)a=b;return a;}
#define pb push_back
#define F first
#define S second
// cf. https://satiseni.hateblo.jp/entry/2019/06/22/122116
struct{
constexpr operator int(){return 1<<30;}
constexpr operator long long(){return 1ll<<60;}
constexpr auto operator-(){
struct{
constexpr operator int(){return -1<<30;}
constexpr operator long long(){return -1ll<<60;}
} res;
return res;
}
} inf;
constexpr int INF=1<<30;
constexpr long long LINF=1ll<<60;
//#define endl "\n"
template<class T> inline T gcd(T a,T b){return b?gcd(b,a%b):a;}
template<class T> inline ll lcm(T a,T b){return (log2(a)-log2(gcd(a,b))+log2(b)<62)?(a/gcd(a,b)*b<LINF?a/gcd(a,b)*b:LINF):LINF;}
const vector<int> dx={1,0,-1,0};
const vector<int> dy={0,1,0,-1};
constexpr double pi=3.141592653589793;
constexpr double PI=3.141592653589793;
#define DEBUG(x) cerr<<__LINE__<<#x<<": "<<x<<endl;
#define ALL(vec) (vec).begin(),(vec).end()
template<typename T> inline void SORT(vector<T>& vec){sort(vec.begin(),vec.end());}
template<typename T> inline void RSORT(vector<T>& vec){sort(vec.begin(),vec.end(),greater<T>());}
template<typename T> inline void ans(T res) {cout<<res<<endl;exit(0);}
template<typename T> inline void yesno(T b) {ans((bool)(b)?"Yes":"No");}
inline void yes(){yesno(1);}
inline void no(){yesno(0);}
ostream &operator<<(ostream &os,const mint &x){
os<<x.val(); return os;
}
template<typename T,typename U> ostream &operator<<(ostream &os,const pair<T,U> &p){
os<<p.first<<" "<<p.second; return os;
}
template<typename T,typename U> istream &operator>>(istream& is,pair<T,U>& p){
is>>p.first>>p.second; return is;
}
template<typename T> ostream &operator<<(ostream &os,const vector<T> &vec){
for(int i=0;i<(int)vec.size();i++){os<<vec[i]<<(i+1!=vec.size()?" ":"");}return os;
}
template<typename T> istream &operator>>(istream& is,vector<T>& vec){
for(int i=0;i<(int)vec.size();i++){is>>vec[i];}return is;
}
template<typename T> void operator+=(vector<T>& vec, int n){
for(T& e:vec) e+=n;
}
template<typename T> void operator*=(vector<T>& vec, int n){
for(T& e:vec) e*=n;
}
template<typename T> vector<T> operator+(const vector<T>& vec0, const vector<T> vec1){
int n=vec0.size();
if(n!=vec1.size()){cerr<<"ERROR: oparator*(const vector<T>&, const vector<T>&): v0.size()!=v1.size()"<<endl;n=min(n,(int)vec1.size());}
vector<T> res;
for(int i=0;i<n;i++) res[i]=vec0[i]+vec1[i];
return res;
}
// 内積
template<typename T> T operator*(const vector<T>& vec0, const vector<T> vec1){
int n=vec0.size();
if(n!=vec1.size()){cerr<<"ERROR: oparator*(const vector<T>&, const vector<T>&): v0.size()!=v1.size()"<<endl;n=min(n,(int)vec1.size());}
T res=0;
for(int i=0;i<n;i++) res+=vec0[i]*vec1[i];
return res;
}
// auto vec=make_v(3,4,5,6) で vvvi vec=vvvi(3,vvi(4,vi(5,6))) と同じような挙動
// cf. https://web.archive.org/web/20200506071437/https://beet-aizu.hatenablog.com/entry/2018/04/08/145516
template<typename T> vector<T> make_v(size_t a,T b){return vector<T>(a,b);}
template<typename... Ts> auto make_v(size_t a,Ts... ts){return vector<decltype(make_v(ts...))>(a,make_v(ts...));}
vi range(int N){vi res(N); for(int i=0;i<N;i++) res[i]=i; return res;}
template<typename T> vector<T> range(vector<T>& v){return v;}
#define FOR_PERM(v,n) for(auto (v)=range(n); (v).size()!=0; next_permutation((v).begin(),(v).end()) || ((v)={},0))
vvi forfor(vi v, int depth, const function<int(int,int)>& fbegin, const function<int(int,int)>& fend, int prev, int sum){
if(depth==0) return {v};
int a=fbegin(prev,sum);
int b=fend(prev,sum);
vvi res;
for(int i=a;i<b;i++){
vi vc=v;
vc.push_back(i);
vvi tmp=forfor(vc,depth-1,fbegin,fend,i,sum+i);
for(vi e:tmp) res.push_back(e);
}
return res;
}
#define FORFOR(v,depth,fbegin,fend) for(vi v:forfor({},depth,[&](int prev,int sum){return fbegin;},[&](int prev,int sum){return fend;},0,0))
//エラトステネスの篩 O(NloglogN)
//N以下の素数全体の集合(vector<T>)を返す
template<typename T> vector<T> primes(T N) {
vector<T> res;
vector<bool> isprime(N+1,true); isprime[0]=isprime[1]=false;
for(T p=2;p<N+1;p++){
if(!isprime[p])continue;
res.push_back(p);
for(T q=p*2;q<=N;q+=p){isprime[q]=false;}
}
return res;
}
// ミラーラビン法
namespace IsPrime{
// A^N mod M
template<class T> T pow_mod(T A, T N, T M) {
T res = 1 % M;
A %= M;
while (N) {
if (N & 1) res = (res * A) % M;
A = (A * A) % M;
N >>= 1;
}
return res;
}
bool MillerRabin(long long N, vector<long long> A) {
long long s = 0, d = N - 1;
while (d % 2 == 0) {
++s;
d >>= 1;
}
for (auto a : A) {
if (N <= a) return true;
long long t, x = pow_mod<__int128_t>(a, d, N);
if (x != 1) {
for (t = 0; t < s; ++t) {
if (x == N - 1) break;
x = __int128_t(x) * x % N;
}
if (t == s) return false;
}
}
return true;
}
bool is_prime(long long N) {
if (N <= 1) return false;
if (N == 2) return true;
if (N % 2 == 0) return false;
if (N < 4759123141LL)
return MillerRabin(N, {2, 7, 61});
else
return MillerRabin(N, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
}
}
using IsPrime::is_prime;
//素因数分解
template<typename T>
vector<pair<T,int>> prime_factories(T N){
vector<pair<T,int>> res;
for(T a=2;a*a<=N;a++){
if(N%a!=0)continue;
int ex=0;
while(N%a==0){
ex++;
N/=a;
}
res.pb({a,ex});
}
if(N!=1)res.pb({N,1});
return res;
}
// 二分探索:
// !f(i-1) && f(i) となるようなiのうちいずれか1つを返すか、
// f(0)==f(N) のとき f(0) なら 0, !f(0) なら N を返す
// 制約 : f(0) ⇒ f(N)
template<typename T>
T binarySearch(T N, function<bool(T)> f){
T min=-1, max=N;
for(;max-min>1;){
T mid=(min+max)/2;
if(f(mid)) max=mid;
else min=mid;
}
return max;
}
// 有理数型
template<typename T> class fraction{
public:
// 分母, 分子
T top, bottom;
fraction(T top_,T bottom_){top=top_;bottom=bottom_;if(bottom==0) ERROR("fraction<int> constractor: 0-devided error");}
friend bool operator<(fraction const& a,fraction const& b){
if(a.bottom>0!=b.bottom>0) return (ll)a.top*(ll)b.bottom<(ll)b.top*(ll)a.bottom;
else return (ll)a.top*(ll)b.bottom>(ll)b.top*(ll)a.bottom;
}
friend bool operator>(fraction const& a,fraction const& b){
if(a.bottom>0!=b.bottom>0) return (ll)a.top*(ll)b.bottom>(ll)b.top*(ll)a.bottom;
else return (ll)a.top*(ll)b.bottom<(ll)b.top*(ll)a.bottom;
}
friend bool operator!=(fraction const& a,fraction const& b){return a>b||a<b;}
friend bool operator>=(fraction const& a,fraction const& b){return a>b||a==b;}
friend bool operator<=(fraction const& a,fraction const& b){return a<b||a==b;}
friend bool operator==(fraction const& a,fraction const& b){return !(a!=b);}
};
// マンハッタン距離
template<class T> T L1_norm(T x0,T y0,T x1,T y1){return abs(x0-x1)+abs(y0-y1);}
// ユークリッド距離
template<class T> double L2_norm(T x0,T y0,T x1,T y1){T dx=x0-x1,dy=y0-y1; return sqrt(dx*dx+dy*dy);}
// kitamasa法 O(k^2*logk)
// 漸化式 a[n+k]=sigma{0<=i<k}(d[i]*a[n+i]) で表される数列aの第n項を求めます
// 引数のaは数列aの第0~k-1項をわたすこと
// remind:漸化式に定数項を含む場合は各項に適切な定数を加えることで定数項を含まない形に変形できます
// cf.https://outline.hatenadiary.jp/entry/2020/07/02/205628
template<typename T,typename U> T kitamasa(const vector<T>& a,const vector<T>& d,U n){
int k = a.size();
if(n < k) return a[n];
function<vector<T>(U)> dfs=[&a,&d,k,&dfs](U n){
if(n == k) return d;
vector<T> res(k);
if(n & 1 || n < k * 2){
vector<T> x = dfs(n - 1);
for(int i = 0; i < k; ++i) res[i] = d[i] * x[k - 1];
for(int i = 0; i + 1 < k; ++i) res[i + 1] += x[i];
}
else{
vector<vector<T>> x(k, vector<T>(k));
x[0] = dfs(n >> 1);
for(int i = 0; i + 1 < k; ++i){
for(int j = 0; j < k; ++j) x[i + 1][j] = d[j] * x[i][k - 1];
for(int j = 0; j + 1 < k; ++j) x[i + 1][j + 1] += x[i][j];
}
for(int i = 0; i < k; ++i){
for(int j = 0; j < k; ++j){
res[j] += x[0][i] * x[i][j];
}
}
}
return res;
};
// a_n を求める
vector<T> x = dfs(n);
T res = 0;
for(int i = 0; i < k; ++i) res += x[i] * a[i];
return res;
}
// cf. https://ei1333.github.io/library/
namespace luzhiled{
#line 2 "math/fps/formal-power-series-friendly-ntt.hpp"
#line 1 "math/fft/number-theoretic-transform-friendly-mod-int.hpp"
/**
* @brief Number Theoretic Transform Friendly ModInt
*/
template <typename Mint>
struct NumberTheoreticTransformFriendlyModInt {
static vector<Mint> roots, iroots, rate3, irate3;
static int max_base;
NumberTheoreticTransformFriendlyModInt() = default;
static void init() {
if (roots.empty()) {
const unsigned mod = Mint::mod();
assert(mod >= 3 && mod % 2 == 1);
auto tmp = mod - 1;
max_base = 0;
while (tmp % 2 == 0) tmp >>= 1, max_base++;
Mint root = 2;
while (root.pow((mod - 1) >> 1) == 1) {
root += 1;
}
assert(root.pow(mod - 1) == 1);
roots.resize(max_base + 1);
iroots.resize(max_base + 1);
rate3.resize(max_base + 1);
irate3.resize(max_base + 1);
roots[max_base] = root.pow((mod - 1) >> max_base);
iroots[max_base] = Mint(1) / roots[max_base];
for (int i = max_base - 1; i >= 0; i--) {
roots[i] = roots[i + 1] * roots[i + 1];
iroots[i] = iroots[i + 1] * iroots[i + 1];
}
{
Mint prod = 1, iprod = 1;
for (int i = 0; i <= max_base - 3; i++) {
rate3[i] = roots[i + 3] * prod;
irate3[i] = iroots[i + 3] * iprod;
prod *= iroots[i + 3];
iprod *= roots[i + 3];
}
}
}
}
static void ntt(vector<Mint> &a) {
init();
const int n = (int)a.size();
assert((n & (n - 1)) == 0);
int h = __builtin_ctz(n);
assert(h <= max_base);
int len = 0;
Mint imag = roots[2];
if (h & 1) {
int p = 1 << (h - 1);
Mint rot = 1;
for (int i = 0; i < p; i++) {
auto r = a[i + p];
a[i + p] = a[i] - r;
a[i] += r;
}
len++;
}
for (; len + 1 < h; len += 2) {
int p = 1 << (h - len - 2);
{ // s = 0
for (int i = 0; i < p; i++) {
auto a0 = a[i];
auto a1 = a[i + p];
auto a2 = a[i + 2 * p];
auto a3 = a[i + 3 * p];
auto a1na3imag = (a1 - a3) * imag;
auto a0a2 = a0 + a2;
auto a1a3 = a1 + a3;
auto a0na2 = a0 - a2;
a[i] = a0a2 + a1a3;
a[i + 1 * p] = a0a2 - a1a3;
a[i + 2 * p] = a0na2 + a1na3imag;
a[i + 3 * p] = a0na2 - a1na3imag;
}
}
Mint rot = rate3[0];
for (int s = 1; s < (1 << len); s++) {
int offset = s << (h - len);
Mint rot2 = rot * rot;
Mint rot3 = rot2 * rot;
for (int i = 0; i < p; i++) {
auto a0 = a[i + offset];
auto a1 = a[i + offset + p] * rot;
auto a2 = a[i + offset + 2 * p] * rot2;
auto a3 = a[i + offset + 3 * p] * rot3;
auto a1na3imag = (a1 - a3) * imag;
auto a0a2 = a0 + a2;
auto a1a3 = a1 + a3;
auto a0na2 = a0 - a2;
a[i + offset] = a0a2 + a1a3;
a[i + offset + 1 * p] = a0a2 - a1a3;
a[i + offset + 2 * p] = a0na2 + a1na3imag;
a[i + offset + 3 * p] = a0na2 - a1na3imag;
}
rot *= rate3[__builtin_ctz(~s)];
}
}
}
static void intt(vector<Mint> &a, bool f = true) {
init();
const int n = (int)a.size();
assert((n & (n - 1)) == 0);
int h = __builtin_ctz(n);
assert(h <= max_base);
int len = h;
Mint iimag = iroots[2];
for (; len > 1; len -= 2) {
int p = 1 << (h - len);
{ // s = 0
for (int i = 0; i < p; i++) {
auto a0 = a[i];
auto a1 = a[i + 1 * p];
auto a2 = a[i + 2 * p];
auto a3 = a[i + 3 * p];
auto a2na3iimag = (a2 - a3) * iimag;
auto a0na1 = a0 - a1;
auto a0a1 = a0 + a1;
auto a2a3 = a2 + a3;
a[i] = a0a1 + a2a3;
a[i + 1 * p] = (a0na1 + a2na3iimag);
a[i + 2 * p] = (a0a1 - a2a3);
a[i + 3 * p] = (a0na1 - a2na3iimag);
}
}
Mint irot = irate3[0];
for (int s = 1; s < (1 << (len - 2)); s++) {
int offset = s << (h - len + 2);
Mint irot2 = irot * irot;
Mint irot3 = irot2 * irot;
for (int i = 0; i < p; i++) {
auto a0 = a[i + offset];
auto a1 = a[i + offset + 1 * p];
auto a2 = a[i + offset + 2 * p];
auto a3 = a[i + offset + 3 * p];
auto a2na3iimag = (a2 - a3) * iimag;
auto a0na1 = a0 - a1;
auto a0a1 = a0 + a1;
auto a2a3 = a2 + a3;
a[i + offset] = a0a1 + a2a3;
a[i + offset + 1 * p] = (a0na1 + a2na3iimag) * irot;
a[i + offset + 2 * p] = (a0a1 - a2a3) * irot2;
a[i + offset + 3 * p] = (a0na1 - a2na3iimag) * irot3;
}
irot *= irate3[__builtin_ctz(~s)];
}
}
if (len >= 1) {
int p = 1 << (h - 1);
for (int i = 0; i < p; i++) {
auto ajp = a[i] - a[i + p];
a[i] += a[i + p];
a[i + p] = ajp;
}
}
if (f) {
Mint inv_sz = Mint(1) / n;
for (int i = 0; i < n; i++) a[i] *= inv_sz;
}
}
static vector<Mint> multiply(vector<Mint> a, vector<Mint> b) {
int need = a.size() + b.size() - 1;
int nbase = 1;
while ((1 << nbase) < need) nbase++;
int sz = 1 << nbase;
a.resize(sz, 0);
b.resize(sz, 0);
ntt(a);
ntt(b);
Mint inv_sz = Mint(1) / sz;
for (int i = 0; i < sz; i++) a[i] *= b[i] * inv_sz;
intt(a, false);
a.resize(need);
return a;
}
};
template <typename Mint>
vector<Mint> NumberTheoreticTransformFriendlyModInt<Mint>::roots =
vector<Mint>();
template <typename Mint>
vector<Mint> NumberTheoreticTransformFriendlyModInt<Mint>::iroots =
vector<Mint>();
template <typename Mint>
vector<Mint> NumberTheoreticTransformFriendlyModInt<Mint>::rate3 =
vector<Mint>();
template <typename Mint>
vector<Mint> NumberTheoreticTransformFriendlyModInt<Mint>::irate3 =
vector<Mint>();
template <typename Mint>
int NumberTheoreticTransformFriendlyModInt<Mint>::max_base = 0;
#line 4 "math/fps/formal-power-series-friendly-ntt.hpp"
template <typename T>
struct FormalPowerSeriesFriendlyNTT : vector<T> {
using vector<T>::vector;
using P = FormalPowerSeriesFriendlyNTT;
using NTT = NumberTheoreticTransformFriendlyModInt<T>;
P pre(int deg) const {
return P(begin(*this), begin(*this) + min((int)this->size(), deg));
}
P rev(int deg = -1) const {
P ret(*this);
if (deg != -1) ret.resize(deg, T(0));
reverse(begin(ret), end(ret));
return ret;
}
void shrink() {
while (this->size() && this->back() == T(0)) this->pop_back();
}
P operator+(const P &r) const { return P(*this) += r; }
P operator+(const T &v) const { return P(*this) += v; }
P operator-(const P &r) const { return P(*this) -= r; }
P operator-(const T &v) const { return P(*this) -= v; }
P operator*(const P &r) const { return P(*this) *= r; }
P operator*(const T &v) const { return P(*this) *= v; }
P operator/(const P &r) const { return P(*this) /= r; }
P operator%(const P &r) const { return P(*this) %= r; }
P &operator+=(const P &r) {
if (r.size() > this->size()) this->resize(r.size());
for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i];
return *this;
}
P &operator-=(const P &r) {
if (r.size() > this->size()) this->resize(r.size());
for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i];
return *this;
}
// https://judge.yosupo.jp/problem/convolution_mod
P &operator*=(const P &r) {
if (this->empty() || r.empty()) {
this->clear();
return *this;
}
auto ret = NTT::multiply(*this, r);
return *this = {begin(ret), end(ret)};
}
P &operator/=(const P &r) {
if (this->size() < r.size()) {
this->clear();
return *this;
}
int n = this->size() - r.size() + 1;
return *this = (rev().pre(n) * r.rev().inv(n)).pre(n).rev(n);
}
P &operator%=(const P &r) {
*this -= *this / r * r;
shrink();
return *this;
}
// https://judge.yosupo.jp/problem/division_of_polynomials
pair<P, P> div_mod(const P &r) {
P q = *this / r;
P x = *this - q * r;
x.shrink();
return make_pair(q, x);
}
P operator-() const {
P ret(this->size());
for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i];
return ret;
}
P &operator+=(const T &r) {
if (this->empty()) this->resize(1);
(*this)[0] += r;
return *this;
}
P &operator-=(const T &r) {
if (this->empty()) this->resize(1);
(*this)[0] -= r;
return *this;
}
P &operator*=(const T &v) {
for (int i = 0; i < (int)this->size(); i++) (*this)[i] *= v;
return *this;
}
P dot(P r) const {
P ret(min(this->size(), r.size()));
for (int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i];
return ret;
}
P operator>>(int sz) const {
if ((int)this->size() <= sz) return {};
P ret(*this);
ret.erase(ret.begin(), ret.begin() + sz);
return ret;
}
P operator<<(int sz) const {
P ret(*this);
ret.insert(ret.begin(), sz, T(0));
return ret;
}
T operator()(T x) const {
T r = 0, w = 1;
for (auto &v : *this) {
r += w * v;
w *= x;
}
return r;
}
P diff() const {
const int n = (int)this->size();
P ret(max(0, n - 1));
for (int i = 1; i < n; i++) ret[i - 1] = (*this)[i] * T(i);
return ret;
}
P integral() const {
const int n = (int)this->size();
P ret(n + 1);
ret[0] = T(0);
for (int i = 0; i < n; i++) ret[i + 1] = (*this)[i] / T(i + 1);
return ret;
}
// https://judge.yosupo.jp/problem/inv_of_formal_power_series
// F(0) must not be 0
P inv(int deg = -1) const {
assert(((*this)[0]) != T(0));
const int n = (int)this->size();
if (deg == -1) deg = n;
P res(deg);
res[0] = {T(1) / (*this)[0]};
for (int d = 1; d < deg; d <<= 1) {
P f(2 * d), g(2 * d);
for (int j = 0; j < min(n, 2 * d); j++) f[j] = (*this)[j];
for (int j = 0; j < d; j++) g[j] = res[j];
NTT::ntt(f);
NTT::ntt(g);
f = f.dot(g);
NTT::intt(f);
for (int j = 0; j < d; j++) f[j] = 0;
NTT::ntt(f);
for (int j = 0; j < 2 * d; j++) f[j] *= g[j];
NTT::intt(f);
for (int j = d; j < min(2 * d, deg); j++) res[j] = -f[j];
}
return res;
}
// https://judge.yosupo.jp/problem/log_of_formal_power_series
// F(0) must be 1
P log(int deg = -1) const {
assert((*this)[0] == T(1));
const int n = (int)this->size();
if (deg == -1) deg = n;
return (this->diff() * this->inv(deg)).pre(deg - 1).integral();
}
// https://judge.yosupo.jp/problem/sqrt_of_formal_power_series
P sqrt(
int deg = -1,
const function<T(T)> &get_sqrt = [](T) { return T(1); }) const {
const int n = (int)this->size();
if (deg == -1) deg = n;
if ((*this)[0] == T(0)) {
for (int i = 1; i < n; i++) {
if ((*this)[i] != T(0)) {
if (i & 1) return {};
if (deg - i / 2 <= 0) break;
auto ret = (*this >> i).sqrt(deg - i / 2, get_sqrt);
if (ret.empty()) return {};
ret = ret << (i / 2);
if ((int)ret.size() < deg) ret.resize(deg, T(0));
return ret;
}
}
return P(deg, 0);
}
auto sqr = T(get_sqrt((*this)[0]));
if (sqr * sqr != (*this)[0]) return {};
P ret{sqr};
T inv2 = T(1) / T(2);
for (int i = 1; i < deg; i <<= 1) {
ret = (ret + pre(i << 1) * ret.inv(i << 1)) * inv2;
}
return ret.pre(deg);
}
P sqrt(const function<T(T)> &get_sqrt, int deg = -1) const {
return sqrt(deg, get_sqrt);
}
// https://judge.yosupo.jp/problem/exp_of_formal_power_series
// F(0) must be 0
P exp(int deg = -1) const {
if (deg == -1) deg = this->size();
assert((*this)[0] == T(0));
P inv;
inv.reserve(deg + 1);
inv.push_back(T(0));
inv.push_back(T(1));
auto inplace_integral = [&](P &F) -> void {
const int n = (int)F.size();
auto mod = T::mod();
while ((int)inv.size() <= n) {
int i = inv.size();
inv.push_back((-inv[mod % i]) * (mod / i));
}
F.insert(begin(F), T(0));
for (int i = 1; i <= n; i++) F[i] *= inv[i];
};
auto inplace_diff = [](P &F) -> void {
if (F.empty()) return;
F.erase(begin(F));
T coeff = 1, one = 1;
for (int i = 0; i < (int)F.size(); i++) {
F[i] *= coeff;
coeff += one;
}
};
P b{1, 1 < (int)this->size() ? (*this)[1] : 0}, c{1}, z1, z2{1, 1};
for (int m = 2; m < deg; m *= 2) {
auto y = b;
y.resize(2 * m);
NTT::ntt(y);
z1 = z2;
P z(m);
for (int i = 0; i < m; ++i) z[i] = y[i] * z1[i];
NTT::intt(z);
fill(begin(z), begin(z) + m / 2, T(0));
NTT::ntt(z);
for (int i = 0; i < m; ++i) z[i] *= -z1[i];
NTT::intt(z);
c.insert(end(c), begin(z) + m / 2, end(z));
z2 = c;
z2.resize(2 * m);
NTT::ntt(z2);
P x(begin(*this), begin(*this) + min<int>(this->size(), m));
inplace_diff(x);
x.push_back(T(0));
NTT::ntt(x);
for (int i = 0; i < m; ++i) x[i] *= y[i];
NTT::intt(x);
x -= b.diff();
x.resize(2 * m);
for (int i = 0; i < m - 1; ++i) x[m + i] = x[i], x[i] = T(0);
NTT::ntt(x);
for (int i = 0; i < 2 * m; ++i) x[i] *= z2[i];
NTT::intt(x);
x.pop_back();
inplace_integral(x);
for (int i = m; i < min<int>(this->size(), 2 * m); ++i)
x[i] += (*this)[i];
fill(begin(x), begin(x) + m, T(0));
NTT::ntt(x);
for (int i = 0; i < 2 * m; ++i) x[i] *= y[i];
NTT::intt(x);
b.insert(end(b), begin(x) + m, end(x));
}
return P{begin(b), begin(b) + deg};
}
// https://judge.yosupo.jp/problem/pow_of_formal_power_series
P pow(int64_t k, int deg = -1) const {
const int n = (int)this->size();
if (deg == -1) deg = n;
if (k == 0) {
P ret(deg, T(0));
ret[0] = T(1);
return ret;
}
for (int i = 0; i < n; i++) {
if (i * k > deg) return P(deg, T(0));
if ((*this)[i] != T(0)) {
T rev = T(1) / (*this)[i];
P ret = (((*this * rev) >> i).log() * k).exp() * ((*this)[i].pow(k));
ret = (ret << (i * k)).pre(deg);
if ((int)ret.size() < deg) ret.resize(deg, T(0));
return ret;
}
}
return *this;
}
P mod_pow(int64_t k, P g) const {
P modinv = g.rev().inv();
auto get_div = [&](P base) {
if (base.size() < g.size()) {
base.clear();
return base;
}
int n = base.size() - g.size() + 1;
return (base.rev().pre(n) * modinv.pre(n)).pre(n).rev(n);
};
P x(*this), ret{1};
while (k > 0) {
if (k & 1) {
ret *= x;
ret -= get_div(ret) * g;
ret.shrink();
}
x *= x;
x -= get_div(x) * g;
x.shrink();
k >>= 1;
}
return ret;
}
// https://judge.yosupo.jp/problem/polynomial_taylor_shift
P taylor_shift(T c) const {
int n = (int)this->size();
vector<T> fact(n), rfact(n);
fact[0] = rfact[0] = T(1);
for (int i = 1; i < n; i++) fact[i] = fact[i - 1] * T(i);
rfact[n - 1] = T(1) / fact[n - 1];
for (int i = n - 1; i > 1; i--) rfact[i - 1] = rfact[i] * T(i);
P p(*this);
for (int i = 0; i < n; i++) p[i] *= fact[i];
p = p.rev();
P bs(n, T(1));
for (int i = 1; i < n; i++) bs[i] = bs[i - 1] * c * rfact[i] * fact[i - 1];
p = (p * bs).pre(n);
p = p.rev();
for (int i = 0; i < n; i++) p[i] *= rfact[i];
return p;
}
};
template <typename Mint>
using FPSNTT = FormalPowerSeriesFriendlyNTT<Mint>;
#line 2 "math/fps/formal-power-series.hpp"
#line 1 "math/fft/fast-fourier-transform.hpp"
namespace FastFourierTransform {
using real = double;
struct C {
real x, y;
C() : x(0), y(0) {}
C(real x, real y) : x(x), y(y) {}
inline C operator+(const C &c) const { return C(x + c.x, y + c.y); }
inline C operator-(const C &c) const { return C(x - c.x, y - c.y); }
inline C operator*(const C &c) const {
return C(x * c.x - y * c.y, x * c.y + y * c.x);
}
inline C conj() const { return C(x, -y); }
};
const real PI = acosl(-1);
int base = 1;
vector<C> rts = {{0, 0}, {1, 0}};
vector<int> rev = {0, 1};
void ensure_base(int nbase) {
if (nbase <= base) return;
rev.resize(1 << nbase);
rts.resize(1 << nbase);
for (int i = 0; i < (1 << nbase); i++) {
rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (nbase - 1));
}
while (base < nbase) {
real angle = PI * 2.0 / (1 << (base + 1));
for (int i = 1 << (base - 1); i < (1 << base); i++) {
rts[i << 1] = rts[i];
real angle_i = angle * (2 * i + 1 - (1 << base));
rts[(i << 1) + 1] = C(cos(angle_i), sin(angle_i));
}
++base;
}
}
void fft(vector<C> &a, int n) {
assert((n & (n - 1)) == 0);
int zeros = __builtin_ctz(n);
ensure_base(zeros);
int shift = base - zeros;
for (int i = 0; i < n; i++) {
if (i < (rev[i] >> shift)) {
swap(a[i], a[rev[i] >> shift]);
}
}
for (int k = 1; k < n; k <<= 1) {
for (int i = 0; i < n; i += 2 * k) {
for (int j = 0; j < k; j++) {
C z = a[i + j + k] * rts[j + k];
a[i + j + k] = a[i + j] - z;
a[i + j] = a[i + j] + z;
}
}
}
}
vector<int64_t> multiply(const vector<int> &a, const vector<int> &b) {
int need = (int)a.size() + (int)b.size() - 1;
int nbase = 1;
while ((1 << nbase) < need) nbase++;
ensure_base(nbase);
int sz = 1 << nbase;
vector<C> fa(sz);
for (int i = 0; i < sz; i++) {
int x = (i < (int)a.size() ? a[i] : 0);
int y = (i < (int)b.size() ? b[i] : 0);
fa[i] = C(x, y);
}
fft(fa, sz);
C r(0, -0.25 / (sz >> 1)), s(0, 1), t(0.5, 0);
for (int i = 0; i <= (sz >> 1); i++) {
int j = (sz - i) & (sz - 1);
C z = (fa[j] * fa[j] - (fa[i] * fa[i]).conj()) * r;
fa[j] = (fa[i] * fa[i] - (fa[j] * fa[j]).conj()) * r;
fa[i] = z;
}
for (int i = 0; i < (sz >> 1); i++) {
C A0 = (fa[i] + fa[i + (sz >> 1)]) * t;
C A1 = (fa[i] - fa[i + (sz >> 1)]) * t * rts[(sz >> 1) + i];
fa[i] = A0 + A1 * s;
}
fft(fa, sz >> 1);
vector<int64_t> ret(need);
for (int i = 0; i < need; i++) {
ret[i] = llround(i & 1 ? fa[i >> 1].y : fa[i >> 1].x);
}
return ret;
}
}; // namespace FastFourierTransform
#line 2 "math/fft/arbitrary-mod-convolution.hpp"
/*
* @brief Arbitrary Mod Convolution(任意mod畳み込み)
*/
template <typename T>
struct ArbitraryModConvolution {
using real = FastFourierTransform::real;
using C = FastFourierTransform::C;
ArbitraryModConvolution() = default;
static vector<T> multiply(const vector<T> &a, const vector<T> &b,
int need = -1) {
if (need == -1) need = a.size() + b.size() - 1;
int nbase = 0;
while ((1 << nbase) < need) nbase++;
FastFourierTransform::ensure_base(nbase);
int sz = 1 << nbase;
vector<C> fa(sz);
for (int i = 0; i < a.size(); i++) {
fa[i] = C(a[i].val() & ((1 << 15) - 1), a[i].val() >> 15);
}
fft(fa, sz);
vector<C> fb(sz);
if (a == b) {
fb = fa;
} else {
for (int i = 0; i < b.size(); i++) {
fb[i] = C(b[i].val() & ((1 << 15) - 1), b[i].val() >> 15);
}
fft(fb, sz);
}
real ratio = 0.25 / sz;
C r2(0, -1), r3(ratio, 0), r4(0, -ratio), r5(0, 1);
for (int i = 0; i <= (sz >> 1); i++) {
int j = (sz - i) & (sz - 1);
C a1 = (fa[i] + fa[j].conj());
C a2 = (fa[i] - fa[j].conj()) * r2;
C b1 = (fb[i] + fb[j].conj()) * r3;
C b2 = (fb[i] - fb[j].conj()) * r4;
if (i != j) {
C c1 = (fa[j] + fa[i].conj());
C c2 = (fa[j] - fa[i].conj()) * r2;
C d1 = (fb[j] + fb[i].conj()) * r3;
C d2 = (fb[j] - fb[i].conj()) * r4;
fa[i] = c1 * d1 + c2 * d2 * r5;
fb[i] = c1 * d2 + c2 * d1;
}
fa[j] = a1 * b1 + a2 * b2 * r5;
fb[j] = a1 * b2 + a2 * b1;
}
fft(fa, sz);
fft(fb, sz);
vector<T> ret(need);
for (int i = 0; i < need; i++) {
int64_t aa = llround(fa[i].x);
int64_t bb = llround(fb[i].x);
int64_t cc = llround(fa[i].y);
aa = T(aa).val(), bb = T(bb).val(), cc = T(cc).val();
ret[i] = aa + (bb << 15) + (cc << 30);
}
return ret;
}
};
#line 4 "math/fps/formal-power-series.hpp"
/**
* @brief Formal Power Series(形式的冪級数)
*/
template <typename T>
struct FormalPowerSeries : vector<T> {
using vector<T>::vector;
using P = FormalPowerSeries;
using Conv = ArbitraryModConvolution<T>;
P pre(int deg) const {
return P(begin(*this), begin(*this) + min((int)this->size(), deg));
}
P rev(int deg = -1) const {
P ret(*this);
if (deg != -1) ret.resize(deg, T(0));
reverse(begin(ret), end(ret));
return ret;
}
void shrink() {
while (this->size() && this->back() == T(0)) this->pop_back();
}
P operator+(const P &r) const { return P(*this) += r; }
P operator+(const T &v) const { return P(*this) += v; }
P operator-(const P &r) const { return P(*this) -= r; }
P operator-(const T &v) const { return P(*this) -= v; }
P operator*(const P &r) const { return P(*this) *= r; }
P operator*(const T &v) const { return P(*this) *= v; }
P operator/(const P &r) const { return P(*this) /= r; }
P operator%(const P &r) const { return P(*this) %= r; }
P &operator+=(const P &r) {
if (r.size() > this->size()) this->resize(r.size());
for (int i = 0; i < r.size(); i++) (*this)[i] += r[i];
return *this;
}
P &operator-=(const P &r) {
if (r.size() > this->size()) this->resize(r.size());
for (int i = 0; i < r.size(); i++) (*this)[i] -= r[i];
return *this;
}
// https://judge.yosupo.jp/problem/convolution_mod
P &operator*=(const P &r) {
if (this->empty() || r.empty()) {
this->clear();
return *this;
}
auto ret = Conv::multiply(*this, r);
return *this = {begin(ret), end(ret)};
}
P &operator/=(const P &r) {
if (this->size() < r.size()) {
this->clear();
return *this;
}
int n = this->size() - r.size() + 1;
return *this = (rev().pre(n) * r.rev().inv(n)).pre(n).rev(n);
}
P &operator%=(const P &r) { return *this -= *this / r * r; }
// https://judge.yosupo.jp/problem/division_of_polynomials
pair<P, P> div_mod(const P &r) {
P q = *this / r;
return make_pair(q, *this - q * r);
}
P operator-() const {
P ret(this->size());
for (int i = 0; i < this->size(); i++) ret[i] = -(*this)[i];
return ret;
}
P &operator+=(const T &r) {
if (this->empty()) this->resize(1);
(*this)[0] += r;
return *this;
}
P &operator-=(const T &r) {
if (this->empty()) this->resize(1);
(*this)[0] -= r;
return *this;
}
P &operator*=(const T &v) {
for (int i = 0; i < this->size(); i++) (*this)[i] *= v;
return *this;
}
P dot(P r) const {
P ret(min(this->size(), r.size()));
for (int i = 0; i < ret.size(); i++) ret[i] = (*this)[i] * r[i];
return ret;
}
P operator>>(int sz) const {
if (this->size() <= sz) return {};
P ret(*this);
ret.erase(ret.begin(), ret.begin() + sz);
return ret;
}
P operator<<(int sz) const {
P ret(*this);
ret.insert(ret.begin(), sz, T(0));
return ret;
}
T operator()(T x) const {
T r = 0, w = 1;
for (auto &v : *this) {
r += w * v;
w *= x;
}
return r;
}
P diff() const {
const int n = (int)this->size();
P ret(max(0, n - 1));
for (int i = 1; i < n; i++) ret[i - 1] = (*this)[i] * T(i);
return ret;
}
P integral() const {
const int n = (int)this->size();
P ret(n + 1);
ret[0] = T(0);
for (int i = 0; i < n; i++) ret[i + 1] = (*this)[i] / T(i + 1);
return ret;
}
// https://judge.yosupo.jp/problem/inv_of_formal_power_series
// F(0) must not be 0
P inv(int deg = -1) const {
assert(((*this)[0]) != T(0));
const int n = (int)this->size();
if (deg == -1) deg = n;
P ret({T(1) / (*this)[0]});
for (int i = 1; i < deg; i <<= 1) {
ret = (ret + ret - ret * ret * pre(i << 1)).pre(i << 1);
}
return ret.pre(deg);
}
// https://judge.yosupo.jp/problem/log_of_formal_power_series
// F(0) must be 1
P log(int deg = -1) const {
assert((*this)[0] == T(1));
const int n = (int)this->size();
if (deg == -1) deg = n;
return (this->diff() * this->inv(deg)).pre(deg - 1).integral();
}
// https://judge.yosupo.jp/problem/sqrt_of_formal_power_series
P sqrt(
int deg = -1,
const function<T(T)> &get_sqrt = [](T) { return T(1); }) const {
const int n = (int)this->size();
if (deg == -1) deg = n;
if ((*this)[0] == T(0)) {
for (int i = 1; i < n; i++) {
if ((*this)[i] != T(0)) {
if (i & 1) return {};
if (deg - i / 2 <= 0) break;
auto ret = (*this >> i).sqrt(deg - i / 2, get_sqrt);
if (ret.empty()) return {};
ret = ret << (i / 2);
if (ret.size() < deg) ret.resize(deg, T(0));
return ret;
}
}
return P(deg, 0);
}
auto sqr = T(get_sqrt((*this)[0]));
if (sqr * sqr != (*this)[0]) return {};
P ret{sqr};
T inv2 = T(1) / T(2);
for (int i = 1; i < deg; i <<= 1) {
ret = (ret + pre(i << 1) * ret.inv(i << 1)) * inv2;
}
return ret.pre(deg);
}
P sqrt(const function<T(T)> &get_sqrt, int deg = -1) const {
return sqrt(deg, get_sqrt);
}
// https://judge.yosupo.jp/problem/exp_of_formal_power_series
// F(0) must be 0
P exp(int deg = -1) const {
if (deg == -1) deg = this->size();
assert((*this)[0] == T(0));
const int n = (int)this->size();
if (deg == -1) deg = n;
P ret({T(1)});
for (int i = 1; i < deg; i <<= 1) {
ret = (ret * (pre(i << 1) + T(1) - ret.log(i << 1))).pre(i << 1);
}
return ret.pre(deg);
}
// https://judge.yosupo.jp/problem/pow_of_formal_power_series
P pow(int64_t k, int deg = -1) const {
const int n = (int)this->size();
if (deg == -1) deg = n;
if (k == 0) {
P ret(deg, T(0));
ret[0] = T(1);
return ret;
}
for (int i = 0; i < n; i++) {
if (i * k > deg) return P(deg, T(0));
if ((*this)[i] != T(0)) {
T rev = T(1) / (*this)[i];
P ret = (((*this * rev) >> i).log() * k).exp() * ((*this)[i].pow(k));
ret = (ret << (i * k)).pre(deg);
if (ret.size() < deg) ret.resize(deg, T(0));
return ret;
}
}
return *this;
}
// https://yukicoder.me/problems/no/215
P mod_pow(int64_t k, P g) const {
P modinv = g.rev().inv();
auto get_div = [&](P base) {
if (base.size() < g.size()) {
base.clear();
return base;
}
int n = base.size() - g.size() + 1;
return (base.rev().pre(n) * modinv.pre(n)).pre(n).rev(n);
};
P x(*this), ret{1};
while (k > 0) {
if (k & 1) {
ret *= x;
ret -= get_div(ret) * g;
ret.shrink();
}
x *= x;
x -= get_div(x) * g;
x.shrink();
k >>= 1;
}
return ret;
}
// https://judge.yosupo.jp/problem/polynomial_taylor_shift
P taylor_shift(T c) const {
int n = (int)this->size();
vector<T> fact(n), rfact(n);
fact[0] = rfact[0] = T(1);
for (int i = 1; i < n; i++) fact[i] = fact[i - 1] * T(i);
rfact[n - 1] = T(1) / fact[n - 1];
for (int i = n - 1; i > 1; i--) rfact[i - 1] = rfact[i] * T(i);
P p(*this);
for (int i = 0; i < n; i++) p[i] *= fact[i];
p = p.rev();
P bs(n, T(1));
for (int i = 1; i < n; i++) bs[i] = bs[i - 1] * c * rfact[i] * fact[i - 1];
p = (p * bs).pre(n);
p = p.rev();
for (int i = 0; i < n; i++) p[i] *= rfact[i];
return p;
}
};
template <typename Mint>
using FPSFFT = FormalPowerSeries<Mint>;
};
using luzhiled::FPSNTT;
using luzhiled::FPSFFT;
//template ends here
//main() is at the top of this file
#endif