QOJ.ac
QOJ
| ID | Problem | Submitter | Result | Time | Memory | Language | File size | Submit time | Judge time |
|---|---|---|---|---|---|---|---|---|---|
| #110352 | #6513. Expression 3 | ffao | TL | 7ms | 4944kb | C++23 | 71.4kb | 2023-06-01 17:13:53 | 2023-06-01 17:13:55 |
Judging History
answer
#ifdef LOCAL
#define _GLIBCXX_DEBUG 1
#define dbg(...) cerr << "LINE(" << __LINE__ << ") -> [" << #__VA_ARGS__ << "]: [", DBG(__VA_ARGS__)
#else
#define dbg(...) 0
#endif
#if 0
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
template<class T>
using ordered_set = __gnu_pbds::tree<T, __gnu_pbds::null_type, std::less<T>, __gnu_pbds::rb_tree_tag,
__gnu_pbds::tree_order_statistics_node_update>;
#endif
#include <optional>
#include <vector>
#include <list>
#include <map>
#include <set>
#include <queue>
#include <stack>
#include <bitset>
#include <algorithm>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <cstring>
#include <random>
#include <chrono>
#include <cassert>
using namespace std;
#define sz(x) (int)(x).size()
#define all(x) begin(x), end(x)
#define FOR(i,a,b) for (int i = (a); i < (b); ++i)
#define F0R(i,a) FOR(i,0,a)
#define REP(i,n) for(int (i)=0;(i)<(int)(n);(i)++)
#define each(a,x) for (auto& a: x)
#define tcT template<class T
#define tcTU tcT, class U
#define tcTUU tcT, class ...U
template<class T> using V = vector<T>;
template<class T, size_t SZ> using AR = array<T,SZ>;
typedef string str;
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<int> vi;
typedef vector<vi> vvi;
template<typename T, typename U> T &ctmax(T &x, const U &y){ return x = max<T>(x, y); }
template<typename T, typename U> T &ctmin(T &x, const U &y){ return x = min<T>(x, y); }
mt19937 rng((unsigned)chrono::steady_clock::now().time_since_epoch().count());
#define ts to_string
str ts(char c) { return str(1,c); }
str ts(bool b) { return b ? "true" : "false"; }
str ts(const char* s) { return (str)s; }
str ts(str s) { return s; }
str ts(vector<bool> v) { str res = "{"; F0R(i,sz(v)) res += char('0'+v[i]); res += "}"; return res; }
template<size_t SZ> str ts(bitset<SZ> b) { str res = ""; F0R(i,SZ) res += char('0'+b[i]); return res; }
template<class A, class B> str ts(pair<A,B> p);
template<class T> str ts(T v) { bool fst = 1; str res = "{"; for (const auto& x: v) {if (!fst) res += ", "; fst = 0; res += ts(x);} res += "}"; return res;}
template<class A, class B> str ts(pair<A,B> p) {return "("+ts(p.first)+", "+ts(p.second)+")"; }
template<class A> void pr(A x) { cout << ts(x); }
template<class H, class... T> void pr(const H& h, const T&... t) { pr(h); pr(t...); }
void ps() { pr("\n"); }
template<class H, class... T> void ps(const H& h, const T&... t) { pr(h); if (sizeof...(t)) pr(" "); ps(t...); }
void DBG() { cerr << "]" << endl; }
template<class H, class... T> void DBG(H h, T... t) {cerr << ts(h); if (sizeof...(t)) cerr << ", "; DBG(t...); }
tcTU> void re(pair<T,U>& p);
tcT> void re(V<T>& v);
tcT, size_t SZ> void re(AR<T,SZ>& a);
tcT> void re(T& x) { cin >> x; }
void re(double& d) { str t; re(t); d = stod(t); }
void re(long double& d) { str t; re(t); d = stold(t); }
tcTUU> void re(T& t, U&... u) { re(t); re(u...); }
tcTU> void re(pair<T,U>& p) { re(p.first,p.second); }
tcT> void re(V<T>& x) { each(a,x) re(a); }
tcT, size_t SZ> void re(AR<T,SZ>& x) { each(a,x) re(a); }
tcT> void rv(int n, V<T>& x) { x.rsz(n); re(x); }
constexpr bool multitest() {return 0;}
void solve();
int main() {
ios_base::sync_with_stdio(false); cin.tie(NULL);
int t = 1;
if (multitest()) cin >> t;
for (; t; t--) solve();
}
#define rng(i,a,b) for(int i=int(a);i<int(b);i++)
#define rep(i,b) rng(i,0,b)
#define gnr(i,a,b) for(int i=int(b)-1;i>=int(a);i--)
tcT> using vc = vector<T>;
tcT> using vvc=vc<vc<T>>;
using uint=unsigned;
using ull=unsigned long long;
#define si(x) int(x.size())
#define per(i,b) gnr(i,0,b)
#define pb push_back
#define eb emplace_back
#define a first
#define b second
#define bg begin()
#define ed end()
struct modinfo{uint mod,root;};
template<modinfo const&ref>
struct modular{
static constexpr uint const &mod=ref.mod;
static modular root(){return modular(ref.root);}
uint v;
//modular(initializer_list<uint>ls):v(*ls.bg){}
modular(ll vv=0){s(vv%mod+mod);}
modular& s(uint vv){
v=vv<mod?vv:vv-mod;
return *this;
}
modular operator-()const{return modular()-*this;}
modular& operator+=(const modular&rhs){return s(v+rhs.v);}
modular&operator-=(const modular&rhs){return s(v+mod-rhs.v);}
modular&operator*=(const modular&rhs){
v=ull(v)*rhs.v%mod;
return *this;
}
modular&operator/=(const modular&rhs){return *this*=rhs.inv();}
modular operator+(const modular&rhs)const{return modular(*this)+=rhs;}
modular operator-(const modular&rhs)const{return modular(*this)-=rhs;}
modular operator*(const modular&rhs)const{return modular(*this)*=rhs;}
modular operator/(const modular&rhs)const{return modular(*this)/=rhs;}
modular pow(int n)const{
modular res(1),x(*this);
while(n){
if(n&1)res*=x;
x*=x;
n>>=1;
}
return res;
}
modular inv()const{return pow(mod-2);}
};
#define USE_FMT
//size of input must be a power of 2
//output of forward fmt is bit-reversed
//output elements are in the range [0,mod*4)
//input of inverse fmt should be bit-reversed
template<class mint>
void inplace_fmt(const int n,mint*const f,bool inv){
static constexpr uint mod=mint::mod;
static constexpr uint mod2=mod*2;
static const int L=30;
static mint g[L],ig[L],p2[L];
if(g[0].v==0){
rep(i,L){
mint w=-mint::root().pow(((mod-1)>>(i+2))*3);
g[i]=w;
ig[i]=w.inv();
p2[i]=mint(1<<i).inv();
}
}
if(!inv){
int b=n;
if(b>>=1){//input:[0,mod)
rep(i,b){
uint x=f[i+b].v;
f[i+b].v=f[i].v+mod-x;
f[i].v+=x;
}
}
if(b>>=1){//input:[0,mod*2)
mint p=1;
for(int i=0,k=0;i<n;i+=b*2){
rng(j,i,i+b){
uint x=(f[j+b]*p).v;
f[j+b].v=f[j].v+mod-x;
f[j].v+=x;
}
p*=g[__builtin_ctz(++k)];
}
}
while(b){
if(b>>=1){//input:[0,mod*3)
mint p=1;
for(int i=0,k=0;i<n;i+=b*2){
rng(j,i,i+b){
uint x=(f[j+b]*p).v;
f[j+b].v=f[j].v+mod-x;
f[j].v+=x;
}
p*=g[__builtin_ctz(++k)];
}
}
if(b>>=1){//input:[0,mod*4)
mint p=1;
for(int i=0,k=0;i<n;i+=b*2){
rng(j,i,i+b){
uint x=(f[j+b]*p).v;
f[j].v=(f[j].v<mod2?f[j].v:f[j].v-mod2);
f[j+b].v=f[j].v+mod-x;
f[j].v+=x;
}
p*=g[__builtin_ctz(++k)];
}
}
}
}else{
int b=1;
if(b<n/2){//input:[0,mod)
mint p=1;
for(int i=0,k=0;i<n;i+=b*2){
rng(j,i,i+b){
ull x=f[j].v+mod-f[j+b].v;
f[j].v+=f[j+b].v;
f[j+b].v=x*p.v%mod;
}
p*=ig[__builtin_ctz(++k)];
}
b<<=1;
}
for(;b<n/2;b<<=1){
mint p=1;
for(int i=0,k=0;i<n;i+=b*2){
rng(j,i,i+b/2){//input:[0,mod*2)
ull x=f[j].v+mod2-f[j+b].v;
f[j].v+=f[j+b].v;
f[j].v=(f[j].v)<mod2?f[j].v:f[j].v-mod2;
f[j+b].v=x*p.v%mod;
}
rng(j,i+b/2,i+b){//input:[0,mod)
ull x=f[j].v+mod-f[j+b].v;
f[j].v+=f[j+b].v;
f[j+b].v=x*p.v%mod;
}
p*=ig[__builtin_ctz(++k)];
}
}
if(b<n){//input:[0,mod*2)
rep(i,b){
uint x=f[i+b].v;
f[i+b].v=f[i].v+mod2-x;
f[i].v+=x;
}
}
mint z=p2[__lg(n)];
rep(i,n)f[i]*=z;
}
}
template<class mint>
void inplace_fmt(vector<mint>&f,bool inv){
inplace_fmt(si(f),f.data(),inv);
}
template<class mint>
void half_fmt(const int n,mint*const f){
static constexpr uint mod=mint::mod;
static constexpr uint mod2=mod*2;
static const int L=30;
static mint g[L],h[L];
if(g[0].v==0){
rep(i,L){
g[i]=-mint::root().pow(((mod-1)>>(i+2))*3);
h[i]=mint::root().pow((mod-1)>>(i+2));
}
}
int b=n;
int lv=0;
if(b>>=1){//input:[0,mod)
mint p=h[lv++];
for(int i=0,k=0;i<n;i+=b*2){
rng(j,i,i+b){
uint x=(f[j+b]*p).v;
f[j+b].v=f[j].v+mod-x;
f[j].v+=x;
}
p*=g[__builtin_ctz(++k)];
}
}
if(b>>=1){//input:[0,mod*2)
mint p=h[lv++];
for(int i=0,k=0;i<n;i+=b*2){
rng(j,i,i+b){
uint x=(f[j+b]*p).v;
f[j+b].v=f[j].v+mod-x;
f[j].v+=x;
}
p*=g[__builtin_ctz(++k)];
}
}
while(b){
if(b>>=1){//input:[0,mod*3)
mint p=h[lv++];
for(int i=0,k=0;i<n;i+=b*2){
rng(j,i,i+b){
uint x=(f[j+b]*p).v;
f[j+b].v=f[j].v+mod-x;
f[j].v+=x;
}
p*=g[__builtin_ctz(++k)];
}
}
if(b>>=1){//input:[0,mod*4)
mint p=h[lv++];
for(int i=0,k=0;i<n;i+=b*2){
rng(j,i,i+b){
uint x=(f[j+b]*p).v;
f[j].v=(f[j].v<mod2?f[j].v:f[j].v-mod2);
f[j+b].v=f[j].v+mod-x;
f[j].v+=x;
}
p*=g[__builtin_ctz(++k)];
}
}
}
}
template<class mint>
void half_fmt(vector<mint>&f){
half_fmt(si(f),f.data());
}
//59501818244292734739283969-1=5.95*10^25 までの値を正しく計算
//最終的な列の大きさが 2^24 までなら動く
//最終的な列の大きさが 2^20 以下のときは，下の 3 つの素数を使ったほうが速い（は？）
//VERIFY: yosupo
namespace arbitrary_convolution{
//constexpr modinfo base0{167772161,3};//2^25 * 5 + 1
//constexpr modinfo base1{469762049,3};//2^26 * 7 + 1
//constexpr modinfo base2{754974721,11};//2^24 * 45 + 1
constexpr modinfo base0{1045430273,3};//2^20 * 997 + 1
constexpr modinfo base1{1051721729,6};//2^20 * 1003 + 1
constexpr modinfo base2{1053818881,7};//2^20 * 1005 + 1
using mint0=modular<base0>;
using mint1=modular<base1>;
using mint2=modular<base2>;
template<class t,class mint>
vc<t> sub(const vc<mint>&x,const vc<mint>&y,bool same=false){
int n=si(x)+si(y)-1;
int s=1;
while(s<n)s*=2;
vc<t> z(s);rep(i,si(x))z[i]=x[i].v;
inplace_fmt(z,false);
if(!same){
vc<t> w(s);rep(i,si(y))w[i]=y[i].v;
inplace_fmt(w,false);
rep(i,s)z[i]*=w[i];
}else{
rep(i,s)z[i]*=z[i];
}
inplace_fmt(z,true);z.resize(n);
return z;
}
template<class mint>
vc<mint> multiply(const vc<mint>&x,const vc<mint>&y,bool same=false){
auto d0=sub<mint0>(x,y,same);
auto d1=sub<mint1>(x,y,same);
auto d2=sub<mint2>(x,y,same);
int n=si(d0);
vc<mint> res(n);
static const mint1 r01=mint1(mint0::mod).inv();
static const mint2 r02=mint2(mint0::mod).inv();
static const mint2 r12=mint2(mint1::mod).inv();
static const mint2 r02r12=r02*r12;
static const mint w1=mint(mint0::mod);
static const mint w2=w1*mint(mint1::mod);
rep(i,n){
ull a=d0[i].v;
ull b=(d1[i].v+mint1::mod-a)*r01.v%mint1::mod;
ull c=((d2[i].v+mint2::mod-a)*r02r12.v+(mint2::mod-b)*r12.v)%mint2::mod;
res[i].v=(a+b*w1.v+c*w2.v)%mint::mod;
}
return res;
}
}
using arbitrary_convolution::multiply;
template<class mint>
struct Poly:public vc<mint>{
template<class...Args>
Poly(Args...args):vc<mint>(args...){}
Poly(initializer_list<mint>init):vc<mint>(all(init)){}
int size()const{
return vc<mint>::size();
}
void ups(int s){
if(size()<s)this->resize(s,0);
}
Poly low(int s)const{
return Poly(this->bg,this->bg+min(max(s,int(1)),size()));
}
Poly rev()const{
auto r=*this;
reverse(all(r));
return r;
}
mint freq(int i)const{
return i<size()?(*this)[i]:0;
}
Poly operator-()const{
Poly res=*this;
for(auto&v:res)v=-v;
return res;
}
Poly& operator+=(const Poly&r){
ups(r.size());
rep(i,r.size())
(*this)[i]+=r[i];
return *this;
}
template<class T>
Poly& operator+=(T t){
(*this)[0]+=t;
return *this;
}
Poly& operator-=(const Poly&r){
ups(r.size());
rep(i,r.size())
(*this)[i]-=r[i];
return *this;
}
template<class T>
Poly& operator-=(T t){
(*this)[0]-=t;
return *this;
}
Poly& operator*=(const Poly&r){
return *this=multiply(*this,r);
}
Poly operator+(const Poly&r)const{return Poly(*this)+=r;}
Poly operator-(const Poly&r)const{return Poly(*this)-=r;}
Poly operator*(const Poly&r)const{return Poly(*this)*=r;}
};
constexpr modinfo base{998244353,3};
constexpr modinfo base1{998244352,0};
using mint=modular<base>;
using mint1=modular<base1>;
#include <limits>
#include <type_traits>
#undef rng
namespace atcoder {
namespace internal {
// @param m `1 <= m`
// @return x mod m
constexpr long long safe_mod(long long x, long long m) {
x %= m;
if (x < 0) x += m;
return x;
}
// Fast modular multiplication by barrett reduction
// Reference: https://en.wikipedia.org/wiki/Barrett_reduction
// NOTE: reconsider after Ice Lake
struct barrett {
unsigned int _m;
unsigned long long im;
// @param m `1 <= m`
explicit barrett(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1) {}
// @return m
unsigned int umod() const { return _m; }
// @param a `0 <= a < m`
// @param b `0 <= b < m`
// @return `a * b % m`
unsigned int mul(unsigned int a, unsigned int b) const {
// [1] m = 1
// a = b = im = 0, so okay
// [2] m >= 2
// im = ceil(2^64 / m)
// -> im * m = 2^64 + r (0 <= r < m)
// let z = a*b = c*m + d (0 <= c, d < m)
// a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im
// c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2
// ((ab * im) >> 64) == c or c + 1
unsigned long long z = a;
z *= b;
#ifdef _MSC_VER
unsigned long long x;
_umul128(z, im, &x);
#else
unsigned long long x =
(unsigned long long)(((unsigned __int128)(z)*im) >> 64);
#endif
unsigned long long y = x * _m;
return (unsigned int)(z - y + (z < y ? _m : 0));
}
};
// @param n `0 <= n`
// @param m `1 <= m`
// @return `(x ** n) % m`
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 = safe_mod(x, m);
while (n) {
if (n & 1) r = (r * y) % _m;
y = (y * y) % _m;
n >>= 1;
}
return r;
}
// Reference:
// M. Forisek and J. Jancina,
// Fast Primality Testing for Integers That Fit into a Machine Word
// @param n `0 <= n`
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);
// @param b `1 <= b`
// @return pair(g, x) s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g
constexpr std::pair<long long, long long> inv_gcd(long long a, long long b) {
a = safe_mod(a, b);
if (a == 0) return {b, 0};
// Contracts:
// [1] s - m0 * a = 0 (mod b)
// [2] t - m1 * a = 0 (mod b)
// [3] s * |m1| + t * |m0| <= b
long long s = b, t = a;
long long m0 = 0, m1 = 1;
while (t) {
long long u = s / t;
s -= t * u;
m0 -= m1 * u; // |m1 * u| <= |m1| * s <= b
// [3]:
// (s - t * u) * |m1| + t * |m0 - m1 * u|
// <= s * |m1| - t * u * |m1| + t * (|m0| + |m1| * u)
// = s * |m1| + t * |m0| <= b
auto tmp = s;
s = t;
t = tmp;
tmp = m0;
m0 = m1;
m1 = tmp;
}
// by [3]: |m0| <= b/g
// by g != b: |m0| < b/g
if (m0 < 0) m0 += b / s;
return {s, m0};
}
} // namespace internal
} // namespace atcoder
namespace atcoder {
namespace internal {
template <class T>
using is_signed_int128 =
typename std::conditional<std::is_same<T, __int128_t>::value ||
std::is_same<T, __int128>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int128 =
typename std::conditional<std::is_same<T, __uint128_t>::value ||
std::is_same<T, unsigned __int128>::value,
std::true_type,
std::false_type>::type;
template <class T>
using make_unsigned_int128 =
typename std::conditional<std::is_same<T, __int128_t>::value,
__uint128_t,
unsigned __int128>;
template <class T>
using is_integral = typename std::conditional<std::is_integral<T>::value ||
is_signed_int128<T>::value ||
is_unsigned_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_signed_int = typename std::conditional<(is_integral<T>::value &&
std::is_signed<T>::value) ||
is_signed_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int =
typename std::conditional<(is_integral<T>::value &&
std::is_unsigned<T>::value) ||
is_unsigned_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using to_unsigned = typename std::conditional<
is_signed_int128<T>::value,
make_unsigned_int128<T>,
typename std::conditional<std::is_signed<T>::value,
std::make_unsigned<T>,
std::common_type<T>>::type>::type;
template <class T>
using is_signed_int_t = std::enable_if_t<is_signed_int<T>::value>;
template <class T>
using is_unsigned_int_t = std::enable_if_t<is_unsigned_int<T>::value>;
template <class T> using to_unsigned_t = typename to_unsigned<T>::type;
} // namespace internal
} // namespace atcoder
namespace atcoder {
namespace internal {
struct modint_base {};
struct static_modint_base : modint_base {};
template <class T> using is_modint = std::is_base_of<modint_base, T>;
template <class T> using is_modint_t = std::enable_if_t<is_modint<T>::value>;
} // namespace internal
template <int m, std::enable_if_t<(1 <= m)>* = nullptr>
struct static_modint : internal::static_modint_base {
using mint = static_modint;
public:
static constexpr int mod() { return m; }
static mint raw(int v) {
mint x;
x._v = v;
return x;
}
static_modint() : _v(0) {}
template <class T, internal::is_signed_int_t<T>* = nullptr>
static_modint(T v) {
long long x = (long long)(v % (long long)(umod()));
if (x < 0) x += umod();
_v = (unsigned int)(x);
}
template <class T, internal::is_unsigned_int_t<T>* = nullptr>
static_modint(T v) {
_v = (unsigned int)(v % umod());
}
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 = *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 (prime) {
assert(_v);
return pow(umod() - 2);
} else {
auto eg = internal::inv_gcd(_v, m);
assert(eg.first == 1);
return eg.second;
}
}
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;
}
private:
unsigned int _v;
static constexpr unsigned int umod() { return m; }
static constexpr bool prime = internal::is_prime<m>;
};
template <int id> struct dynamic_modint : internal::modint_base {
using mint = dynamic_modint;
public:
static int mod() { return (int)(bt.umod()); }
static void set_mod(int m) {
assert(1 <= m);
bt = internal::barrett(m);
}
static mint raw(int v) {
mint x;
x._v = v;
return x;
}
dynamic_modint() : _v(0) {}
template <class T, internal::is_signed_int_t<T>* = nullptr>
dynamic_modint(T v) {
long long x = (long long)(v % (long long)(mod()));
if (x < 0) x += mod();
_v = (unsigned int)(x);
}
template <class T, internal::is_unsigned_int_t<T>* = nullptr>
dynamic_modint(T v) {
_v = (unsigned int)(v % mod());
}
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 += mod() - rhs._v;
if (_v >= umod()) _v -= umod();
return *this;
}
mint& operator*=(const mint& rhs) {
_v = bt.mul(_v, rhs._v);
return *this;
}
mint& operator/=(const mint& rhs) { return *this = *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 {
auto eg = internal::inv_gcd(_v, mod());
assert(eg.first == 1);
return eg.second;
}
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;
}
private:
unsigned int _v;
static internal::barrett bt;
static unsigned int umod() { return bt.umod(); }
};
template <int id> internal::barrett dynamic_modint<id>::bt(998244353);
using modint998244353 = static_modint<998244353>;
using modint1000000007 = static_modint<1000000007>;
using modint = dynamic_modint<-1>;
namespace internal {
template <class T>
using is_static_modint = std::is_base_of<internal::static_modint_base, T>;
template <class T>
using is_static_modint_t = std::enable_if_t<is_static_modint<T>::value>;
template <class> struct is_dynamic_modint : public std::false_type {};
template <int id>
struct is_dynamic_modint<dynamic_modint<id>> : public std::true_type {};
template <class T>
using is_dynamic_modint_t = std::enable_if_t<is_dynamic_modint<T>::value>;
} // namespace internal
} // namespace atcoder
namespace atcoder {
long long pow_mod(long long x, long long n, int m) {
assert(0 <= n && 1 <= m);
if (m == 1) return 0;
internal::barrett bt((unsigned int)(m));
unsigned int r = 1, y = (unsigned int)(internal::safe_mod(x, m));
while (n) {
if (n & 1) r = bt.mul(r, y);
y = bt.mul(y, y);
n >>= 1;
}
return r;
}
long long inv_mod(long long x, long long m) {
assert(1 <= m);
auto z = internal::inv_gcd(x, m);
assert(z.first == 1);
return z.second;
}
// (rem, mod)
std::pair<long long, long long> crt(const std::vector<long long>& r,
const std::vector<long long>& m) {
assert(r.size() == m.size());
int n = int(r.size());
// Contracts: 0 <= r0 < m0
long long r0 = 0, m0 = 1;
for (int i = 0; i < n; i++) {
assert(1 <= m[i]);
long long r1 = internal::safe_mod(r[i], m[i]), m1 = m[i];
if (m0 < m1) {
std::swap(r0, r1);
std::swap(m0, m1);
}
if (m0 % m1 == 0) {
if (r0 % m1 != r1) return {0, 0};
continue;
}
long long g, im;
std::tie(g, im) = internal::inv_gcd(m0, m1);
long long u1 = (m1 / g);
if ((r1 - r0) % g) return {0, 0};
long long x = (r1 - r0) / g % u1 * im % u1;
r0 += x * m0;
m0 *= u1; // -> lcm(m0, m1)
if (r0 < 0) r0 += m0;
}
return {r0, m0};
}
} // namespace atcoder
namespace suisen {
// ! utility
template <typename ...Types>
using constraints_t = std::enable_if_t<std::conjunction_v<Types...>, std::nullptr_t>;
template <bool cond_v, typename Then, typename OrElse>
constexpr decltype(auto) constexpr_if(Then&& then, OrElse&& or_else) {
if constexpr (cond_v) {
return std::forward<Then>(then);
} else {
return std::forward<OrElse>(or_else);
}
}
// ! function
template <typename ReturnType, typename Callable, typename ...Args>
using is_same_as_invoke_result = std::is_same<std::invoke_result_t<Callable, Args...>, ReturnType>;
template <typename F, typename T>
using is_uni_op = is_same_as_invoke_result<T, F, T>;
template <typename F, typename T>
using is_bin_op = is_same_as_invoke_result<T, F, T, T>;
template <typename Comparator, typename T>
using is_comparator = std::is_same<std::invoke_result_t<Comparator, T, T>, bool>;
// ! integral
template <typename T, typename = constraints_t<std::is_integral<T>>>
constexpr int bit_num = std::numeric_limits<std::make_unsigned_t<T>>::digits;
template <typename T, unsigned int n>
struct is_nbit { static constexpr bool value = bit_num<T> == n; };
template <typename T, unsigned int n>
static constexpr bool is_nbit_v = is_nbit<T, n>::value;
// ?
template <typename T>
struct safely_multipliable {};
template <>
struct safely_multipliable<int> { using type = long long; };
template <>
struct safely_multipliable<long long> { using type = __int128_t; };
template <>
struct safely_multipliable<unsigned int> { using type = unsigned long long; };
template <>
struct safely_multipliable<unsigned long int> { using type = __uint128_t; };
template <>
struct safely_multipliable<unsigned long long> { using type = __uint128_t; };
template <>
struct safely_multipliable<float> { using type = float; };
template <>
struct safely_multipliable<double> { using type = double; };
template <>
struct safely_multipliable<long double> { using type = long double; };
template <typename T>
using safely_multipliable_t = typename safely_multipliable<T>::type;
template <typename T, typename = void>
struct rec_value_type {
using type = T;
};
template <typename T>
struct rec_value_type<T, std::void_t<typename T::value_type>> {
using type = typename rec_value_type<typename T::value_type>::type;
};
template <typename T>
using rec_value_type_t = typename rec_value_type<T>::type;
} // namespace suisen
#include <cstdint>
#include <iterator>
#include <tuple>
namespace suisen {
namespace internal::montgomery {
template <typename Int, typename DInt>
struct Montgomery {
private:
static constexpr uint32_t bits = std::numeric_limits<Int>::digits;
static constexpr Int mask = ~Int(0);
Int N, N2, R2;
// RR = R * R (mod N)
Int RR;
public:
constexpr Montgomery() = default;
explicit constexpr Montgomery(Int N) : N(N), N2(calcN2(N)), R2(calcR2(N, N2)), RR(calcRR(N)) {
assert(N & 1);
}
// @returns t * R (mod N)
constexpr Int make(Int t) const {
return reduce(static_cast<DInt>(t) * RR);
}
// @returns T * R^(-1) (mod N)
constexpr Int reduce(DInt T) const {
DInt m = modR(static_cast<DInt>(modR(T)) * N2);
DInt t = divR(T + m * N);
return t >= N ? t - N : t;
}
constexpr Int add(Int A, Int B) const {
return (A += B) >= N ? A - N : A;
}
constexpr Int sub(Int A, Int B) const {
return (A -= B) < 0 ? A + N : A;
}
constexpr Int mul(Int A, Int B) const {
return reduce(static_cast<DInt>(A) * B);
}
constexpr Int div(Int A, Int B) const {
return reduce(static_cast<DInt>(A) * inv(B));
}
constexpr Int inv(Int A) const; // TODO: Implement
constexpr Int pow(Int A, long long b) const {
Int P = make(1);
for (; b; b >>= 1) {
if (b & 1) P = mul(P, A);
A = mul(A, A);
}
return P;
}
private:
static constexpr Int divR(DInt t) { return t >> bits; }
static constexpr Int modR(DInt t) { return t & mask; }
static constexpr Int calcN2(Int N) {
DInt invN = N; // = N^{-1} (mod 2^2)
for (uint32_t cur_bits = 2; cur_bits < bits; cur_bits *= 2) {
invN = modR(invN * modR(2 - N * invN));
}
assert(modR(N * invN) == 1);
return modR(-invN);
}
static constexpr Int calcR2(Int N, Int N2) {
// R * R2 - N * N2 = 1
// => R2 = (1 + N * N2) / R
return divR(1 + static_cast<DInt>(N) * N2);
}
static constexpr Int calcRR(Int N) {
return -DInt(N) % N;
}
};
} // namespace internal::montgomery
using Montgomery32 = internal::montgomery::Montgomery<uint32_t, uint64_t>;
using Montgomery64 = internal::montgomery::Montgomery<uint64_t, __uint128_t>;
} // namespace suisen
namespace suisen::miller_rabin {
namespace internal {
constexpr uint64_t THRESHOLD_1 = 341531ULL;
constexpr uint64_t BASE_1[]{ 9345883071009581737ULL };
constexpr uint64_t THRESHOLD_2 = 1050535501ULL;
constexpr uint64_t BASE_2[]{ 336781006125ULL, 9639812373923155ULL };
constexpr uint64_t THRESHOLD_3 = 350269456337ULL;
constexpr uint64_t BASE_3[]{ 4230279247111683200ULL, 14694767155120705706ULL, 16641139526367750375ULL };
constexpr uint64_t THRESHOLD_4 = 55245642489451ULL;
constexpr uint64_t BASE_4[]{ 2ULL, 141889084524735ULL, 1199124725622454117ULL, 11096072698276303650ULL };
constexpr uint64_t THRESHOLD_5 = 7999252175582851ULL;
constexpr uint64_t BASE_5[]{ 2ULL, 4130806001517ULL, 149795463772692060ULL, 186635894390467037ULL, 3967304179347715805ULL };
constexpr uint64_t THRESHOLD_6 = 585226005592931977ULL;
constexpr uint64_t BASE_6[]{ 2ULL, 123635709730000ULL, 9233062284813009ULL, 43835965440333360ULL, 761179012939631437ULL, 1263739024124850375ULL };
constexpr uint64_t BASE_7[]{ 2U, 325U, 9375U, 28178U, 450775U, 9780504U, 1795265022U };
template <auto BASE, std::size_t SIZE>
constexpr bool miller_rabin(uint64_t n) {
if (n == 2 or n == 3 or n == 5 or n == 7) return true;
if (n <= 1 or n % 2 == 0 or n % 3 == 0 or n % 5 == 0 or n % 7 == 0) return false;
if (n < 121) return true;
const uint32_t s = __builtin_ctzll(n - 1); // >= 1
const uint64_t d = (n - 1) >> s;
const Montgomery64 mg{ n };
const uint64_t one = mg.make(1), minus_one = mg.make(n - 1);
for (std::size_t i = 0; i < SIZE; ++i) {
uint64_t a = BASE[i] % n;
if (a == 0) continue;
uint64_t Y = mg.pow(mg.make(a), d);
if (Y == one) continue;
for (uint32_t r = 0;; ++r, Y = mg.mul(Y, Y)) {
// Y = a^(d 2^r)
if (Y == minus_one) break;
if (r == s - 1) return false;
}
}
return true;
}
}
template <typename T, std::enable_if_t<std::is_integral_v<T>, std::nullptr_t> = nullptr>
constexpr bool is_prime(T n) {
if constexpr (std::is_signed_v<T>) {
assert(n >= 0);
}
const std::make_unsigned_t<T> n_unsigned = n;
assert(n_unsigned <= std::numeric_limits<uint64_t>::max()); // n < 2^64
using namespace internal;
if (n_unsigned < THRESHOLD_1) return miller_rabin<BASE_1, 1>(n_unsigned);
if (n_unsigned < THRESHOLD_2) return miller_rabin<BASE_2, 2>(n_unsigned);
if (n_unsigned < THRESHOLD_3) return miller_rabin<BASE_3, 3>(n_unsigned);
if (n_unsigned < THRESHOLD_4) return miller_rabin<BASE_4, 4>(n_unsigned);
if (n_unsigned < THRESHOLD_5) return miller_rabin<BASE_5, 5>(n_unsigned);
if (n_unsigned < THRESHOLD_6) return miller_rabin<BASE_6, 6>(n_unsigned);
return miller_rabin<BASE_7, 7>(n_unsigned);
}
} // namespace suisen::miller_rabin
#include <vector>
namespace suisen::internal::sieve {
constexpr std::uint8_t K = 8;
constexpr std::uint8_t PROD = 2 * 3 * 5;
constexpr std::uint8_t RM[K] = { 1, 7, 11, 13, 17, 19, 23, 29 };
constexpr std::uint8_t DR[K] = { 6, 4, 2, 4, 2, 4, 6, 2 };
constexpr std::uint8_t DF[K][K] = {
{ 0, 0, 0, 0, 0, 0, 0, 1 }, { 1, 1, 1, 0, 1, 1, 1, 1 },
{ 2, 2, 0, 2, 0, 2, 2, 1 }, { 3, 1, 1, 2, 1, 1, 3, 1 },
{ 3, 3, 1, 2, 1, 3, 3, 1 }, { 4, 2, 2, 2, 2, 2, 4, 1 },
{ 5, 3, 1, 4, 1, 3, 5, 1 }, { 6, 4, 2, 4, 2, 4, 6, 1 },
};
constexpr std::uint8_t DRP[K] = { 48, 32, 16, 32, 16, 32, 48, 16 };
constexpr std::uint8_t DFP[K][K] = {
{ 0, 0, 0, 0, 0, 0, 0, 8 }, { 8, 8, 8, 0, 8, 8, 8, 8 },
{ 16, 16, 0, 16, 0, 16, 16, 8 }, { 24, 8, 8, 16, 8, 8, 24, 8 },
{ 24, 24, 8, 16, 8, 24, 24, 8 }, { 32, 16, 16, 16, 16, 16, 32, 8 },
{ 40, 24, 8, 32, 8, 24, 40, 8 }, { 48, 32, 16, 32, 16, 32, 48, 8 },
};
constexpr std::uint8_t MASK[K][K] = {
{ 0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80 }, { 0x02, 0x20, 0x10, 0x01, 0x80, 0x08, 0x04, 0x40 },
{ 0x04, 0x10, 0x01, 0x40, 0x02, 0x80, 0x08, 0x20 }, { 0x08, 0x01, 0x40, 0x20, 0x04, 0x02, 0x80, 0x10 },
{ 0x10, 0x80, 0x02, 0x04, 0x20, 0x40, 0x01, 0x08 }, { 0x20, 0x08, 0x80, 0x02, 0x40, 0x01, 0x10, 0x04 },
{ 0x40, 0x04, 0x08, 0x80, 0x01, 0x10, 0x20, 0x02 }, { 0x80, 0x40, 0x20, 0x10, 0x08, 0x04, 0x02, 0x01 },
};
constexpr std::uint8_t OFFSET[K][K] = {
{ 0, 1, 2, 3, 4, 5, 6, 7, },
{ 1, 5, 4, 0, 7, 3, 2, 6, },
{ 2, 4, 0, 6, 1, 7, 3, 5, },
{ 3, 0, 6, 5, 2, 1, 7, 4, },
{ 4, 7, 1, 2, 5, 6, 0, 3, },
{ 5, 3, 7, 1, 6, 0, 4, 2, },
{ 6, 2, 3, 7, 0, 4, 5, 1, },
{ 7, 6, 5, 4, 3, 2, 1, 0, },
};
constexpr std::uint8_t mask_to_index(const std::uint8_t bits) {
switch (bits) {
case 1 << 0: return 0;
case 1 << 1: return 1;
case 1 << 2: return 2;
case 1 << 3: return 3;
case 1 << 4: return 4;
case 1 << 5: return 5;
case 1 << 6: return 6;
case 1 << 7: return 7;
default: assert(false);
}
}
} // namespace suisen::internal::sieve
namespace suisen {
template <unsigned int N>
class SimpleSieve {
private:
static constexpr unsigned int siz = N / internal::sieve::PROD + 1;
static std::uint8_t flag[siz];
public:
SimpleSieve() {
using namespace internal::sieve;
flag[0] |= 1;
unsigned int k_max = (unsigned int) std::sqrt(N + 2) / PROD;
for (unsigned int kp = 0; kp <= k_max; ++kp) {
for (std::uint8_t bits = ~flag[kp]; bits; bits &= bits - 1) {
const std::uint8_t mp = mask_to_index(bits & -bits), m = RM[mp];
unsigned int kr = kp * (PROD * kp + 2 * m) + m * m / PROD;
for (std::uint8_t mq = mp; kr < siz; kr += kp * DR[mq] + DF[mp][mq], ++mq &= 7) {
flag[kr] |= MASK[mp][mq];
}
}
}
}
std::vector<int> prime_list(unsigned int max_val = N) const {
using namespace internal::sieve;
std::vector<int> res { 2, 3, 5 };
res.reserve(max_val / 25);
for (unsigned int i = 0, offset = 0; i < siz and offset < max_val; ++i, offset += PROD) {
for (uint8_t f = ~flag[i]; f;) {
uint8_t g = f & -f;
res.push_back(offset + RM[mask_to_index(g)]);
f ^= g;
}
}
while (res.size() and (unsigned int) res.back() > max_val) res.pop_back();
return res;
}
bool is_prime(const unsigned int p) const {
using namespace internal::sieve;
switch (p) {
case 2: case 3: case 5: return true;
default:
switch (p % PROD) {
case RM[0]: return ((flag[p / PROD] >> 0) & 1) == 0;
case RM[1]: return ((flag[p / PROD] >> 1) & 1) == 0;
case RM[2]: return ((flag[p / PROD] >> 2) & 1) == 0;
case RM[3]: return ((flag[p / PROD] >> 3) & 1) == 0;
case RM[4]: return ((flag[p / PROD] >> 4) & 1) == 0;
case RM[5]: return ((flag[p / PROD] >> 5) & 1) == 0;
case RM[6]: return ((flag[p / PROD] >> 6) & 1) == 0;
case RM[7]: return ((flag[p / PROD] >> 7) & 1) == 0;
default: return false;
}
}
}
};
template <unsigned int N>
std::uint8_t SimpleSieve<N>::flag[SimpleSieve<N>::siz];
template <unsigned int N>
class Sieve {
private:
static constexpr unsigned int base_max = (N + 1) * internal::sieve::K / internal::sieve::PROD;
static unsigned int pf[base_max + internal::sieve::K];
public:
Sieve() {
using namespace internal::sieve;
pf[0] = 1;
unsigned int k_max = ((unsigned int) std::sqrt(N + 1) - 1) / PROD;
for (unsigned int kp = 0; kp <= k_max; ++kp) {
const int base_i = kp * K, base_act_i = kp * PROD;
for (int mp = 0; mp < K; ++mp) {
const int m = RM[mp], i = base_i + mp;
if (pf[i] == 0) {
unsigned int act_i = base_act_i + m;
unsigned int base_k = (kp * (PROD * kp + 2 * m) + m * m / PROD) * K;
for (std::uint8_t mq = mp; base_k <= base_max; base_k += kp * DRP[mq] + DFP[mp][mq], ++mq &= 7) {
pf[base_k + OFFSET[mp][mq]] = act_i;
}
}
}
}
}
bool is_prime(const unsigned int p) const {
using namespace internal::sieve;
switch (p) {
case 2: case 3: case 5: return true;
default:
switch (p % PROD) {
case RM[0]: return pf[p / PROD * K + 0] == 0;
case RM[1]: return pf[p / PROD * K + 1] == 0;
case RM[2]: return pf[p / PROD * K + 2] == 0;
case RM[3]: return pf[p / PROD * K + 3] == 0;
case RM[4]: return pf[p / PROD * K + 4] == 0;
case RM[5]: return pf[p / PROD * K + 5] == 0;
case RM[6]: return pf[p / PROD * K + 6] == 0;
case RM[7]: return pf[p / PROD * K + 7] == 0;
default: return false;
}
}
}
int prime_factor(const unsigned int p) const {
using namespace internal::sieve;
switch (p % PROD) {
case 0: case 2: case 4: case 6: case 8:
case 10: case 12: case 14: case 16: case 18:
case 20: case 22: case 24: case 26: case 28: return 2;
case 3: case 9: case 15: case 21: case 27: return 3;
case 5: case 25: return 5;
case RM[0]: return pf[p / PROD * K + 0] ? pf[p / PROD * K + 0] : p;
case RM[1]: return pf[p / PROD * K + 1] ? pf[p / PROD * K + 1] : p;
case RM[2]: return pf[p / PROD * K + 2] ? pf[p / PROD * K + 2] : p;
case RM[3]: return pf[p / PROD * K + 3] ? pf[p / PROD * K + 3] : p;
case RM[4]: return pf[p / PROD * K + 4] ? pf[p / PROD * K + 4] : p;
case RM[5]: return pf[p / PROD * K + 5] ? pf[p / PROD * K + 5] : p;
case RM[6]: return pf[p / PROD * K + 6] ? pf[p / PROD * K + 6] : p;
case RM[7]: return pf[p / PROD * K + 7] ? pf[p / PROD * K + 7] : p;
default: assert(false);
}
}
/**
* Returns a vector of `{ prime, index }`.
*/
std::vector<std::pair<int, int>> factorize(unsigned int n) const {
assert(0 < n and n <= N);
std::vector<std::pair<int, int>> prime_powers;
while (n > 1) {
int p = prime_factor(n), c = 0;
do { n /= p, ++c; } while (n % p == 0);
prime_powers.emplace_back(p, c);
}
return prime_powers;
}
/**
* Returns the divisors of `n`.
* It is NOT guaranteed that the returned vector is sorted.
*/
std::vector<int> divisors(unsigned int n) const {
assert(0 < n and n <= N);
std::vector<int> divs { 1 };
for (auto [prime, index] : factorize(n)) {
int sz = divs.size();
for (int i = 0; i < sz; ++i) {
int d = divs[i];
for (int j = 0; j < index; ++j) {
divs.push_back(d *= prime);
}
}
}
return divs;
}
};
template <unsigned int N>
unsigned int Sieve<N>::pf[Sieve<N>::base_max + internal::sieve::K];
} // namespace suisen
namespace suisen::fast_factorize {
namespace internal {
template <typename T>
constexpr int floor_log2(T n) {
int i = 0;
while (n) n >>= 1, ++i;
return i - 1;
}
template <typename T, std::enable_if_t<std::is_integral_v<T>, std::nullptr_t> = nullptr>
T pollard_rho(const T n) {
using M = safely_multipliable_t<T>;
const T m = T(1) << (floor_log2(n) / 5);
static std::mt19937_64 rng{std::random_device{}()};
std::uniform_int_distribution<T> dist(0, n - 1);
// const Montgomery64 mg{n};
while (true) {
T c = dist(rng);
auto f = [&](T x) -> T { return (M(x) * x + c) % n; };
T x, y = 2, ys, q = 1, g = 1;
for (T r = 1; g == 1; r <<= 1) {
x = y;
for (T i = 0; i < r; ++i) y = f(y);
for (T k = 0; k < r and g == 1; k += m) {
ys = y;
for (T i = 0; i < std::min(m, r - k); ++i) y = f(y), q = M(q) * (x > y ? x - y : y - x) % n;
g = std::gcd(q, n);
}
}
if (g == n) {
g = 1;
while (g == 1) ys = f(ys), g = std::gcd(x > ys ? x - ys : ys - x, n);
}
if (g < n) {
if (miller_rabin::is_prime(g)) return g;
if (T d = n / g; miller_rabin::is_prime(d)) return d;
return pollard_rho(g);
}
}
}
}
template <typename T, std::enable_if_t<std::is_integral_v<T>, std::nullptr_t> = nullptr>
std::vector<std::pair<T, int>> factorize(T n) {
static constexpr int threshold = 1000000;
static Sieve<threshold> sieve;
std::vector<std::pair<T, int>> res;
if (n <= threshold) {
for (auto [p, q] : sieve.factorize(n)) res.emplace_back(p, q);
return res;
}
if ((n & 1) == 0) {
int q = 0;
do ++q, n >>= 1; while ((n & 1) == 0);
res.emplace_back(2, q);
}
for (T p = 3; p * p <= n; p += 2) {
if (p >= 101 and n >= 1 << 20) {
while (n > 1) {
if (miller_rabin::is_prime(n)) {
res.emplace_back(std::exchange(n, 1), 1);
} else {
p = internal::pollard_rho(n);
int q = 0;
do ++q, n /= p; while (n % p == 0);
res.emplace_back(p, q);
}
}
break;
}
if (n % p == 0) {
int q = 0;
do ++q, n /= p; while (n % p == 0);
res.emplace_back(p, q);
}
}
if (n > 1) res.emplace_back(n, 1);
return res;
}
} // namespace suisen::fast_factorize
namespace suisen {
namespace internal::order_prime_mod {
template <int id>
struct mint64 {
static uint64_t mod() { return _mod; }
static void set_mod(uint64_t new_mod) { mint64<id>::_mod = new_mod; }
mint64() : _val(0) {}
mint64(long long val) : _val(safe_mod(val)) {}
uint64_t val() { return _val; }
friend mint64& operator*=(mint64& x, const mint64& y) {
x._val = __uint128_t(x._val) * y._val % _mod;
return x;
}
friend mint64 operator*(mint64 x, const mint64& y) {
x *= y;
return x;
}
mint64 pow(long long b) const {
assert(b >= 0);
mint64 p = *this, res = 1;
for (; b; b >>= 1) {
if (b & 1) res *= p;
p *= p;
}
return res;
}
private:
static inline uint64_t _mod;
uint64_t _val;
static uint64_t safe_mod(long long val) { return (val %= _mod) < 0 ? val + _mod : val; }
};
}
template <typename T, std::enable_if_t<std::is_integral_v<T>, std::nullptr_t> = nullptr>
struct OrderMod {
using U = std::make_unsigned_t<T>;
OrderMod() = default;
OrderMod(T m) : _mod(m) {
auto factorized = fast_factorize::factorize<T>(_mod);
_is_prime = factorized.size() == 1;
_lambda = _carmichael(factorized);
_phi = _totient(factorized);
for (auto [p, q] : fast_factorize::factorize<T>(_lambda)) {
U r = 1;
for (int i = 0; i < q; ++i) r *= p;
_fac_lambda.emplace_back(p, q, r);
}
}
bool is_primitive_root(U a) const {
if (_mod < 1ULL << 32) {
using mint = atcoder::dynamic_modint<1000000000>;
U old_mod = mint::mod();
mint::set_mod(_mod);
bool res = _is_primitive_root_impl<mint>(a);
mint::set_mod(old_mod);
return res;
} else {
using mint = internal::order_prime_mod::mint64<1000000000>;
U old_mod = mint::mod();
mint::set_mod(_mod);
bool res = _is_primitive_root_impl<mint>(a);
mint::set_mod(old_mod);
return res;
}
}
T primitive_root() const {
assert(_lambda == _phi);
if (_mod < 1ULL << 32) {
return _primitive_root_impl<std::mt19937>();
} else {
return _primitive_root_impl<std::mt19937_64>();
}
}
T operator()(U a) const {
if (_mod < 1ULL << 32) {
using mint = atcoder::dynamic_modint<1000000000>;
U old_mod = mint::mod();
mint::set_mod(_mod);
T res = _order_impl<mint>(a);
mint::set_mod(old_mod);
return res;
} else {
using mint = internal::order_prime_mod::mint64<1000000000>;
U old_mod = mint::mod();
mint::set_mod(_mod);
T res = _order_impl<mint>(a);
mint::set_mod(old_mod);
return res;
}
}
T mod() const {
return _mod;
}
T totient() const {
return _phi;
}
T carmichael() const {
return _lambda;
}
bool is_prime() const {
return _is_prime;
}
std::vector<T> carmichael_prime_factors() const {
std::vector<T> res;
for (const auto &e : _fac_lambda) res.push_back(std::get<0>(e));
return res;
}
private:
U _mod;
U _phi;
U _lambda;
bool _is_prime;
std::vector<std::tuple<U, int, U>> _fac_lambda;
static T _carmichael(const std::vector<std::pair<T, int>>& factorized) {
T lambda = 1;
for (auto [p, ep] : factorized) {
T phi = p - 1;
int exponent = ep - (1 + (p == 2 and ep >= 3));
for (int i = 0; i < exponent; ++i) phi *= p;
lambda = std::lcm(lambda, phi);
}
return lambda;
}
static T _totient(const std::vector<std::pair<T, int>>& factorized) {
T t = 1;
for (const auto& [p, ep] : factorized) {
t *= p - 1;
for (int i = 0; i < ep - 1; ++i) t *= p;
}
return t;
}
template <typename mint>
bool _is_primitive_root_impl(U a) const {
if (_lambda != _phi) return false;
if (_mod == 2) return a % 2 == 1;
const int k = _fac_lambda.size();
U x = _lambda;
for (const auto& [p, q, pq] : _fac_lambda) x /= p;
mint b = mint(a).pow(x);
if (k == 1) return b.val() != 1;
auto dfs = [&](auto dfs, const int l, const int r, const mint val) -> bool {
const int m = (l + r) >> 1;
U lp = 1;
for (int i = m; i < r; ++i) lp *= std::get<0>(_fac_lambda[i]);
mint lval = val.pow(lp);
if (m - l == 1) {
if (lval.val() == 1) return false;
} else {
if (not dfs(dfs, l, m, lval)) return false;
}
U rp = 1;
for (int i = l; i < m; ++i) rp *= std::get<0>(_fac_lambda[i]);
mint rval = val.pow(rp);
if (r - m == 1) {
if (rval.val() == 1) return false;
} else {
if (not dfs(dfs, m, r, rval)) return false;
}
return true;
};
return dfs(dfs, 0, k, b);
}
template <typename Rng>
T _primitive_root_impl() const {
if (_mod == 2) return 1;
Rng rng{ std::random_device{}() };
while (true) {
U a = rng() % (_mod - 2) + 2;
while (not _is_prime and std::gcd(a, _mod) != 1) {
a = rng() % (_mod - 2) + 2;
}
if (is_primitive_root(a)) return a;
}
}
template <typename mint>
U _order_impl(U a) const {
if (_mod == 2) return a % 2 == 1;
const int k = _fac_lambda.size();
U res = 1;
auto update = [&](U p, mint val) {
while (val.val() != 1) {
val = val.pow(p);
res *= p;
}
};
if (k == 1) {
update(std::get<0>(_fac_lambda.front()), a);
return res;
}
auto dfs = [&](auto dfs, const int l, const int r, const mint val) -> void {
const int m = (l + r) >> 1;
U lp = 1;
for (int i = m; i < r; ++i) lp *= std::get<2>(_fac_lambda[i]);
mint lval = val.pow(lp);
if (m - l == 1) {
update(std::get<0>(_fac_lambda[l]), lval);
} else {
dfs(dfs, l, m, lval);
}
U rp = 1;
for (int i = l; i < m; ++i) rp *= std::get<2>(_fac_lambda[i]);
mint rval = val.pow(rp);
if (r - m == 1) {
update(std::get<0>(_fac_lambda[m]), rval);
} else {
dfs(dfs, m, r, rval);
}
};
dfs(dfs, 0, k, a);
return res;
}
};
} // namespace suisen
#include <unordered_map>
namespace suisen {
namespace internal::discrete_logarithm {
int floor_log2(int m) {
int res = 0;
while (1 << (res + 1) <= m) ++res;
return res;
}
}
template <typename mint>
int discrete_logarithm_coprime(mint x, mint y) {
const int m = mint::mod();
if (y == 1 or m == 1) return 0;
if (x == 0) return y == 0 ? 1 : -1;
const int p = ::sqrt(m) + 1;
mint a = mint(x).inv().pow(p);
std::unordered_map<int, int> mp;
mint pow_x = 1;
for (int j = 0; j < p; ++j, pow_x *= x) mp.try_emplace(pow_x.val(), j);
mint ya = y;
for (int i = 0; i < p; ++i, ya *= a) {
if (auto it = mp.find(ya.val()); it != mp.end()) return p * i + it->second;
}
return -1;
}
template <typename mint>
int discrete_logarithm(mint x, mint y) {
using namespace internal::discrete_logarithm;
const int m = mint::mod();
if (m == 1) return 0;
const int d = floor_log2(m);
mint pow_x = 1;
for (int i = 0; i < d; ++i, pow_x *= x) if (pow_x == y) return i;
const int g = std::gcd(pow_x.val(), m);
if (y.val() % g) return -1;
mint::set_mod(m / g);
const int t = discrete_logarithm_coprime(mint(x.val()), mint(y.val()) * mint(pow_x.val()).inv());
mint::set_mod(m);
return t != -1 ? d + t : -1;
}
} // namespace suisen
namespace suisen {
namespace internal::index_calculus {
template <typename mint>
struct SystemOfLinearEquations {
SystemOfLinearEquations() = default;
SystemOfLinearEquations(int target_var_id) : target_var_id(target_var_id) {}
void append(std::vector<mint> r, mint bi) {
int ti = target_var_id + 1;
for (int i = 0; i <= target_var_id; ++i) if (r[i] != 0) {
ti = i;
break;
}
if (ti > target_var_id) return;
const int l = A.size();
for (int i = 0; i < l; ++i) {
eliminate(top_pos[i], A[i], b[i], ti, r, bi);
if (ti > target_var_id) return;
}
A.push_back(std::move(r));
top_pos.push_back(ti);
b.push_back(bi);
}
void append(const std::vector<std::pair<int, mint>>& ri, mint bi) {
std::vector<mint> r(target_var_id + 1);
for (auto [j, v] : ri) r[j] = v;
append(r, bi);
}
std::optional<mint> target() {
if (top_pos.empty() or top_pos.back() != target_var_id) return std::nullopt;
auto [g, inv_c] = atcoder::internal::inv_gcd(A.back().back().val(), mint::mod());
return g == 1 ? std::make_optional(b.back() * inv_c) : std::nullopt;
}
void change_target() {
const int l = A.size();
int t = -1;
for (int i = l - 1; i >= 0; --i) {
if (A[i][target_var_id] != 0) {
t = i;
break;
}
}
SystemOfLinearEquations<mint> nxt(target_var_id);
for (int i = l - 1; i > t; --i) {
nxt.append(A[i], b[i]);
}
for (int i = t - 1; i >= 0; --i) {
if (A[i][target_var_id] != 0) {
const auto T = euclid(A[t][target_var_id].val(), A[i][target_var_id].val());
for (int col = top_pos[i]; col <= target_var_id; ++col) {
apply_euclid(T, A[t][col], A[i][col]);
}
apply_euclid(T, b[t], b[i]);
}
nxt.append(A[i], b[i]);
}
*this = std::move(nxt);
}
private:
int target_var_id;
std::vector<std::vector<mint>> A;
std::vector<int> top_pos;
std::vector<mint> b;
// Euclidean Algorithm
// A * [a, b] = [gcd(a,b), 0]
static std::array<std::array<int, 2>, 2> euclid(int a, int b) {
int x = 1, y = 0, z = 0, w = 1;
while (b) {
int p = a / b, q = a % b;
x -= p * z, std::swap(x, z);
y -= p * w, std::swap(y, w);
a = b, b = q;
}
return { x, y, z, w };
}
static void apply_euclid(const std::array<std::array<int, 2>, 2>& A, mint& x, mint& y) {
mint x2 = x, y2 = y;
x = A[0][0] * x2 + A[0][1] * y2;
y = A[1][0] * x2 + A[1][1] * y2;
}
static void eliminate(int& ti, std::vector<mint>& ri, mint& bi, int& tj, std::vector<mint>& rj, mint& bj) {
if (ti != tj) {
if (ti > tj) std::swap(ti, tj), std::swap(ri, rj), std::swap(bi, bj);
return;
}
const int siz = ri.size();
assert(int(rj.size()) == siz);
const auto T = euclid(ri[ti].val(), rj[tj].val());
for (int col = ti; col < siz; ++col) {
apply_euclid(T, ri[col], rj[col]);
}
apply_euclid(T, bi, bj);
while (tj < siz and rj[tj] == 0) ++tj;
}
};
}
// mod prime only.
struct IndexCalculus {
private:
using mint = atcoder::dynamic_modint<73495793>;
using mint2 = atcoder::dynamic_modint<73495794>;
public:
IndexCalculus() = default;
IndexCalculus(const OrderMod<int>& ord) : _p(ord.mod()), _ord(ord) {};
int discrete_log(int a, int b, int p) {
int old_mod = mint::mod(), old_mod2 = mint2::mod();
mint::set_mod(p), mint2::set_mod(p - 1);
int res = discrete_log_impl(a, b, p);
mint::set_mod(old_mod), mint2::set_mod(old_mod2);
return res;
}
private:
static constexpr uint32_t primes[60]{
2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
};
// inv_primes[i] = ceil(2^32 / primes[i])
static constexpr uint32_t inv_primes[60]{
2147483648, 1431655766, 858993460, 613566757, 390451573, 330382100, 252645136, 226050911, 186737709, 148102321,
138547333, 116080198, 104755300, 99882961, 91382283, 81037119, 72796056, 70409300, 64103990, 60492498,
58835169, 54366675, 51746594, 48258060, 44278014, 42524429, 41698712, 40139882, 39403370, 38008561,
33818641, 32786010, 31350127, 30899046, 28825284, 28443493, 27356480, 26349493, 25718368, 24826401,
23994231, 23729102, 22486740, 22253717, 21801865, 21582751, 20355296, 19259944, 18920561, 18755316,
18433337, 17970575, 17821442, 17111424, 16711936, 16330675, 15966422, 15848588, 15505298, 15284582,
};
int _p = -1;
OrderMod<int> _ord;
int _a = -1;
internal::index_calculus::SystemOfLinearEquations<mint2> _eq;
static int safe_mod(int a, int p) { return ((a %= p) < 0) ? a + p : a; }
static int discrete_log_naive(mint a, mint b) {
mint pow_a = 1;
for (int res = 0; res < int(mint::mod()) - 1; ++res) {
if (pow_a == b) return res;
pow_a *= a;
}
return -1;
}
static std::optional<std::vector<std::pair<int, mint2>>> try_factorize(const int B, uint32_t v) {
assert(B <= 60);
std::vector<std::pair<int, mint2>> res;
for (int i = 0; i < B; ++i) {
const uint32_t m = primes[B - i - 1], im = inv_primes[B - i - 1];
if (uint32_t q = (uint64_t(v) * im) >> 32; v == q * m) {
// q = floor(v/m)
int c = 0;
do v = q, q = (uint64_t(v) * im) >> 32, ++c; while (v == q * m);
res.emplace_back(i, mint2::raw(c));
}
}
return v == 1 ? std::make_optional(std::move(res)) : std::nullopt;
}
int discrete_log_impl(int a, int b, int p) {
a = safe_mod(a, p), b = safe_mod(b, p);
if (b == 0) return a == 0 ? 1 : -1;
if (b == 1) return 0;
if (a == 0) return -1;
if (a == 1) return -1;
if (p <= 64) {
return discrete_log_naive(a, b);
}
const int B = ::pow(p - 1, 1 / (2 * ::sqrt(::log(p - 1) / ::log(::log(p - 1)))));
assert(B <= 60);
if (p != _p) {
_p = p, _ord = OrderMod<int>(p);
_a = a, _eq = internal::index_calculus::SystemOfLinearEquations<mint2>(B);
} else if (a != _a) {
_a = a, _eq = internal::index_calculus::SystemOfLinearEquations<mint2>(B);
} else {
_eq.change_target();
}
assert(_ord.is_prime());
const int ord_a = _ord(a), ord_b = _ord(b);
if (ord_a % ord_b) return -1;
static std::mt19937 rng{};
while (true) {
const int k = rng() % ord_a;
auto opt_fac = try_factorize(B, (mint(a).pow(k) * b).val());
if (not opt_fac) continue;
auto& fac = *opt_fac;
fac.emplace_back(B, -1);
_eq.append(std::move(fac), k);
if (auto opt_res = _eq.target()) return opt_res->val() % _ord(a);
}
}
};
namespace internal::pohlig_hellman {
using mint = atcoder::dynamic_modint<73495793>;
int naive(int a, int b, int n) {
mint pow_a = 1;
for (int res = 0; res < n; ++res) {
if (pow_a == b) return res;
pow_a *= a;
}
return -1;
}
int pohlig_hellman_prime_power_bsgs(int a, int b, int n, int p, int q) {
const mint gamma = mint(a).pow(n / p);
const mint inv_a = mint(a).inv();
int x = 0;
for (int k = 0, pl = 1, pr = n; k < q; ++k) {
pr /= p;
int d = suisen::discrete_logarithm_coprime(gamma, (inv_a.pow(x) * b).pow(pr));
if (d == -1) return -1;
x += pl * d;
pl *= p;
}
return x;
}
int pohlig_hellman_prime_power_index_calculus(int a, int b, int n, int p, int q, const OrderMod<int>& ord) {
const mint gamma = mint(a).pow(n / p);
const mint inv_a = mint(a).inv();
IndexCalculus ic{ ord };
int x = 0;
for (int k = 0, pl = 1, pr = n; k < q; ++k) {
pr /= p;
int d = ic.discrete_log(gamma.val(), (inv_a.pow(x) * b).pow(pr).val(), ord.mod());
x += pl * d;
pl *= p;
}
return x;
}
int pohlig_hellman(int a, int b, int p, int q, const OrderMod<int> &ord) {
// gcd(a, p) = 1
if (b == 1) return 0;
if (a == 1) return b == 1 ? 0 : -1;
if (b % p == 0) return -1;
if (q == 1) return IndexCalculus{ ord }.discrete_log(a, b, p);
const int pq = ord.mod();
const int n = ord(a);
if (n <= 64) return naive(a, b, pq);
if (n % ord(b)) return -1;
const std::vector<std::pair<int, int>> fac_n = [&] {
std::vector<std::pair<int, int>> res;
int n2 = n;
for (int p : ord.carmichael_prime_factors()) if (n2 % p == 0) {
int& c = res.emplace_back(p, 0).second;
do n2 /= p, ++c; while (n2 % p == 0);
}
return res;
}();
std::vector<long long> rs, ms;
for (auto [p, q] : fac_n) {
int pq = 1;
for (int i = 0; i < q; ++i) pq *= p;
int na = mint(a).pow(n / pq).val(), nb = mint(b).pow(n / pq).val();
int x;
if (pq <= 64) {
x = naive(na, nb, pq);
} else if (ord.is_prime()) {
x = pohlig_hellman_prime_power_index_calculus(na, nb, pq, p, q, ord);
} else {
x = pohlig_hellman_prime_power_bsgs(na, nb, pq, p, q);
}
if (x == -1) return -1;
rs.push_back(x), ms.push_back(pq);
}
return atcoder::crt(rs, ms).first;
}
int discrete_logarithm_crt(int a, int b, int m) {
assert(std::gcd(a, m) == 1);
long long r = 0, md = 1;
for (auto [p, q] : fast_factorize::factorize(m)) {
int pq = 1;
for (int i = 0; i < q; ++i) pq *= p;
mint::set_mod(pq);
OrderMod ord(pq);
const long long r2 = pohlig_hellman(a % pq, b % pq, p, q, ord);
if (r2 == -1) return -1;
const long long md2 = ord(a);
std::tie(r, md) = atcoder::crt({ r, r2 }, { md, md2 });
if (md == 0) return -1;
}
return r;
}
int floor_log2(int m) {
int res = 0;
while (1 << (res + 1) <= m) ++res;
return res;
}
int discrete_logarithm(int a, int b, int m) {
if (m == 1 or b == 1) return 0;
if (a == 0) return b == 0 ? 1 : b == 1 ? 0 : -1;
if (a == 1) return b == 1 ? 0 : -1;
const int d = floor_log2(m);
mint pow_x = 1;
for (int i = 0; i < d; ++i, pow_x *= a) if (pow_x == b) return i;
const int g = std::gcd(pow_x.val(), m);
if (b % g) return -1;
mint::set_mod(m / g);
const int t = discrete_logarithm_crt(mint(a).val(), (mint(b) / pow_x.val()).val(), m / g);
return t != -1 ? d + t : -1;
}
}
int discrete_logarithm_fast(int a, int b, int m) {
if ((a %= m) < 0) a += m;
if ((b %= m) < 0) b += m;
int old_mod = internal::pohlig_hellman::mint::mod();
internal::pohlig_hellman::mint::set_mod(m);
int res = internal::pohlig_hellman::discrete_logarithm(a, b, m);
internal::pohlig_hellman::mint::set_mod(old_mod);
return res;
}
} // namespace suisen
str ts(mint s) { return ts(s.v); }
str ts(mint1 s) { return ts(s.v); }
void solve() {
int n; re(n);
vector<mint> v(n);
for (int i = 0; i < n; i++) {
int x; re(x); v[i] = x;
}
string s;
re(s);
Poly<mint1> stars(n-1);
for (int i = 0; i < n-1; i++) stars[i]=(s[i]=='-');
Poly<mint1> dlog(n-1);
for (int i = 1; i < n; i++) dlog[i-1] = -suisen::discrete_logarithm_fast(3,i,base.mod);
for (int i = n-2; i >= 2; i--) dlog[i] -= dlog[i-2];
dlog[0] += suisen::discrete_logarithm_fast(3,base.mod-1,base.mod);
Poly<mint1> den = stars * dlog;
mint ans = 0;
for (int i = 0; i < n; i++) {
mint prod;
if (i>=2 && s[i-2]=='-') prod=0;
else if (i == 0) prod = 1;
else prod = mint(3).pow(den[i-1].v);
//if (i>=1) dbg(i, prod, den[i-1]);
mint prob = (prod + 1) / 2;
ans += (prob*2 - 1) * v[i];
}
for (int i = 1; i <= n-1; i++) {
ans *= mint(i);
}
ps(ans.v);
}
Details
Tip: Click on the bar to expand more detailed information
Test #1:
score: 100
Accepted
time: 6ms
memory: 4944kb
input:
4 9 1 4 1 -+-
output:
46
result:
ok 1 number(s): "46"
Test #2:
score: 0
Accepted
time: 7ms
memory: 4868kb
input:
5 1 2 3 4 5 +-+-
output:
998244313
result:
ok 1 number(s): "998244313"
Test #3:
score: -100
Time Limit Exceeded
input:
100000 664815434 205025136 871445392 797947979 379688564 336946672 231295524 401655676 526374414 670533644 156882283 372427821 700299596 166140732 677498490 44858761 185182210 559696133 813911251 842364231 681916958 114039865 222372111 784286397 437994571 152137641 650875922 613727135 209302742 5321...