結果

問題 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
権限があれば一括ダウンロードができます

ソースコード

diff #
プレゼンテーションモードにする

#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;
}
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
0