結果
| 問題 | No.1962 Not Divide |
| ユーザー |
 PCTprobability
|
| 提出日時 | 2022-03-04 15:50:48 |
| 言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 196 ms / 2,000 ms |
| コード長 | 16,934 bytes |
| コンパイル時間 | 5,950 ms |
| コンパイル使用メモリ | 278,892 KB |
| 最終ジャッジ日時 | 2025-01-28 04:37:10 |
|
ジャッジサーバーID （参考情報） |
judge3 / judge2 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 3 |
| other | AC * 21 |
ソースコード
#include <bits/stdc++.h>
using namespace std;
#if __has_include(<atcoder/all>)
#include <atcoder/all>
using namespace atcoder;
using mint = modint998244353;
#endif
using ll = long long;
using ld = long double;
using ull = unsigned long long;
#define endl "\n"
typedef pair<int,int> Pii;
#define REP(i, n) for (int i = 0; i < (n)a; ++i)
#define REP3(i, m, n) for (int i = (m); (i) < int(n); ++ (i))
#define FOR(i,a,b) for(ll i=a;i<=(ll)(b);i++)
#define rep(i,a,b) for(int i=(int)(a);i<(int)(b);i++)
#define ALL(x) begin(x), end(x)
#define PB push_back
#define rrep(i,a,b) for(int i=a;i>=b;i--)
#define fore(i,a) for(auto &i:a)
#define all(s) (s).begin(),(s).end()
#define drep2(i, m, n) for (int i = (m)-1; i >= (n); --i)
#define drep(i, n) drep2(i, n, 0)
#define rever(vec) reverse(vec.begin(), vec.end())
#define sor(vec) sort(vec.begin(), vec.end())
#define fi first
#define se second
#define pb push_back
#define P pair<ll,ll>
#define PQminll priority_queue<ll, vector<ll>, greater<ll>>
#define PQmaxll priority_queue<ll,vector<ll>,less<ll>>
#define PQminP priority_queue<P, vector<P>, greater<P>>
#define PQmaxP priority_queue<P,vector<P>,less<P>>
#define NP next_permutation
//typedef string::const_iterator State;
//class ParseError {};
//const ll mod = 1000000009;
//const ll mod = 998244353;
const ll mod = 1000000007;
const ll inf = 4100000000000000000ll;
const ld eps = ld(0.000000001);
const long double pi = 3.141592653589793;
template<class T>void vcin(vector<T> &n){for(int i=0;i<int(n.size());i++) cin>>n[i];}
template<class T,class K>void vcin(vector<T> &n,vector<K> &m){for(int i=0;i<int(n.size());i++) cin>>n[i]>>m[i];}
template<class T>void vcout(vector<T> &n){for(int i=0;i<int(n.size());i++){cout<<n[i]<<" ";}cout<<endl;}
template<class T>void vcin(vector<vector<T>> &n){for(int i=0;i<int(n.size());i++){for(int j=0;j<int(n[i].size());j++){cin>>n[i][j];}}}
template<class T>void vcout(vector<vector<T>> &n){for(int i=0;i<int(n.size());i++){for(int j=0;j<int(n[i].size());j++){cout<<n[i][j]<<" ";}cout<<endl;}cout<<endl;}
void yes(bool a){cout<<(a?"yes":"no")<<endl;}
void YES(bool a){cout<<(a?"YES":"NO")<<endl;}
void Yes(bool a){cout<<(a?"Yes":"No")<<endl;}
void possible(bool a){ cout<<(a?"possible":"impossible")<<endl; }
void Possible(bool a){ cout<<(a?"Possible":"Impossible")<<endl; }
void POSSIBLE(bool a){ cout<<(a?"POSSIBLE":"IMPOSSIBLE")<<endl; }
template<class T>void print(T a){cout<<a<<endl;}
template<class T,class F>void print(pair<T,F> a){cout<<a.fi<<" "<<a.se<<endl;}
template<class T>bool chmax(T &a, const T &b) { if (a<b) { a=b; return 1; } return 0;}
template<class T>bool chmin(T &a, const T &b) { if (b<a) { a=b; return 1; } return 0;}
template<class T> void ifmin(T t,T u){if(t>u){cout<<-1<<endl;}else{cout<<t<<endl;}}
template<class T> void ifmax(T t,T u){if(t>u){cout<<-1<<endl;}else{cout<<t<<endl;}}
ll fastgcd(ll u,ll v){ll shl=0;while(u&&v&&u!=v){bool eu=!(u&1);bool ev=!(v&1);if(eu&&ev){++shl;u>>=1;v>>=1;}else if(eu&&!ev){u>>=1;}else if(!eu&&ev){v>>=1;}else if(u>=v){u=(u-v)>>1;}else{ll tmp=u;u=(v-u)>>1;v=tmp;}}return !u?v<<shl:u<<shl;}
ll modPow(ll a, ll n, ll mod) { if(mod==1) return 0;ll ret = 1; ll p = a % mod; while (n) { if (n & 1) ret = ret * p % mod; p = p * p % mod; n >>= 1; } return ret; }
vector<ll> divisor(ll x){ vector<ll> ans; for(ll i = 1; i * i <= x; i++){ if(x % i == 0) {ans.push_back(i); if(i*i!=x){ ans.push_back(x / ans[i]);}}}sor(ans); return ans; }
ll pop(ll x){return __builtin_popcountll(x);}
ll poplong(ll x){ll y=0;while(x){x/=2;y++;}return y;}
void cincout() {
ios::sync_with_stdio(false);
std::cin.tie(nullptr);
cout<< fixed << setprecision(20);
}
template<class T>
vector<T> NTT(vector<T> a,vector<T> b){
ll nmod=T::mod();
int n=a.size();
int m=b.size();
vector<int> x1(n);
vector<int> y1(m);
for(int i=0;i<n;i++){
ll tmp1,tmp2,tmp3;
tmp1=a[i].val();
x1[i]=tmp1;
}
for(int i=0;i<m;i++){
ll tmp1,tmp2,tmp3;
tmp1=b[i].val();
y1[i]=tmp1;
}
auto z1=convolution<167772161>(x1,y1);
auto z2=convolution<469762049>(x1,y1);
auto z3=convolution<1224736769>(x1,y1);
vector<T> res(n+m-1);
ll m1=167772161;
ll m2=469762049;
ll m3=1224736769;
ll m1m2=104391568;
ll m1m2m3=721017874;
ll mm12=m1*m2%nmod;
for(int i=0;i<n+m-1;i++){
int v1=(z2[i]-z1[i])*m1m2%m2;
if(v1<0) v1+=m2;
int v2=(z3[i]-(z1[i]+v1*m1)%m3)*m1m2m3%m3;
if(v2<0) v2+=m3;
res[i]=(z1[i]+v1*m1+v2*mm12);
}
return res;
}
template<class T>
struct FormalPowerSeries:vector<T>{
using vector<T>::vector;
using F=FormalPowerSeries;
F &operator=(const vector<T> &g){
int n=g.size();
int m=(*this).size();
(*this).resize(n);
for(int i=0;i<n;i++) (*this)[i]=g[i];
return (*this);
}
F &operator=(const F &g){
int n=g.size();
int m=(*this).size();
(*this).resize(n);
for(int i=0;i<n;i++) (*this)[i]=g[i];
return (*this);
}
F &operator-(){
for(int i=0;i<(*this).size();i++) (*this)[i]*=-1;
return (*this);
}
F &operator+=(const F &g){
int n=(*this).size();
int m=g.size();
if(n<m) (*this).resize(m);
for(int i=0;i<m;i++) (*this)[i]+=g[i];
return (*this);
}
F &operator+=(const T &r){
if((*this).size()==0) (*this).resize(1);
(*this)[0]+=r;
return (*this);
}
F &operator-=(const F &g){
int n=(*this).size();
int m=g.size();
if(n<m) (*this).resize(m);
for(int i=0;i<m;i++) (*this)[i]-=g[i];
return (*this);
}
F &operator-=(const T &r){
if((*this).size()==0) (*this).resize(1);
(*this)[0]-=r;
return (*this);
}
F &operator*=(const F &g){
(*this)=convolution((*this),g);
return (*this);
}
F &operator*=(const T &r){
for(int i=0;i<(*this).size();i++) (*this)[i]*=r;
return (*this);
}
F &operator/=(const F &g){
int n=(*this).size();
(*this)=convolution((*this),g.inv());
(*this).resize(n);
return (*this);
}
F &operator/=(const T &r){
r=r.inv();
for(int i=0;i<(*this).size();i++) (*this)[i]*=r;
return (*this);
}
F &operator<<=(const int d) {
int n=(*this).size();
(*this).insert((*this).begin(),d,0);
return *this;
}
F &operator>>=(const int d) {
int n=(*this).size();
(*this).erase((*this).begin(),(*this).begin()+min(n, d));
return *this;
}
F operator*(const T &g) const { return F(*this)*=g;}
F operator-(const T &g) const { return F(*this)-=g;}
F operator+(const T &g) const { return F(*this)+=g;}
F operator/(const T &g) const { return F(*this)/=g;}
F operator*(const F &g) const { return F(*this)*=g;}
F operator-(const F &g) const { return F(*this)-=g;}
F operator+(const F &g) const { return F(*this)+=g;}
F operator/(const F &g) const { return F(*this)/=g;}
F operator%(const F &g) const { return F(*this)%=g;}
F operator<<(const int d) const { return F(*this)<<=d;}
F operator>>(const int d) const { return F(*this)>>=d;}
F pre(int sz) const {
return F(begin(*this), begin(*this) + min((int)this->size(), sz));
}
F inv(int deg=-1) const {
int n=(*this).size();
if(deg==-1) deg=n;
assert(n>0&&(*this)[0]!=T(0));
F g(1);
g[0]=(*this)[0].inv();
while(g.size()<deg){
int m=g.size();
F f(begin(*this),begin(*this)+min(n,2*m));
F r(g);
f.resize(2*m);
r.resize(2*m);
internal::butterfly(f);
internal::butterfly(r);
for(int i=0;i<2*m;i++) f[i]*=r[i];
internal::butterfly_inv(f);
f.erase(f.begin(),f.begin()+m);
f.resize(2*m);
internal::butterfly(f);
for(int i=0;i<2*m;i++) f[i]*=r[i];
internal::butterfly_inv(f);
T in=T(2*m).inv();
in*=-in;
for(int i=0;i<m;i++) f[i]*=in;
g.insert(g.end(),f.begin(),f.begin()+m);
}
return g.pre(deg);
}
T eval(const T &a){
T x=1;
T ret=0;
for(int i=0;i<(*this).size();i++){
ret+=(*this)[i]*x;
x*=a;
}
return ret;
}
void onemul(const int d,const T c){
int n=(*this).size();
for(int i=n-d-1;i>=0;i--){
(*this)[i+d]+=(*this)[i]*c;
}
}
void onediv(const int d,const T c){
int n=(*this).size();
for(int i=0;i<n-d;i++){
(*this)[i+d]-=(*this)[i]*c;
}
}
F diff() const {
int n=(*this).size();
F ret(n);
for(int i=1;i<n;i++) ret[i-1]=(*this)[i]*i;
ret[n-1]=0;
return ret;
}
F integral() const {
int n=(*this).size(),mod =T::mod();
vector<T> inv(n);
inv[1]=1;
for(int i=2;i<n;i++) inv[i]=T(mod)-inv[mod%i]*(mod/i);
F ret(n);
for(int i=n-2;i>=0;i--) ret[i+1]=(*this)[i]*inv[i+1];
ret[0]=0;
return ret;
}
F log(int deg=-1) const {
int n=(*this).size();
if(deg==-1) deg=n;
assert((*this)[0]==T(1));
return ((*this).diff()*(*this).inv(deg)).pre(deg).integral();
}
F exp(int deg=-1) const {
int n=(*this).size();
if(deg==-1) deg=n;
assert(n==0||(*this)[0]==0);
F Inv;
Inv.reserve(deg);
Inv.push_back(T(0));
Inv.push_back(T(1));
auto inplace_integral = [&](F& f) -> void {
const int n = (int)f.size();
int mod=T::mod();
while(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 = [](F &f) -> void {
if(f.empty()) return;
f.erase(begin(f));
T coeff=1,one=1;
for(int i=0;i<f.size();i++){
f[i]*=coeff;
coeff++;
}
};
F b{1,1<(int)(*this).size()?(*this)[1]:0},c{1},z1,z2{1,1};
for(int m=2;m<=deg;m<<=1){
auto y=b;
y.resize(2*m);
internal::butterfly(y);
z1=z2;
F z(m);
for(int i=0;i<m;i++) z[i]=y[i]*z1[i];
internal::butterfly_inv(z);
T si=T(m).inv();
for(int i=0;i<m;i++) z[i]*=si;
fill(begin(z),begin(z)+m/2,T(0));
internal::butterfly(z);
for(int i=0;i<m;i++) z[i]*=-z1[i];
internal::butterfly_inv(z);
for(int i=0;i<m;i++) z[i]*=si;
c.insert(end(c),begin(z)+m/2,end(z));
z2=c;
z2.resize(2*m);
internal::butterfly(z2);
F x(begin((*this)),begin((*this))+min<int>((*this).size(),m));
x.resize(m);
inplace_diff(x);
x.push_back(T(0));
internal::butterfly(x);
for(int i=0;i<m;i++) x[i]*=y[i];
internal::butterfly_inv(x);
for(int i=0;i<m;i++) x[i]*=si;
x-=b.diff();
x.resize(2*m);
for(int i=0;i<m-1;i++) x[m+i]=x[i],x[i]=T(0);
internal::butterfly(x);
for(int i=0;i<2*m;i++) x[i]*=z2[i];
internal::butterfly_inv(x);
T si2=T(m<<1).inv();
for(int i=0;i<2*m;i++) x[i]*=si2;
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));
internal::butterfly(x);
for(int i=0;i<2*m;i++) x[i]*=y[i];
internal::butterfly_inv(x);
for(int i=0;i<2*m;i++) x[i]*=si2;
b.insert(end(b),begin(x)+m,end(x));
}
return b.pre(deg);
}
F pow(ll m){
int n=(*this).size();
int x=0;
while(x<(*this).size()&&(*this)[x]==T(0)){
x++;
}
if(x*m>=n){
F ret(n);
return ret;
}
F f(n-x);
T y=(*this)[x];
for(int i=x;i<n;i++) f[i-x]=(*this)[i]/y;
f=f.log();
for(int i=0;i<f.size();i++) f[i]*=m;
f=f.exp();
y=y.pow(m);
for(int i=0;i<f.size();i++) f[i]*=y;
F ret(n);
for(int i=x*m;i<n;i++) ret[i]=f[i-x*m];
return ret;
}
F shift(T c){
int n=(*this).size();
int mod=T::mod();
vector<T> inv(n+1);
inv[1]=1;
for(int i=2;i<=n;i++) inv[i]=mod-inv[mod%i]*(mod/i);
T x=1;
for(int i=0;i<n;i++){
(*this)[i]*=x;
x*=(i+1);
}
F g(n);
T y=1;
T now=1;
for(int i=0;i<n;i++){
g[n-i-1]=now*y;
now*=c;
y*=inv[i+1];
}
auto tmp=convolution(g,(*this));
T z=1;
for(int i=0;i<n;i++){
(*this)[i]=tmp[n+i-1]*z;
z*=inv[i+1];
}
return (*this);
}
pair<F,F> division(F g){
F f=(*this);
int n=f.size();
int m=g.size();
if(n<m){
F p(0);
return {p,f};
}
F p(n-m+1),q(n-m+1);
for(int i=0;i<n-m+1;i++) p[i]=f[n-i-1];
for(int i=0;i<n-m+1&&i<m;i++) q[i]=g[m-i-1];
p/=q;
for(int i=0;i<(n-m+1)/2;i++) swap(p[i],p[(n-m+1)-i-1]);
g.resize(n);
g*=p;
for(int i=0;i<n;i++) f[i]-=g[i];
int v=n-m+1,u=0;
for(int i=0;i<n;i++) if(f[i].val()) chmax(u,i+1);
p.resize(v);
f.resize(u);
return {p,f};
}
vector<T> multieva(vector<T> p){
int m=p.size();
int n=(*this).size();
int M=1;
int l=0;
while(M<m){
M*=2;
l++;
}
p.resize(M);
swap(m,M);
vector<vector<F>> g(l+1);
g[0].resize(m);
for(int i=0;i<m;i++){
g[0][i].resize(2);
g[0][i][0]=-p[i];
g[0][i][1]=1;
}
for(int i=0;i<l;i++){
g[i+1].resize(m>>(i+1));
for(int j=0;j<(m>>(i+1));j++) g[i+1][j]=g[i][2*j]*g[i][2*j+1];
}
g[l][0]=(*this).division(g[l][0]).se;
for(int i=l;i>=1;i--){
for(int j=0;j<(m>>(i-1));j++){
g[i-1][j]=g[i][j/2].division(g[i-1][j]).se;
}
}
for(int i=0;i<M;i++) if(g[0][i].size()==0) g[0][i].resize(1);
vector<T> ret(M);
for(int i=0;i<M;i++) ret[i]=g[0][i][0];
return ret;
}
};
template<class T>
void GaussJordan(vector<vector<T>> &A,bool is_extended = false){
ll m=A.size(),n=A[0].size();
ll rank=0;
for(int i=0;i<n;i++){
if(is_extended&&i==n-1) break;
ll p=-1;
for(int j=rank;j<m;j++){
if(A[j][i]!=T(0)){
p=j;
break;
}
}
if(p==-1) continue;
swap(A[p],A[rank]);
auto k=A[rank][i];
for(int i2=0;i2<n;i2++){
A[rank][i2]/=k;
}
for(int j=0;j<m;j++){
if(j!=rank&&A[j][i]!=T(0)){
auto fac=A[j][i];
for(int i2=0;i2<n;i2++){
A[j][i2]-=A[rank][i2]*fac;
}
}
}
rank++;
}
}
template<class T>
void linear_equation(vector<vector<T>> a, vector<T> b, vector<T> &res) {
ll m=a.size(),n=a[0].size();
vector<vector<T>> M(m,vector<T>(n+1));
for(int i=0;i<m;i++){
for(int j=0;j<n;j++){
M[i][j]=a[i][j];
}
M[i][n]=b[i];
}
GaussJordan(M,true);
res.assign(n,0);
for(int i=0;i<n;i++) res[i]=M[i][n];
}
template<class F>
pair<F,F> Characteristic_equation(const F &a) {
using T=typename F::value_type;
ll n=a.size();
ll p=n/2;
ll u=p+(p+1);
vector<vector<T>> f(u,vector<T>(u));
f[0][0]=1;
for(int i=1;i<=p;i++){
f[i][i-1]=-1;
}
for(int i=p;i<u;i++){
ll t=0;
for(int j=1+i-p;j<u;j++){
f[j][i]=a[t];
t++;
}
}
vector<T> b(u);
b[0]=1;
vector<T> res(u);
linear_equation(f,b,res);
F X(p),Y(p+1);
for(int i=0;i<p;i++) X[i]=res[i];
for(int j=p;j<res.size();j++) Y[j-p]=res[j];
return {X,Y};
}
template <class T>
T getK(FormalPowerSeries<T> p, FormalPowerSeries<T> q,ll k){
if(p.size()==0) return 0;
if(k==0) return p[0]/q[0];
if(p.size()>=q.size()){
p=p.division(q).se;
}
if(k<0) return T(0);
ll d=q.size();
while(k){
auto qn=q;
for(int i=1;i<d;i+=2) qn[i]*=-1;
p*=qn;
q*=qn;
for(int i=0;i<d-1;i++){
p[i]=p[(i<<1)|(k&1)];
}
for(int i=0;i<d;i++){
q[i]=q[(i<<1)];
}
p.resize(d-1);
q.resize(d);
k/=2;
}
return p[0];
}
using fps=FormalPowerSeries<modint998244353>;
using mint = modint998244353;
struct fpsfraction{
fps mol,den;
fpsfraction(fps a,fps b):mol(a),den(b){}
fpsfraction &operator=(const fpsfraction &g){
(*this).mol=g.mol;
(*this).den=g.den;
return (*this);
}
fpsfraction &operator*=(const fpsfraction &g){
(*this).mol*=g.mol;
(*this).den*=g.den;
return (*this);
}
fpsfraction &operator/=(const fpsfraction &g){
(*this).mol*=g.den;
(*this).den*=g.mol;
return (*this);
}
fpsfraction &operator+=(const fpsfraction &g){
fps f;
f+=g.mol*(*this).den;
f+=g.den*(*this).mol;
(*this).mol=f;
(*this).den*=g.den;
return (*this);
}
fpsfraction &operator-=(const fpsfraction &g){
fps f;
f-=g.mol*(*this).den;
f+=g.den*(*this).mol;
(*this).mol=f;
(*this).den*=g.den;
return (*this);
}
fpsfraction operator*(const fpsfraction &g) const { return fpsfraction(*this)*=g;}
fpsfraction operator-(const fpsfraction &g) const { return fpsfraction(*this)-=g;}
fpsfraction operator+(const fpsfraction &g) const { return fpsfraction(*this)+=g;}
fpsfraction operator/(const fpsfraction &g) const { return fpsfraction(*this)/=g;}
};
int main(){
cincout();
ll n,m;
cin>>n>>m;
fpsfraction f(fps({0}),fps({1}));
for(int i=1;i<=m;i++){
fpsfraction p(fps({1}),fps({1,-1}));
fps g(i+1);
g[0]=1;
g[i]=-1;
fpsfraction q(fps({1}),g);
p-=q;
q=p;
q+=fpsfraction(fps({1}),fps({1}));
p/=q;
f+=p;
}
fpsfraction g(fps({1}),fps({1}));
g-=f;
swap(g.mol,g.den);
cout<<getK(g.mol,g.den,n).val()<<endl;
}