結果
| 問題 | No.1559 Next Rational |
| コンテスト | |
| ユーザー |
 kakao745
|
| 提出日時 | 2025-12-24 14:16:01 |
| 言語 | C++23 (gcc 13.3.0 + boost 1.89.0) |
| 結果 |
AC
|
| 実行時間 | 2 ms / 2,000 ms |
| コード長 | 34,639 bytes |
| 記録 | |
| コンパイル時間 | 8,758 ms |
| コンパイル使用メモリ | 364,804 KB |
| 実行使用メモリ | 7,848 KB |
| 最終ジャッジ日時 | 2025-12-24 14:16:12 |
| 合計ジャッジ時間 | 10,400 ms |
|
ジャッジサーバーID （参考情報） |
judge4 / judge1 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 3 |
| other | AC * 15 |
ソースコード
#ifdef LOCAL
#include "my_header.h"
#else
#define PRAGMA_OPTIMIZE(s) _Pragma(#s)
//#pragma GCC target ("avx,avx2")//四則演算
PRAGMA_OPTIMIZE(GCC optimize("Ofast"))
PRAGMA_OPTIMIZE(GCC optimize("unroll-loops"))//ループ
//#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,tune=native")//浮動小数点
#include <bits/stdc++.h>
using namespace std;
#if __has_include(<atcoder/all>)
#include <atcoder/all>
using namespace atcoder;
using minta=modint998244353;
using mintb=modint1000000007;
std::ostream&operator<<(std::ostream& os,minta v){os << v.val();return os;}
std::istream&operator>>(std::istream&is,minta &v){long long t;is >> t;v=t;return is;}
#endif
/*//多倍長整数
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/cpp_int.hpp>
namespace mp = boost::multiprecision;
// 任意長整数型
using Bint = mp::cpp_int;
//仮数部が10進数で1024桁の浮動小数点数型(TLEしたら小さくする)
using Real = mp::number<mp::cpp_dec_float<1024>>;
*/
#define rep(i,n) for(long long i=0;i<(long long)n;i++)
#define reps(i,n) for(long long i=1;i<=(long long)n;i++)
#define repi(i,n) for(int i=0;i<(int)n;i++)
#define loop(i,l,r) for(long long i=(long long)l;i<=(long long)r;i++)
#define loopi(i,l,r) for(int i=(int)l;i<=(int)r;i++)
#define drep(i,n) for(long long i=(long long)n-1;i>=0;i--)
#define drepi(i,n) for(int i=(int)n-1;i>=0;i--)
#define dreps(i,n) for(int i=(int)n;i>=1;i--)
#define dloop(i,l,r) for(long long i=(long long)l;i>=(long long)r;i--)
#define dloopi(i,l,r) for(int i=(int)l;i>=(int)r;i--)
#define all(v) v.begin(), v.end()
#define rall(v) v.rbegin(), v.rend()
#define yn(x) cout << (x? "Yes":"No") << endl;
#define cou(x) cout << x << endl;
#define emp emplace_back
#define mp make_pair
const long long moda=998244353LL;
const long long modb=1000000007LL;
const int kaz=1000000005;
const long long yab=2500000000000000000LL;
const long long aho =-yab;
const long double eps=1.0e-14L;
const long double pi=acosl(-1.0L);
using ll=long long;
using st=string;
using P=pair<ll,ll>;
using tup=tuple<ll,ll,ll>;
using vi=vector<ll>;
using vin=vector<int>;
using vc=vector<char>;
using vb=vector<bool>;
using vd=vector<double>;
using vs=vector<string>;
using vp=vector<P>;
using sp=set<P>;
using si=set<ll>;
using vvi=vector<vector<ll>>;
using vvin=vector<vin>;
using vvc=vector<vc>;
using vvb=vector<vb>;
using vvvi=vector<vvi>;
using vvvin=vector<vvin>;
const int dx[4]={0,1,0,-1};
const int dy[4]={1,0,-1,0};
const vector<int> ex = {-1, -1, -1, 0, 0, 1, 1, 1};
const vector<int> ey = {-1, 0, 1, -1, 1, -1, 0, 1};
template<typename T1,typename T2>ostream&operator<<(ostream&os,pair<T1,T2> v){os << '(' << v.first << ',' << v.second << ')';return os;}
template<typename T1,typename T2>istream&operator>>(istream&is,pair<T1,T2> &v){is >> v.first >> v.second;return is;}
template<typename T>istream&operator>>(istream&is,vector<T>&v){for(T&in:v)is>>in;return is;}
template<typename T>ostream&operator<<(ostream&os,vector<T>v){rep(i,v.size())os<<v[i]<<(i+1!=v.size()?" ":"\n");return os;}
template<typename T1,typename T2>
void co(bool x,T1 y,T2 z){
if(x)cout << y << endl;
else cout << z << endl;
}
long long isqrt(long long n){
long long ok=0,ng=1000000000;//1e9
while(ng-ok>1){
long long mid=(ng+ok)/2;
if(mid*mid<=n)ok=mid;
else ng=mid;
}
return ok;
}
template<typename T>
bool chmax(T &a, T b){
if(a<b){
a=b;
return true;
}
return false;
}
template <typename T, typename U>
T ceil(T x, U y) {
return (x > 0 ? (x + y - 1) / y : x / y);
}
template <typename T, typename U>
T floor(T x, U y) {
return (x > 0 ? x / y : (x - y + 1) / y);
}
template<typename T,typename U>
pair<pair<T,T>,T> unit(T x, T y,U k){
T s=floor(x-1,k);
T t=floor(y,k);
return make_pair(make_pair(s+1,t),1);
if(s==t)return make_pair(make_pair(0,0),-1);
}
pair<pair<ll,ll>,int> jufuku(ll a,ll b,ll c,ll d){
//a<=b,c<=dが保証されているとする
if(a>c){
swap(a,c);
swap(b,d);
}
if(c>b)return make_pair(make_pair(0,0),-1);
return make_pair(make_pair(c,min(b,d)),1);
}
template<typename T>
bool chmin(T &a, T b){
if(a>b){
a=b;
return true;
}
return false;
}
template<typename T>
void her(vector<T> &a){
for(auto &g:a)g--;
}
ll mypow(ll x,ll y,ll MOD){
if(MOD==-1){
MOD=9223372036854775807LL;
}
ll ret=1;
while(y>0){
if(y&1)ret=ret*x%MOD;
x=x*x%MOD;
y>>=1;
}
return ret;
}
struct UnionFind {
vector<int> par,siz,mi,ma;
UnionFind(int n) : par(n,-1), siz(n,1) {mi.resize(n);ma.resize(n);iota(mi.begin(),mi.end(),0);iota(ma.begin(),ma.end(),0); }
int root(int x) {
if(par[x]==-1)return x;
else return par[x]=root(par[x]);
}
bool same(int x, int y) {
return root(x)==root(y);
}
bool marge(int x, int y) {
int rx = root(x), ry = root(y);
if (rx==ry) return false;
if(siz[rx]<siz[ry])swap(rx,ry);
par[ry] = rx;
siz[rx] += siz[ry];
mi[rx]=min(mi[rx],mi[ry]);
ma[rx]=max(ma[rx],ma[ry]);
return true;
}
int size(int x) {
return siz[root(x)];
}
int mini(int x){
return mi[root(x)];
}
int maxi(int x){
return ma[root(x)];
}
};
vector<long long> fac,finv,invs;
// テーブルを作る前処理
void COMinit(int MAX,int MOD) {
fac.resize(MAX);
finv.resize(MAX);
invs.resize(MAX);
fac[0] = fac[1] = 1;
finv[0] = finv[1] = 1;
invs[1] = 1;
for (int i = 2; i < MAX; i++){
fac[i] = fac[i - 1] * i % MOD;
invs[i] = MOD - invs[MOD%i] * (MOD / i) % MOD;
finv[i] = finv[i - 1] * invs[i] % MOD;
}
}
// 二項係数計算
long long binom(int n, int k,int MOD){
if (n < k) return 0;
if (n < 0 || k < 0) return 0;
return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD;
}
#endif
//けんちょんさんの
// modint
//NTTのmulのあれは絶対にやろう！
template<int MOD> struct Fp {
// inner value
long long val;
// constructor
constexpr Fp() : val(0) { }
constexpr Fp(long long v) : val(v % MOD) {
if (val < 0) val += MOD;
}
constexpr long long get() const { return val; }
constexpr int get_mod() const { return MOD; }
// arithmetic operators
constexpr Fp operator + () const { return Fp(*this); }
constexpr Fp operator - () const { return Fp(0) - Fp(*this); }
constexpr Fp operator + (const Fp &r) const { return Fp(*this) += r; }
constexpr Fp operator - (const Fp &r) const { return Fp(*this) -= r; }
constexpr Fp operator * (const Fp &r) const { return Fp(*this) *= r; }
constexpr Fp operator / (const Fp &r) const { return Fp(*this) /= r; }
constexpr Fp& operator += (const Fp &r) {
val += r.val;
if (val >= MOD) val -= MOD;
return *this;
}
constexpr Fp& operator -= (const Fp &r) {
val -= r.val;
if (val < 0) val += MOD;
return *this;
}
constexpr Fp& operator *= (const Fp &r) {
val = val * r.val % MOD;
return *this;
}
constexpr Fp& operator /= (const Fp &r) {
long long a = r.val, b = MOD, u = 1, v = 0;
while (b) {
long long t = a / b;
a -= t * b, swap(a, b);
u -= t * v, swap(u, v);
}
val = val * u % MOD;
if (val < 0) val += MOD;
return *this;
}
constexpr Fp pow(long long n) const {
Fp res(1), mul(*this);
while (n > 0) {
if (n & 1) res *= mul;
mul *= mul;
n >>= 1;
}
return res;
}
constexpr Fp inv() const {
Fp res(1), div(*this);
return res / div;
}
constexpr Fp sqrt(){
if(this->val < 2) return *this;
if(this->pow((MOD-1)>>1).val != 1) return Fp(0) ;
Fp b = 1 , one = 1 ;
while (b.pow((MOD-1) >> 1) == 1) b += one;
long long m = MOD-1, e = 0;
while (m % 2 == 0) m >>= 1, e += 1;
Fp x = this->pow((m - 1) >> 1);
Fp y = (*this) * x * x;
x *= (*this);
Fp z = b.pow(m);
while (y.val != 1) {
int j = 0;
Fp t = y;
while (t != one) j += 1, t *= t;
z = z.pow(1LL << (e-j-1));
x *= z; z *= z; y *= z; e = j;
}
return x;
}
// other operators
constexpr bool operator == (const Fp &r) const {
return this->val == r.val;
}
constexpr bool operator != (const Fp &r) const {
return this->val != r.val;
}
constexpr Fp& operator ++ () {
++val;
if (val >= MOD) val -= MOD;
return *this;
}
constexpr Fp& operator -- () {
if (val == 0) val += MOD;
--val;
return *this;
}
constexpr Fp operator ++ (int) const {
Fp res = *this;
++*this;
return res;
}
constexpr Fp operator -- (int) const {
Fp res = *this;
--*this;
return res;
}
friend constexpr istream& operator >> (istream &is, Fp<MOD> &x) {
is >> x.val;
x.val %= MOD;
if (x.val < 0) x.val += MOD;
return is;
}
friend constexpr ostream& operator << (ostream &os, const Fp<MOD> &x) {
return os << x.val;
}
friend constexpr Fp<MOD> pow(const Fp<MOD> &r, long long n) {
return r.pow(n);
}
friend constexpr Fp<MOD> inv(const Fp<MOD> &r) {
return r.inv();
}
};
namespace NTT {
long long modpow(long long a, long long n, int mod) {
long long res = 1;
while (n > 0) {
if (n & 1) res = res * a % mod;
a = a * a % mod;
n >>= 1;
}
return res;
}
long long modinv(long long a, int mod) {
long long b = mod, u = 1, v = 0;
while (b) {
long long t = a / b;
a -= t * b, swap(a, b);
u -= t * v, swap(u, v);
}
u %= mod;
if (u < 0) u += mod;
return u;
}
int calc_primitive_root(int mod) {
if (mod == 2) return 1;
if (mod == 167772161) return 3;
if (mod == 469762049) return 3;
if (mod == 754974721) return 11;
if (mod == 998244353) return 3;
int divs[20] = {};
divs[0] = 2;
int cnt = 1;
long long x = (mod - 1) / 2;
while (x % 2 == 0) x /= 2;
for (long long i = 3; i * i <= x; i += 2) {
if (x % i == 0) {
divs[cnt++] = i;
while (x % i == 0) x /= i;
}
}
if (x > 1) divs[cnt++] = x;
for (int g = 2;; g++) {
bool ok = true;
for (int i = 0; i < cnt; i++) {
if (modpow(g, (mod - 1) / divs[i], mod) == 1) {
ok = false;
break;
}
}
if (ok) return g;
}
}
int get_fft_size(int N, int M) {
int size_a = 1, size_b = 1;
while (size_a < N) size_a <<= 1;
while (size_b < M) size_b <<= 1;
return max(size_a, size_b) << 1;
}
// number-theoretic transform
template<class mint> void trans(vector<mint> &v, bool inv = false) {
if (v.empty()) return;
int N = (int)v.size();
int MOD = v[0].get_mod();
int PR = calc_primitive_root(MOD);
static bool first = true;
static vector<long long> vbw(30), vibw(30);
if (first) {
first = false;
for (int k = 0; k < 30; ++k) {
vbw[k] = modpow(PR, (MOD - 1) >> (k + 1), MOD);
vibw[k] = modinv(vbw[k], MOD);
}
}
for (int i = 0, j = 1; j < N - 1; j++) {
for (int k = N >> 1; k > (i ^= k); k >>= 1);
if (i > j) swap(v[i], v[j]);
}
for (int k = 0, t = 2; t <= N; ++k, t <<= 1) {
long long bw = vbw[k];
if (inv) bw = vibw[k];
for (int i = 0; i < N; i += t) {
mint w = 1;
for (int j = 0; j < t/2; ++j) {
int j1 = i + j, j2 = i + j + t/2;
mint c1 = v[j1], c2 = v[j2] * w;
v[j1] = c1 + c2;
v[j2] = c1 - c2;
w *= bw;
}
}
}
if (inv) {
long long invN = modinv(N, MOD);
for (int i = 0; i < N; ++i) v[i] = v[i] * invN;
}
}
// for garner
static constexpr int MOD0 = 754974721;
static constexpr int MOD1 = 167772161;
static constexpr int MOD2 = 469762049;
using mint0 = Fp<MOD0>;
using mint1 = Fp<MOD1>;
using mint2 = Fp<MOD2>;
static const mint1 imod0 = 95869806; // modinv(MOD0, MOD1);
static const mint2 imod1 = 104391568; // modinv(MOD1, MOD2);
static const mint2 imod01 = 187290749; // imod1 / MOD0;
// small case (T = mint, long long)
template<class T> vector<T> naive_mul(const vector<T> &A, const vector<T> &B) {
if (A.empty() || B.empty()) return {};
int N = (int)A.size(), M = (int)B.size();
vector<T> res(N + M - 1);
for (int i = 0; i < N; ++i)
for (int j = 0; j < M; ++j)
res[i + j] += A[i] * B[j];
return res;
}
// mul by convolution
template<class mint> vector<mint> mul(const vector<mint> &A, const vector<mint> &B) {
if (A.empty() || B.empty()) return {};
int N = (int)A.size(), M = (int)B.size();
if (min(N, M) < 30) return naive_mul(A, B);
int MOD = A[0].get_mod();
int size_fft = get_fft_size(N, M);
if (MOD == 998244353) {
vector<mint> a(size_fft), b(size_fft), c(size_fft);
for (int i = 0; i < N; ++i) a[i] = A[i];
for (int i = 0; i < M; ++i) b[i] = B[i];
trans(a), trans(b);
vector<mint> res(size_fft);
for (int i = 0; i < size_fft; ++i) res[i] = a[i] * b[i];
trans(res, true);
res.resize(min(N + M - 1,2500000));
return res;
}
vector<mint0> a0(size_fft, 0), b0(size_fft, 0), c0(size_fft, 0);
vector<mint1> a1(size_fft, 0), b1(size_fft, 0), c1(size_fft, 0);
vector<mint2> a2(size_fft, 0), b2(size_fft, 0), c2(size_fft, 0);
for (int i = 0; i < N; ++i)
a0[i] = A[i].val, a1[i] = A[i].val, a2[i] = A[i].val;
for (int i = 0; i < M; ++i)
b0[i] = B[i].val, b1[i] = B[i].val, b2[i] = B[i].val;
trans(a0), trans(a1), trans(a2), trans(b0), trans(b1), trans(b2);
for (int i = 0; i < size_fft; ++i) {
c0[i] = a0[i] * b0[i];
c1[i] = a1[i] * b1[i];
c2[i] = a2[i] * b2[i];
}
trans(c0, true), trans(c1, true), trans(c2, true);
mint mod0 = MOD0, mod01 = mod0 * MOD1;
vector<mint> res(N + M - 1);
for (int i = 0; i < N + M - 1; ++i) {
int y0 = c0[i].val;
int y1 = (imod0 * (c1[i] - y0)).val;
int y2 = (imod01 * (c2[i] - y0) - imod1 * y1).val;
res[i] = mod01 * y2 + mod0 * y1 + y0;
}
return res;
}
};
template<class T> struct BiCoef {
vector<T> fact_, inv_, finv_;
constexpr BiCoef() {}
constexpr BiCoef(int n) noexcept : fact_(n, 1), inv_(n, 1), finv_(n, 1) {
init(n);
}
constexpr void init(int n) noexcept {
fact_.assign(n, 1), inv_.assign(n, 1), finv_.assign(n, 1);
int MOD = fact_[0].get_mod();
for(int i = 2; i < n; i++){
fact_[i] = fact_[i-1] * i;
inv_[i] = -inv_[MOD%i] * (MOD/i);
finv_[i] = finv_[i-1] * inv_[i];
}
}
constexpr T com(int n, int k) const noexcept {
if (n < k || n < 0 || k < 0) return 0;
return fact_[n] * finv_[k] * finv_[n-k];
}
constexpr T fact(int n) const noexcept {
if (n < 0) return 0;
return fact_[n];
}
constexpr T inv(int n) const noexcept {
if (n < 0) return 0;
return inv_[n];
}
constexpr T finv(int n) const noexcept {
if (n < 0) return 0;
return finv_[n];
}
};
// Formal Power Series
template<typename mint> struct FPS : vector<mint> {
using vector<mint>::vector;
// constructor
constexpr FPS(const vector<mint> &r) : vector<mint>(r) {}
// core operator
constexpr FPS pre(int siz) const {
return FPS(begin(*this), begin(*this) + min((int)this->size(), siz));
}
constexpr FPS rev() const {
FPS res = *this;
reverse(begin(res), end(res));
return res;
}
constexpr FPS& normalize() {
while (!this->empty() && this->back() == 0) this->pop_back();
return *this;
}
// basic operator
constexpr FPS operator - () const noexcept {
FPS res = (*this);
for (int i = 0; i < (int)res.size(); ++i) res[i] = -res[i];
return res;
}
constexpr FPS operator + (const mint &v) const { return FPS(*this) += v; }
constexpr FPS operator + (const FPS &r) const { return FPS(*this) += r; }
constexpr FPS operator - (const mint &v) const { return FPS(*this) -= v; }
constexpr FPS operator - (const FPS &r) const { return FPS(*this) -= r; }
constexpr FPS operator * (const mint &v) const { return FPS(*this) *= v; }
constexpr FPS operator * (const FPS &r) const { return FPS(*this) *= r; }
constexpr FPS operator / (const mint &v) const { return FPS(*this) /= v; }
constexpr FPS operator / (const FPS &r) const { return FPS(*this) /= r; }
constexpr FPS operator % (const FPS &r) const { return FPS(*this) %= r; }
constexpr FPS operator << (int x) const { return FPS(*this) <<= x; }
constexpr FPS operator >> (int x) const { return FPS(*this) >>= x; }
constexpr FPS& operator += (const mint &v) {
if (this->empty()) this->resize(1);
(*this)[0] += v;
return *this;
}
constexpr FPS& operator += (const FPS &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->normalize();
}
constexpr FPS& operator -= (const mint &v) {
if (this->empty()) this->resize(1);
(*this)[0] -= v;
return *this;
}
constexpr FPS& operator -= (const FPS &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->normalize();
}
constexpr FPS& operator *= (const mint &v) {
for (int i = 0; i < (int)this->size(); ++i) (*this)[i] *= v;
return *this;
}
constexpr FPS& operator *= (const FPS &r) {
return *this = NTT::mul((*this), r);
}
constexpr FPS& operator*=(vector<pair<int, mint>> g) {
using T=mint;
int n = (*this).size();
auto [d, c] = g.front();
if (d == 0) g.erase(g.begin());
else c = 0;
drep(i, n) {
(*this)[i] *= c;
for (auto &[j, b] : g) {
if (j > i) break;
(*this)[i] += (*this)[i-j] * b;
}
}
return *this;
}
constexpr FPS& operator/=(vector<pair<int, mint>> g) {
using T=mint;
int n = (*this).size();
auto [d, c] = g.front();
assert(d == 0 && c != T(0));
T ic = c.inv();
g.erase(g.begin());
rep(i, n) {
for (auto &[j, b] : g) {
if (j > i) break;
(*this)[i] -= (*this)[i-j] * b;
}
(*this)[i] *= ic;
}
return *this;
}
constexpr FPS& operator /= (const mint &v) {
assert(v != 0);
mint iv = v.inv();
for (int i = 0; i < (int)this->size(); ++i) (*this)[i] *= iv;
return *this;
}
// division, r must be normalized (r.back() must not be 0)
constexpr FPS& operator /= (const FPS &r) {
assert(!r.empty());
assert(r.back() != 0);
this->normalize();
if (this->size() < r.size()) {
this->clear();
return *this;
}
int need = (int)this->size() - (int)r.size() + 1;
*this = (rev().pre(need) * r.rev().inv(need)).pre(need).rev();
return *this;
}
constexpr FPS& operator %= (const FPS &r) {
assert(!r.empty());
assert(r.back() != 0);
this->normalize();
FPS q = (*this) / r;
return *this -= q * r;
}
constexpr FPS& operator <<= (int x) {
FPS res(x, 0);
res.insert(res.end(), begin(*this), end(*this));
return *this = res;
}
constexpr FPS& operator >>= (int x) {
FPS res;
res.insert(res.end(), begin(*this) + x, end(*this));
return *this = res;
}
constexpr mint eval(const mint &v) {
mint res = 0;
for (int i = (int)this->size()-1; i >= 0; --i) {
res *= v;
res += (*this)[i];
}
return res;
}
// advanced operation
// df/dx
constexpr FPS diff() const {
int n = (int)this->size();
FPS res(n-1);
for (int i = 1; i < n; ++i) res[i-1] = (*this)[i] * i;
return res;
}
// \int f dx
constexpr FPS integral() const {
int n = (int)this->size();
FPS res(n+1, 0);
for (int i = 0; i < n; ++i) res[i+1] = (*this)[i] / (i+1);
return res;
}
// inv(f), f[0] must not be 0
constexpr FPS inv(int deg) const {
assert((*this)[0] != 0);
if (deg < 0) deg = (int)this->size();
FPS res({mint(1) / (*this)[0]});
for (int i = 1; i < deg; i <<= 1) {
res = (res + res - res * res * pre(i << 1)).pre(i << 1);
}
res.resize(deg);
return res;
}
constexpr FPS inv() const {
return inv((int)this->size());
}
constexpr FPS inv_sparse(vector<pair<int, mint>> f){
mint f0 = f[0].second ;
mint g0 = mint(1)/f0 ;
f.erase(f.begin()) ;
const int N = (*this).size() ;
FPS g(N) ;
g[0] = g0 ;
for(int n=1; n < N ;n++){
mint rhs = 0 ;
for(auto &&[i, fi]:f){
if(i <= n) rhs -= fi * g[n-i] ;
}
g[n] = rhs * g0 ;
}
return *this = g ;
}
// log(f) = \int f'/f dx, f[0] must be 1
constexpr FPS log(int deg) const {
assert((*this)[0]==1);
FPS res = (diff() * inv(deg)).integral();
res.resize(deg);
return res;
}
constexpr FPS log() const {
return log((int)this->size());
}
constexpr FPS log_sparse(){
FPS f=(*this);
assert(f[0]==1);
const int n=f.size();
const int mod=f[0].get_mod();
vector<pair<int,mint>>f_sub;
for(int i=1;i<n;i++)if(f[i]!=0)f_sub.emplace_back(i,f[i]);
FPS g(n),inv(n);
for(auto[a,c]:f_sub)if(a)g[a]=c;
if(n>1)inv[1]=1;
for(int i=1;i<n;i++){
mint sum=0;
for(auto[a,c]:f_sub)if(a&&a<=i)sum+=c*(i-a)*g[i-a];
if(i>=2){
inv[i]=-inv[mod%i]*(mod/i);
sum*=inv[i];
}
g[i]-=sum;
}
return g;
}
// exp(f), f[0] must be 0
constexpr FPS exp(int deg) const {
assert((*this)[0] == 0);
FPS res(1, 1);
for (int i = 1; i < deg; i <<= 1) {
res = res * (pre(i<<1) - res.log(i<<1) + 1).pre(i<<1);
}
res.resize(deg);
return res;
}
constexpr FPS exp() const {
return exp((int)this->size());
}
constexpr FPS exp_sparse(const vector<pair<int,mint>>& f_sub, int n){
int mod=f_sub[0].second.get_mod();
FPS g(n),inv(n);
g[0]=1;
if(n>1)inv[1]=1;
for(int i=1;i<n;i++){
for(auto[a,c]:f_sub)if(a<=i)g[i]+=c*g[i-a]*a;
if(i>=2){
inv[i]=-inv[mod%i]*(mod/i);
g[i]*=inv[i];
}
}
return g;
}
//pow(f) = exp(e * log f)
constexpr FPS pow(long long e, int deg) const {
if (e == 0) {
FPS res(deg, 0);
res[0] = 1;
return res;
}
long long i = 0;
while (i < (int)this->size() && (*this)[i] == 0) ++i;
if (i == (int)this->size() || i > (deg - 1) / e) return FPS(deg, 0);
mint k = (*this)[i];
FPS res = ((((*this) >> i) / k).log(deg) * e).exp(deg) * mint(k).pow(e) << (e * i);
res.resize(deg);
return res;
}
constexpr FPS pow(long long e) const {
return pow(e, (int)this->size());
}
constexpr FPS pow_sparse(const vector<pair<int,mint>>&f_sub, int n, long long k){
FPS g(n);
if(f_sub.empty()){
if(k==0)g[0]=1;
return g;
}
auto[a0,c0]=f_sub[0];
if(a0>0){
if(k>(n-1)/a0)return FPS(n);
vector<pair<int,mint>>g_sub=f_sub;
for(auto&[a,c]:g_sub)a-=a0;
auto h=pow_sparse(g_sub,n-k*a0,k);
for(int i=0;i<n-k*a0;i++)g[i+k*a0]=h[i];
return g;
}
const long long mod=(*this)[0].get_mod();
FPS inv(n);
inv[1]=1;
for(int i=2;i<n;i++)inv[i]=-inv[mod%i]*(mod/i);
g[0]=c0.pow(k);
mint cef=c0.inv();
for(int i=1;i<n;i++){
for(auto[a,c]:f_sub)if(a&&i-a>=0)g[i]+=c*g[i-a]*(mint(k+1)*a-i);
g[i]*=inv[i]*cef;
}
return g;
}
// sqrt(f), f[0] must be 1
constexpr FPS sqrt_base(int deg) const {
assert((*this)[0] == 1);
mint inv2 = mint(1) / 2;
FPS res(1, 1);
for (int i = 1; i < deg; i <<= 1) {
res = (res + pre(i << 1) * res.inv(i << 1)).pre(i << 1);
for (mint &x : res) x *= inv2;
}
res.resize(deg);
return res;
}
constexpr FPS sqrt_sparse() {
FPS f = *this ;
int n = f.size();
int di = n;
for (int i = n - 1; i >= 0; --i)if (f[i] != 0) di = i;
if (di == n) return f;
if (di & 1) return {};
mint y = f[di];
mint x = y.sqrt();
if (x * x != y) return {};
mint c = mint(1) / y;
FPS g(n - di);
for (int i = 0; i < n - di; ++i) g[i] = f[di + i] * c;
g = g.sqrt_base();
for (int i = 0; i < int(g.size()); ++i) g[i] *= x; g.resize(n);
for (int i = n - 1; i >= 0; --i) {
if (i >= di / 2) g[i] = g[i - di / 2];
else g[i] = 0;
}
return *this = g;
}
constexpr FPS sqrt_base() const {
return sqrt_base((int)this->size());
}
constexpr FPS subset_sum(const vector<int>& s,int t) {
vector<mint> a(t+1);
for(auto g:s)a[g]+=1;
vector<mint> invmemo(t+1);
reps(i,t){
mint k=i;
invmemo[i]=k.inv();
}
FPS g(t+1);
reps(i,t){
for(int j=1;j*i<=t;j++){
if(j&1)g[j*i]+=invmemo[j]*a[i];
else g[j*i]-=invmemo[j]*a[i];
}
}
g=g.exp();
return g;
}
constexpr FPS subset_sum2(const vector<int>& s,int t) {
vector<mint> a(t+1);
for(auto g:s)a[g]+=1;
vector<mint> invmemo(t+1);
reps(i,t){
mint k=i;
invmemo[i]=k.inv();
}
FPS g(t+1);
reps(i,t){
for(int j=1;j*i<=t;j++){
g[j*i]+=invmemo[j]*a[i];
}
}
g=g.exp();
return g;
}
// friend operators
friend constexpr FPS diff(const FPS &f) { return f.diff(); }
friend constexpr FPS integral(const FPS &f) { return f.integral(); }
friend constexpr FPS inv(const FPS &f, int deg) { return f.inv(deg); }
friend constexpr FPS inv(const FPS &f) { return f.inv((int)f.size()); }
friend constexpr FPS log(const FPS &f, int deg) { return f.log(deg); }
friend constexpr FPS log(const FPS &f) { return f.log((int)f.size()); }
friend constexpr FPS exp(const FPS &f, int deg) { return f.exp(deg); }
friend constexpr FPS exp(const FPS &f) { return f.exp((int)f.size()); }
friend constexpr FPS pow(const FPS &f, long long e, int deg) { return f.pow(e, deg); }
friend constexpr FPS pow(const FPS &f, long long e) { return f.pow(e, (int)f.size()); }
friend constexpr FPS sqrt_base(const FPS &f, int deg) { return f.sqrt_base(deg); }
friend constexpr FPS sqrt_base(const FPS &f) { return f.sqrt_base((int)f.size()); }
};
//------------------------------//
// Polynomial, FPS algorithms
//------------------------------//
//高速ゼータ変換
template<typename T>
vector<T> zeta(vector<T> x,int type){
int n=x.size();
for(int b=1;b<n;b<<=1){
for(int j=0;j<n;j++){
switch(type){
case 1:
if(!(b&j))x[j]+=x[b|j];
break;
case 2:
if(!(b&j))x[b|j]+=x[j];
break;
case 3:
if(!(b&j))x[j]-=x[b|j];
break;
case 4:
if(!(b&j))x[b|j]-=x[j];
break;
}
}
}
return x;
}
template<typename mint>
FPS<mint> convolution_or(FPS<mint> &a,FPS<mint> &b){
FPS<mint> A=zeta(a,2),B=zeta(b,2);
A*=B;
return zeta(A,4);
};
template<typename mint>
FPS<mint> subset_convolution(FPS<mint> &a,FPS<mint> &b){
vector<int> pc(max(a.size(),b.size()));
for(int i=1;i<pc.size();i++)pc[i]=pc[i-(i&(-i))]+1;
auto lift=[&](FPS<mint> &a){
vector<FPS<mint>> A(a.size(),FPS<mint>(__builtin_ctz(a.size())+1,0));
for(int i=0;i<a.size();i++){
A[i][pc[i]]=a[i];
}
return A;
};
auto unlift=[&](vector<FPS<mint>> &A){
FPS<mint> a(A.size());
for(int i=0;i<a.size();i++){
a[i]=A[i][pc[i]];
}
return a;
};
vector<FPS<mint>> A=lift(a),B=lift(b);
A=zeta(A,2),B=zeta(B,2);
auto prod=[&](vector<FPS<mint>>& A,vector<FPS<mint>>& B) {
int n = A.size(), d = __builtin_ctz(n);
for(int i=0;i<n;i++)B[i].resize(d+1);
for(int i=0;i<n;i++)A[i].resize(d+1);
for (int i = 0; i < n; i++) {
FPS<mint> c(d+3);
for (int j = 0; j <= d; j++) {
for (int k = 0; k <= d - j; k++) {
c[j + k] +=
A[i][j] *
B[i]
[k];
}
}
A[i].swap(c);
}
return;
};
prod(A,B);
A=zeta(A,4);
return unlift(A);
}
//convolution all
template <typename mint>
FPS<mint> convolution_all(vector<FPS<mint>>& polys) {
if (polys.size() == 0) return {mint(1)};
while (1) {
int n = polys.size();
if (n == 1) break;
int m = (n > 0 ? (n+ 1) / 2 : n / 2);
rep(i, m) {
if (2 * i + 1 == n) {
polys[i] = polys[2 * i];
} else {
polys[i] = polys[2 * i]*polys[2 * i + 1];
}
}
polys.resize(m);
}
return polys[0];
}
//普通の除算としてのFPS
template<typename mint>
pair<FPS<mint>,FPS<mint>> div_fps(FPS<mint> &f,FPS<mint> g){
while(!g.empty()&&g.back().val==0){
g.pop_back();
}
if(g.empty()){
cerr<<"ERROR : g is 0"<<endl;
exit(9);
}
if(f.size()<g.size()){
return {{},f};
}
reverse(all(g));
int n=f.size(),m=g.size();
FPS<mint> invg=g.inv(n-m+1);
reverse(all(f));
FPS<mint> t=f*invg;
t.resize(n-m+1);
FPS<mint> q(n-m+1);
rep(i,n-m+1){
q[n-m-i]=t[i];
}
reverse(all(f));
reverse(all(g));
t=g*q;
FPS<mint> r(m-1);
rep(i,m-1){
r[i]=f[i]-t[i];
}
while(!r.empty()&&r.back().val==0){
r.pop_back();
}
return {q,r};
}
template<typename mint>
FPS<mint> shift(FPS<mint>a,int t,int m=-1){
/*input: f(0),f(1),...f(n-1),t,m
return:f(c),f(c+1),...f(c+m-1)
order;O((n+m)log(n+m))
*/
int n=a.size();
if(m==-1){
return shift(a,t,n);
}else if(m==0){
return FPS<mint>();
}else if(m<n){
auto p=shift(a,t,n);
p.resize(m);
return p;
}else if(m>n){
int z=max(n,m-n);
auto p=shift(a,t,z);
auto q=shift(a,t+z,n);
p.insert(p.end(),q.begin(),q.end());
p.resize(m);
return p;
}
//前計算
FPS<mint> fact(n,1),expx(n,1),expr(n,1),v(n,1);
loop(i,1,n-1)fact[i]=fact[i-1]*i;
expx[n-1]=fact[n-1].inv();
for(int i=n-2;i>=0;--i)expx[i]=expx[i+1]*(i+1);
loop(i,0,n-1)expr[i]=expx[i]*(i%2?-1:1);
loop(i,1,n-1)v[i]=v[i-1]*(t-i+1);
loop(i,0,n-1)v[i]*=expx[i];
// 下降冪表示に変換
loop(i,0,n-1)a[i]*=expx[i];
a*=expr;
a.resize(n);
// shift
loop(i,0,n-1)a[i]*=fact[i];
reverse(all(a));
a*=v;
a.resize(n);
reverse(all(a));
loop(i,0,n-1)a[i]*=expx[i];
// 標本点に変換
a*=expx;
a.resize(n);
loop(i,0,n-1)a[i]*=fact[i];
// loop(i,0,n-1){
// cout<<v[i].val<<" \n"[i==n-1];
// }
return a;
}
// Bostan-Mori
// find [x^N] P(x)/Q(x), O(K log K log N)
// deg(Q(x)) = K, deg(P(x)) < K
template<typename mint> mint BostanMori(const FPS<mint> &P, const FPS<mint> &Q, long long N) {
assert(!P.empty() && !Q.empty());
if (N == 0 || Q.size() == 1) return P[0] / Q[0];
int qdeg = (int)Q.size();
FPS<mint> P2{P}, minusQ{Q};
P2.resize(qdeg - 1);
for (int i = 1; i < (int)Q.size(); i += 2) minusQ[i] = -minusQ[i];
P2 *= minusQ;
FPS<mint> Q2 = Q * minusQ;
FPS<mint> S(qdeg - 1), T(qdeg);
for (int i = 0; i < (int)S.size(); ++i) {
S[i] = (N % 2 == 0 ? P2[i * 2] : P2[i * 2 + 1]);
}
for (int i = 0; i < (int)T.size(); ++i) {
T[i] = Q2[i * 2];
}
return BostanMori(S, T, N >> 1);
}
//no composition
//0~n partition number
template <typename mint>
FPS<mint> partition_number(int N) {
ll M = sqrt(N) + 10;
FPS<mint> f(N + 1);
for(int x=-M;x<M;x++) {
ll d = x * (3 * x - 1) / 2;
if (d > N) continue;
f[d] += (x % 2 == 0 ? 1 : -1);
}
return f.inv();
}
//はい 数式を入れたら 漸化式が返ってくる 例 1 1 2 3 5 8->1 -1 -1(返ってくるのが1,-1,-1なのは注意)2d+1項以上いれないと作動しない（はず）
template <typename mint>
vector<mint> BerlekampMassey(const vector<mint> &s) {
const int N = (int)s.size();
vector<mint> b, c;
b.reserve(N + 1);
c.reserve(N + 1);
b.push_back(mint(1));
c.push_back(mint(1));
mint y = mint(1);
for (int ed = 1; ed <= N; ed++) {
int l = int(c.size()), m = int(b.size());
mint x = 0;
for (int i = 0; i < l; i++) x += c[i] * s[ed - l + i];
b.emplace_back(mint(0));
m++;
if (x == mint(0)) continue;
mint freq = x / y;
if (l < m) {
auto tmp = c;
c.insert(begin(c), m - l, mint(0));
for (int i = 0; i < m; i++) c[m - 1 - i] -= freq * b[m - 1 - i];
b = tmp;
y = x;
} else {
for (int i = 0; i < m; i++) c[l - 1 - i] -= freq * b[m - 1 - i];
}
}
reverse(begin(c), end(c));
return c;
}
template <typename mint>
mint nth_term(long long n, vector<mint> &s) {
FPS<mint> bm = BerlekampMassey<mint>(s);
FPS<mint> ss=s;
ss*=bm;
ss.resize(bm.size()-1);
return BostanMori(ss,bm,n);
}
//no composition halfgcd (f^k mod g) (x^k mod f) (f^k mod (1-x^n))(Π(1-x^a_i)^-1)
using mint=Fp<modb>;
int main(){
cin.tie(nullptr);
ll n,a,b,k;
cin >> n >> a >> b >> k;
vector<mint> s(100);
s[1]=a,s[2]=b;
for(int i=3;i<=99;i++)s[i]=(s[i-1]*s[i-1]+k)/s[i-2];
cout << nth_term(n,s) << "\n";
}