#776888 (C++17) No.2005 Sum of Power Sums

提出ソース
結果

問題	No.2005 Sum of Power Sums
ユーザー	ytqm3
提出日時	2022-07-08 23:46:42
言語	C++17 (gcc 12.3.0 + boost 1.83.0)
結果	TLE
実行時間	-
コード長	17,407 bytes
コンパイル時間	6,535 ms
コンパイル使用メモリ	314,824 KB
実行使用メモリ	125,436 KB
最終ジャッジ日時	2024-06-09 04:24:09
合計ジャッジ時間	16,092 ms
ジャッジサーバーID （参考情報）	judge1 / judge2
このコードへのチャレンジ
（要ログイン）
テストケース

テストケース表示
入力	結果	実行時間実行使用メモリ
testcase_00	AC	170 ms 34,560 KB
testcase_01	AC	172 ms 29,112 KB
testcase_02	AC	197 ms 29,300 KB
testcase_03	AC	168 ms 29,056 KB
testcase_04	AC	166 ms 28,964 KB
testcase_05	AC	173 ms 29,184 KB
testcase_06	AC	178 ms 29,180 KB
testcase_07	AC	174 ms 29,136 KB
testcase_08	AC	186 ms 29,092 KB
testcase_09	AC	197 ms 29,304 KB
testcase_10	AC	182 ms 29,056 KB
testcase_11	AC	205 ms 29,568 KB
testcase_12	AC	171 ms 29,028 KB
testcase_13	AC	208 ms 29,556 KB
testcase_14	TLE	-
testcase_15	TLE	-
testcase_16	--	-
testcase_17	--	-
testcase_18	--	-
testcase_19	--	-
testcase_20	--	-
権限があれば一括ダウンロードができます
ソースコード

raw source code
#include <bits/stdc++.h>
#if __has_include(<atcoder/all>)
#include <atcoder/all>
#endif

namespace ttl{

using namespace std;
using f80=long double;
using i64=int64_t;
using u64=uint64_t;

template<typename T> void Scan_(T& a){
  cin>>a;
}
template<typename T,typename U> void Scan_(pair<T,U>& a){
  Scan_(a.first),Scan_(a.second);
}
template<typename T> void Scan_(vector<T>& a){
  for(auto& v:a){
    Scan_(v);
  }
}
template<typename T> void Scan_(vector<vector<T>>& a){
  for(auto& v:a){
    for(auto& u:v){
      Scan_(u);
    }
  }
}
void Scan(){}
template<typename T,class... Args> void Scan(T& n,Args&... args){
  Scan_(n),Scan(args...);
}
template<typename T> void Print_(T a){
  cout<<a;
}
template<typename T,typename U> void Print_(pair<T,U> a){
  Print_(a.first),cout<<" ";Print_(a.second);
}
void Print_(f80 a){
  printf("%.10Lf",a);
}
template<typename T> void Print(vector<T> a){
  for(size_t i=0;i<a.size();++i){
    Print_(a[i]);
    cout<<" \n"[i==a.size()-1];
  }
}
template<typename T> void Print(vector<vector<T>> a){
  for(auto& v:a){
    for(size_t i=0;i<v.size();++i){
      Print_(v[i]);
      cout<<" \n"[i==v.size()-1];
    }
  }
}
template<typename T> void Print(T a){
  Print_(a);
  cout<<"\n";
}
template<typename T,class... Args> void Print(T a,Args... args){
  Print_(a),cout<<" ",Print(args...);
}

//しおむすびありがとう
template<class Head,class... Tail> struct MultiDimVec{
  using type=vector<typename MultiDimVec<Tail...>::type>;
};
template<class T> struct MultiDimVec<T>{
  using type=T;
};
template<class T,class Arg> vector<T> MakeVec(int n,Arg&& arg){
  return vector<T>(n,arg);
}
template<class T,class... Args> typename MultiDimVec<Args...,T>::type MakeVec(int n,Args&&... args){
  return typename MultiDimVec<Args...,T>::type(n,MakeVec<T>(args...));
}

template<typename T> T Sum(vector<T> a){
  return accumulate(a.begin(),a.end(),T(0));
}

template<typename T> T Rev_(T a){
  reverse(a.begin(),a.end());
  return a;
}
template<typename T> void Rev(T& a){
  reverse(a.begin(),a.end());
}

template<typename T> T Sort_(T a){
  sort(a.begin(),a.end());
  return a;
}
template<typename T> void Sort(T& a){
  sort(a.begin(),a.end());
}
template<typename T> T RSort_(T a){
  sort(a.rbegin(),a.rend());
  return a;
}
template<typename T> void RSort(T& a){
  sort(a.rbegin(),a.rend());
}

template<typename T> T Max(vector<T> a){
  return *max_element(a.begin(),a.end());
}

template<typename T> T Min(vector<T> a){
  return *min_element(a.begin(),a.end());
}

template<typename T> void ChMax(T& a,T b){
  a=max(a,b);
}

template<typename T> void ChMin(T& a,T b){
  a=min(a,b);
}

i64 Tr(i64 n){
  return n*(n+1)/2;
}

i64 PopCnt(u64 k){
  return __builtin_popcountll(k);
}

i64 SqrtF(i64 n){
  i64 ok=0,ng=1e9+5;
  while(std::abs(ok-ng)>1){
    i64 mid=(ok+ng)/2;
    (mid*mid<=n?ok:ng)=mid;
  }
  return ok;
}

i64 FDiv(i64 a,i64 b){
  if(b<0){
    a*=-1,b*=-1;
  }
  if(a<0){
    return -(-a+b-1)/b;
  }
  return a/b;
}

i64 CDiv(i64 a,i64 b){
  if(b<0){
    a*=-1,b*=-1;
  }
  if(a<0){
    return -(-a)/b;
  }
  return (a+b-1)/b;
}

vector<i64> LISSize(vector<i64> A){
  int N=A.size();
  vector<i64> dp(N+1,2e18),res(N+1);
  dp[0]=-1;
  for(int i=0;i<N;++i){
    auto j=(lower_bound(dp.begin(),dp.end(),A[i])-dp.begin())-1;
    dp[j+1]=A[i];
    res[i+1]=max(res[i],i64(j+1));
  }
  return res;
}

#if __has_include(<atcoder/all>)
i64 CountInverse(vector<i64> A){
  int N=A.size();
  auto B=A;
  sort(B.begin(),B.end());
  B.erase(unique(B.begin(),B.end()),B.end());
  map<i64,i64> mp;
  for(size_t i=0;i<B.size();++i){
    mp[B[i]]=i;
  }
  for(int i=0;i<N;++i){
    A[i]=mp[A[i]];
  }
  atcoder::fenwick_tree<i64> fwt(N);
  i64 ans=0;
  for(int i=0;i<N;++i){
    fwt.add(A[i],1);
    ans+=fwt.sum(A[i]+1,N);
  }
  return ans;
}
#endif

struct ESieve{
  int n;
  vector<i64> lpf;
  ESieve(int n_):n(n_),lpf(n_+1,-1){
    for(i64 p=2;p<=n;++p){
      if(lpf[p]!=-1){
        continue;
      }
      for(i64 q=p;q<=n;q+=p){
        if(lpf[q]==-1){
          lpf[q]=p;
        }
      }
    }
  }
  vector<pair<i64,i64>> operator()(int m){
    vector<i64> v;
    while(m!=1){
      v.emplace_back(lpf[m]);
      m/=lpf[m];
    }
    if(v.size()==0){
      return {};
    }
    vector<pair<i64,i64>> res;
    res.emplace_back(v[0],1);
    for(size_t i=1;i<v.size();++i){
      if(v[i-1]!=v[i]){
        res.emplace_back(v[i],1);
      }
      else{
        res.back().second++;
      }
    }
    return res;
  }
};

vector<int> Dist(vector<vector<int>> G,int v){
  int N=G.size();
  vector<int> dst(N,-1);
  queue<int> q;
  dst[v]=0;
  q.emplace(v);
  while(q.size()){
    int t=q.front();
    q.pop();
    for(auto u:G[t]){
      if(dst[u]==-1){
        dst[u]=dst[t]+1;
        q.emplace(u);
      }
    }
  }
  return dst;
}

void SpreadGrid(vector<string>& S,int h,int w,char c){
  int H=S.size(),W=S[0].size();
  auto res=MakeVec<string>(h,"");
  for(int i=0;i<h;++i){
    for(int j=0;j<w;++j){
      if(i<H && j<W){
        res[i]+=S[i][j];
      }
      else{
        res[i]+=c;
      }
    }
  }
  S=res;
}

bool CheckPrime(i64 n){
  if(n<2){
    return 0;
  }
  for(i64 i=2;i*i<=n;++i){
    if(n%i==0){
      return 0;
    }
  }
  return 1;
}

vector<pair<i64,i64>> PrimeFact(i64 n){
  vector<pair<i64,i64>> res;
  for(i64 i=2;i*i<=n;++i){
    if(n%i!=0){
      continue;
    }
    i64 ex=0;
    while(n%i==0){
      ex++,n/=i;
    }
    res.emplace_back(i,ex);
  }
  if(n!=1){
    res.emplace_back(n,1);
  }
  return res;
}

vector<i64> EnumDiv(i64 n){
  vector<i64> res;
  for(i64 i=1;i*i<=n;++i){
    if(n%i!=0){
      continue;
    }
    res.emplace_back(i);
    if(i*i!=n){
      res.emplace_back(n/i);
    }
  }
  sort(res.begin(),res.end());
  return res;
}

template<typename T> vector<pair<T,i64>> RunLenEnc(vector<T> a){
  int n=a.size();
  vector<pair<T,i64>> res;
  T now=a[0];
  int l=1;
  for(int i=1;i<n;++i){
    if(a[i-1]==a[i]){
      l++;
    }
    else{
      res.emplace_back(now,l);
      now=a[i],l=1;
    }
  }
  res.emplace_back(now,l);
  return res;
}

vector<pair<char,i64>> RunLenEnc(string a){
  int n=a.size();
  vector<pair<char,i64>> res;
  char now=a[0];
  int l=1;
  for(int i=1;i<n;++i){
    if(a[i-1]==a[i]){
      l++;
    }
    else{
      res.emplace_back(now,l);
      now=a[i],l=1;
    }
  }
  res.emplace_back(now,l);
  return res;
}

template<typename T> struct Comb{
  vector<T> fac,ifac;
  Comb(int mx=3000000):fac(mx+1,1),ifac(mx+1,1){
    for(int i=1;i<=mx;++i){
      fac[i]=fac[i-1]*i;
    }
    ifac[mx]/=fac[mx];
    for(int i=mx;i>0;--i){
      ifac[i-1]=ifac[i]*i;
    }
  }
  T operator()(int n,int k){
    if(n<0||k<0||n<k){
      return 0;
    }
    return fac[n]*ifac[k]*ifac[n-k];
  }
};

}

using namespace ttl;
template<u64 mod> using ModInt=atcoder::static_modint<mod>;

template<u64 mod> struct FPS{
  using mint=ModInt<mod>;
  vector<mint> val;
  FPS():val({0}){}
  FPS(mint t):val({t}){}
  FPS(int siz):val(max(1,siz)){}
  FPS(initializer_list<mint> init):val(init){}
  FPS(vector<mint> init):val(init){}
  mint &operator[](int i){
    return val[i];
  }
  int size(){
    return val.size();
  }
  void resize(int siz){
    val.resize(siz);
  }
  FPS& resize_(int siz){
    auto tmp=val;
    tmp.resize(siz);
    return (*this)=tmp;
  }
  FPS operator-(){
    for(mint& v:val){
      v=mint(0)-v;
    }
    return (*this);
  }
  FPS& operator+=(mint rhs){
    val[0]+=rhs;
    return (*this);
  }
  FPS& operator-=(mint rhs){
    val[0]-=rhs;
    return (*this);
  }
  FPS& operator*=(mint rhs){
    for(auto& v:val){
      v*=rhs;
    }
    return (*this);
  }
  FPS& operator/=(mint rhs){
    for(auto& v:val){
      v/=rhs;
    }
    return (*this);
  }
  FPS operator+(mint rhs){
    return FPS(*this)+=rhs;
  }
  FPS operator-(mint rhs){
    return FPS(*this)-=rhs;
  }
  FPS operator*(mint rhs){
    return FPS(*this)*=rhs;
  }
  FPS operator/(mint rhs){
    return FPS(*this)/=rhs;
  }
  FPS& operator+=(FPS rhs){
    resize(max(this->size(),rhs.size()));
    for(int i=0;i<int(rhs.size());++i){
      (*this)[i]+=rhs[i];
    }
    return (*this);
  }
  FPS& operator-=(FPS rhs){
    resize(max(this->size(),rhs.size()));
    for(int i=0;i<int(rhs.size());++i){
      (*this)[i]-=rhs[i];
    }
    return (*this);
  }
  FPS& operator*=(FPS rhs){
    val=atcoder::convolution(val,rhs.val);
    return (*this);
  }
  FPS& operator/=(FPS rhs){
    return (*this)*=rhs.inv();
  }
  FPS& operator>>=(int k){
    if(int(val.size())<=k){
      return (*this)={0};
    }
    FPS res=val;
    res.val.erase(res.val.begin(),res.val.begin()+k);
    return (*this)=res;
  }
  FPS& operator<<=(int k){
    FPS res=val;
    res.val.insert(res.val.begin(),k,mint(0));
    return (*this)=res;
  }
  FPS operator+(FPS rhs){
    return FPS(*this)+=rhs;
  }
  FPS operator-(FPS rhs){
    return FPS(*this)-=rhs;
  }
  FPS operator*(FPS rhs){
    return FPS(*this)*=rhs;
  }
  FPS operator/(FPS rhs){
    return FPS(*this)/=rhs;
  }
  FPS operator%(FPS rhs){
    return FPS(*this)%=rhs;
  }
  FPS operator<<(int k){
    return FPS(*this)<<=k;
  }
  FPS operator>>(int k){
    return FPS(*this)>>=k;
  }
  FPS shrink(){
    for(int i=val.size()-1;i>0;--i){
      if(val[i]==0){
        val.pop_back();
      }
      else{
        break;
      }
    }
    return (*this);
  }
  FPS diff_(){
    if(val.size()==1){
      return (*this)={0};
    }
    FPS f(val.size()-1);
    for(size_t i=1;i<val.size();++i){
      f[i-1]=val[i]*i;
    }
    return f;
  }
  FPS integral_(){
    FPS f(val.size()+1);
    for(size_t i=0;i<val.size();++i){
      f[i+1]=val[i]/(i+1);
    }
    return f;
  }
  FPS inv_(int mx=-1){
    if(mx==-1){
      mx=val.size();
    }
    if(val[0]==0){
      assert(0);
    }
    FPS g({mint(1)/val[0]});
    int now=1;
    while(now<mx){
      now<<=1;
      FPS t=(*this);
      t.resize(now);
      t*=g;
      t=-t+mint(2);
      g*=t;
      g.resize(now);
    }
    g.resize(mx);
    return g;
  }
  FPS exp_(int mx=-1){
    if(mx==-1){
      mx=val.size();
    }
    if(val[0]!=0){
      assert(0);
    }
    FPS g(mint(1));
    int now=1;
    while(now<mx){
      now<<=1;
      FPS t=(*this);
      t.resize(now);
      g*=t-g.log_(now)+mint(1);
      g.resize(now);
    }
    g.resize(mx);
    return g;
  }
  FPS log_(int mx=-1){
    if(mx==-1){
      mx=val.size();
    }
    if(val[0]!=1){
      assert(0);
    }
    auto f=(*this);
    f.resize(mx);
    return (f.diff_()/f).integral_().resize_(mx);
  }
  FPS pow_(i64 k,int mx=-1){
    if(mx==-1){
      mx=val.size();
    }
    i64 t=0;
    for(auto v:val){
      if(v==0){
        t++;
      }
      else{
        break;
      }
    }
    auto f=(*this)>>t;
    if(f[0]==0){
      f={0},f.resize(mx);
      return f;
    }
    mint c=f[0];
    f/=c;
    (f.log(mx)*=mint(k)).exp(mx)*=(c.pow(k));
    if(t*k<=mx){
      f<<=t*k;
      f.resize(mx);
      return f;
    }
    else{
      f={0},f.resize(mx);
      return f;
    }
  }
  FPS& diff(){
    return (*this)=diff_();
  }
  FPS& integral(){
    return (*this)=integral_();
  }
  FPS& inv(int mx=-1){
    return (*this)=inv_(mx);
  }
  FPS& exp(int mx=-1){
    return (*this)=exp_(mx);
  }
  FPS& log(int mx=-1){
    return (*this)=log_(mx);
  }
  FPS& pow(i64 k,int mx=-1){
    return (*this)=pow_(k,mx);
  }
  void TaylorShift(mint c){
    Comb<mint> C(val.size());
    vector<mint> g(val.size());
    mint now=1;
    for(size_t i=0;i<val.size();++i){
      val[i]*=C.dat[i];
      g[i]=now*C.idat[i];
      now*=c;
    }
    reverse(val.begin(),val.end());
    g=atcoder::convolution(val,g);
    g.resize(val.size());
    reverse(g.begin(),g.end());
    for(size_t i=0;i<val.size();++i){
      g[i]*=C.idat[i];
    }
    val=g;
  }
};

template<u64 mod> struct Poly{
  typedef ModInt<mod> mint;
  vector<mint> val;
  Poly():val({0}){}
  Poly(mint t):val({t}){}
  Poly(int siz):val(max(1,siz)){}
  Poly(initializer_list<mint> init):val(init){}
  Poly(vector<mint> init):val(init){}
  mint &operator[](int i){
    return val[i];
  }
  int size(){
    return val.size();
  }
  void resize(int siz){
    val.resize(siz);
  }
  void shrink(){
    for(int i=val.size()-1;i>0;--i){
      if(val[i]==0){
        val.pop_back();
      }
      else{
        return;
      }
    }
  }
  Poly operator+(Poly rhs){
    return Poly(*this)+=rhs;
  }
  Poly operator-(Poly rhs){
    return Poly(*this)-=rhs;
  }
  Poly operator*(Poly rhs){
    return Poly(*this)*=rhs;
  }
  Poly operator/(Poly rhs){
    return Poly(*this)/=rhs;
  }
  Poly operator%(Poly rhs){
    return Poly(*this)%=rhs;
  }
  Poly operator-(){
    for(mint& v:val){
      v=mint(0)-v;
    }
    return (*this);
  }
  Poly operator+=(Poly rhs){
    resize(max(this->size(),rhs.size()));
    for(int i=0;i<int(rhs.size());++i){
      (*this)[i]+=rhs[i];
    }
    shrink();
    return (*this);
  }
  Poly operator-=(Poly rhs){
    resize(max(this->size(),rhs.size()));
    for(int i=0;i<int(rhs.size());++i){
      (*this)[i]-=rhs[i];
    }
    shrink();
    return (*this);
  }
  Poly operator*=(Poly rhs){
    val=atcoder::convolution(val,rhs.val);
    return (*this);
  }
  Poly operator/=(Poly rhs){
    if(val.size()<rhs.size()){
      val.resize(0);
      return (*this);
    }
    int rsiz=val.size()-rhs.size()+1;
    reverse(val.begin(),val.end());
    reverse(rhs.val.begin(),rhs.val.end());
    val.resize(rsiz),rhs.inv(rsiz);
    (*this)*=rhs;
    val.resize(rsiz);
    reverse(val.begin(),val.end());
    return (*this);
  }
  Poly operator%=(Poly rhs){
    if(val.size()<rhs.size()){
      return (*this);
    }
    (*this)-=(*this)/rhs*rhs;
    val.resize(rhs.size()-1);
    shrink();
    return (*this);
  }
  mint eval(mint a){
    mint t=1,res=0;
    int n=(*this).size();
    for(int i=0;i<n;++i){
      res+=val[i]*t;
      t*=a;
    }
    return res;
  }
  void diffrent(){
    Poly f(val.size()-1);
    for(int i=1;i<val.size();++i){
      f[i-1]=val[i]*i;
    }
    (*this)=f;
  }
  void integral(){
    Poly f(val.size()+1);
    for(int i=0;i<val.size();++i){
      f[i+1]=val[i]/(i+1);
    }
    (*this)=f;
  }
  void inv(int mx){
    Poly g({mint(1)/val[0]});
    int now=1;
    while(now<mx){
      now<<=1;
      Poly t=(*this);
      t.resize(now);
      t*=g,t=-t;
      t[0]+=2,g*=t;
      g.resize(now);
    }
    g.resize(mx);
    (*this)=g;
  }
  vector<mint> MultiEval(int K){
    int siz=1;
    while(siz<K){
      siz<<=1;
    }
    vector<Poly> t(siz*2-1,{1});
    for(int i=0;i<=K;++i){
      t[i+siz-1]={-i,1};
    }
    for(int i=siz-2;i>=0;--i){
      t[i]=t[2*i+1]*t[2*i+2];
    }
    vector<Poly> g(siz*2-1);
    g[0]=(*this)%t[0];
    for(int i=1;i<siz*2-1;++i){
      g[i]=g[(i-1)/2]%t[i];
    }
    vector<mint> res(K+1);
    for(int i=0;i<=K;++i){
      res[i]=g[i+siz-1][0];
    }
    return res;
  }
  void LagrangeInterpolation(vector<mint> x,vector<mint> y){
    int N=x.size();
    vector<Poly> h(N);
    for(int i=0;i<N;++i){
      h[i]={-x[i],1};
    }
    auto g=Product(h);
    Poly<mod> g_=g;
    g_.diffrent();
    vector<mint> t=g_.MultiEval(x);
    for(int i=0;i<N;++i){
      t[i]=y[i]/t[i];
    }
    vector<Poly<mod>> den(2*N-1),num(2*N-1,{1});
    for(int i=0;i<N;++i){
      den[i+N-1]={-x[i],1};
      num[i+N-1]={t[i]};
    }
    for(int i=N-2;i>=0;--i){
      den[i]=den[2*i+1]*den[2*i+2];
      num[i]=num[2*i+1]*den[2*i+2]+num[2*i+2]*den[2*i+1];
    }
    (*this)=num[0];
  }
};


template<typename T> T Product(vector<T> a){
  int siz=1;
  while(siz<int(a.size())){
    siz<<=1;
  }
  vector<T> res(siz*2-1,{1});
  for(int i=0;i<int(a.size());++i){
    res[i+siz-1]=a[i];
  }
  for(int i=siz-2;i>=0;--i){
    res[i]=res[2*i+1]*res[2*i+2];
  }
  return res[0];
}

template<u64 mod> vector<ModInt<mod>> StirlingNumber2(int N){
  typedef ModInt<mod> mint;
  FPS<mod> f(N+1),g(N+1);
  mint fact=1;
  for(int i=0;i<=N;++i){
    f[i]=(mint(i).pow(N))/fact;
    g[i]=(mint(-1).pow(i))/fact;
    fact*=i+1;
  }
  f*=g;
  vector<mint> res(N+1);
  for(int i=0;i<=N;++i){
    res[i]=f[i];
  }
  return res;
}

template<typename mint> mint lagrange(vector<mint> A,int N,mint T){
  if(T.val()<=N){
    return A[T.val()];
  }
  vector<mint> Q_i(N+1,1);
  for(int i=1;i<=N;++i){
    Q_i[0]*=-i;
  }
  for(int i=1;i<=N;++i){
    Q_i[i]=Q_i[i-1]/(i-N-1)*i;
  }
  vector<mint> c(N+1);
  for(int i=0;i<=N;++i){
    c[i]=A[i]/Q_i[i];
  }
  mint prod=1;
  for(int i=0;i<=N;++i){
    prod*=T-i;
  }
  mint res=0;
  for(int i=0;i<=N;++i){
    res+=c[i]*prod/(T-i);
  }
  return res;
}

int main(){
  constexpr u64 mod=998244353;
  using mint=ModInt<mod>;
  Comb<ModInt<998244353>> Cb;
  //f := sum[j=1,5000] cnt[j]*(M-x)^j
  //h(x) := sum[i=0,x] Cb(N-1+i,i)*f(i) は deg(f)+N 次多項式である
  //すべての k=0,1,...,deg(f)+N について h(k) がわかればよい
  //すべての k=0,1,...,deg(f)+N について g(k)=Cb(N-1+k,k)*f(k) がわかればよい
  //多点評価
  i64 N,M;
  Scan(N,M);
  vector<int> cnt(5001);
  for(int i=0;i<N;++i){
    int k;
    scanf("%d",&k);
    cnt[k]++;
  }
  Poly<mod> f,g={M,-1};
  for(int i=1;i<=5000;++i){
    f+=g*mint(cnt[i]);
    Poly<mod> g_(g.size()+1);
    for(int j=0;j<g.size();++j){
      g_[j]+=g[j]*M;
      g_[j+1]-=g[j];
    }
    g=g_;
  }
  vector<Poly<mod>> J(N);
  J[0]={1};
  for(int i=1;i<=N-1;++i){
    J[i]={1,mint(1)/i};
  }
  auto P=Product(J)*f;
  auto O=P.MultiEval(5000+N);
  for(int i=1;i<=5000+N;++i){
    O[i]+=O[i-1];
  }
  Print(lagrange(O,5000+N,mint(M)).val());
}
yukicoder

結果

テストケース

ソースコード
