#include using namespace std; using ll = long long; using uint = unsigned int; using ull = unsigned long long; #define rep(i,n) for(int i=0;i=0;i--) #define per1(i,n) for(int i=int(n);i>0;i--) #define all(c) c.begin(),c.end() #define si(x) int(x.size()) #define pb push_back #define eb emplace_back #define fs first #define sc second template using V = vector; template using VV = vector>; template bool chmax(T& x, U y){ if(x bool chmin(T& x, U y){ if(y void mkuni(V& v){sort(all(v));v.erase(unique(all(v)),v.end());} template int lwb(const V& v, const T& a){return lower_bound(all(v),a) - v.begin();} template V Vec(size_t a) { return V(a); } template auto Vec(size_t a, Ts... ts) { return V(ts...))>(a, Vec(ts...)); } template ostream& operator<<(ostream& o,const pair &p){ return o<<"("< ostream& operator<<(ostream& o,const vector &vc){ o<<"{"; for(const T& v:vc) o< D divFloor(D a, D b){ return a / b - (((a ^ b) < 0 && a % b != 0) ? 1 : 0); } template D divCeil(D a, D b) { return a / b + (((a ^ b) > 0 && a % b != 0) ? 1 : 0); } /* 2021/04/14 大幅変更 poly 基本, MultipointEval, Interpolate */ template struct ModInt{ using uint = unsigned int; using ll = long long; using ull = unsigned long long; constexpr static uint mod = mod_; uint v; ModInt():v(0){} ModInt(ll _v):v(normS(_v%mod+mod)){} explicit operator bool() const {return v!=0;} static uint normS(const uint &x){return (x [0 , mod-1] static ModInt make(const uint &x){ModInt m; m.v=x; return m;} ModInt operator+(const ModInt& b) const { return make(normS(v+b.v));} ModInt operator-(const ModInt& b) const { return make(normS(v+mod-b.v));} ModInt operator-() const { return make(normS(mod-v)); } ModInt operator*(const ModInt& b) const { return make((ull)v*b.v%mod);} ModInt operator/(const ModInt& b) const { return *this*b.inv();} ModInt& operator+=(const ModInt& b){ return *this=*this+b;} ModInt& operator-=(const ModInt& b){ return *this=*this-b;} ModInt& operator*=(const ModInt& b){ return *this=*this*b;} ModInt& operator/=(const ModInt& b){ return *this=*this/b;} ModInt& operator++(int){ return *this=*this+1;} ModInt& operator--(int){ return *this=*this-1;} template friend ModInt operator+(T a, const ModInt& b){ return (ModInt(a) += b);} template friend ModInt operator-(T a, const ModInt& b){ return (ModInt(a) -= b);} template friend ModInt operator*(T a, const ModInt& b){ return (ModInt(a) *= b);} template friend ModInt operator/(T a, const ModInt& b){ return (ModInt(a) /= b);} ModInt pow(ll p) const { if(p<0) return inv().pow(-p); ModInt a = 1; ModInt x = *this; while(p){ if(p&1) a *= x; x *= x; p >>= 1; } return a; } ModInt inv() const { // should be prime return pow(mod-2); } // ll extgcd(ll a,ll b,ll &x,ll &y) const{ // ll p[]={a,1,0},q[]={b,0,1}; // while(*q){ // ll t=*p/ *q; // rep(i,3) swap(p[i]-=t*q[i],q[i]); // } // if(p[0]<0) rep(i,3) p[i]=-p[i]; // x=p[1],y=p[2]; // return p[0]; // } // ModInt inv() const { // ll x,y; // extgcd(v,mod,x,y); // return make(normS(x+mod)); // } bool operator==(const ModInt& b) const { return v==b.v;} bool operator!=(const ModInt& b) const { return v!=b.v;} bool operator<(const ModInt& b) const { return v>(istream &o,ModInt& x){ ll tmp; o>>tmp; x=ModInt(tmp); return o; } friend ostream& operator<<(ostream &o,const ModInt& x){ return o<; //using mint = ModInt<1000000007>; V fact,ifact,invs; mint Choose(int a,int b){ if(b<0 || a=0;i--) ifact[i] = ifact[i+1] * (i+1); rep1(i,N-1) invs[i] = fact[i-1] * ifact[i]; } // inplace_fmt (without bit rearranging) // fft: // a[rev(i)] <- \sum_j \zeta^{ij} a[j] // invfft: // a[i] <- (1/n) \sum_j \zeta^{-ij} a[rev(j)] // These two are inversions. // !!! CHANGE IF MOD is unusual !!! const int ORDER_2_MOD_MINUS_1 = 23; // ord_2 (mod-1) const mint PRIMITIVE_ROOT = 3; // primitive root of (Z/pZ)* void fft(V& a){ static constexpr uint mod = mint::mod; static constexpr uint mod2 = mod + mod; static const int H = ORDER_2_MOD_MINUS_1; static const mint root = PRIMITIVE_ROOT; static mint magic[H-1]; int n = si(a); assert(!(n & (n-1))); assert(n >= 1); assert(n <= 1<>(i+2))*3); magic[i] = w; } } int m = n; if(m >>= 1){ rep(i,m){ uint v = a[i+m].v; // < M a[i+m].v = a[i].v + mod - v; // < 2M a[i].v += v; // < 2M } } if(m >>= 1){ mint p = 1; for(int h=0,s=0; s>= 1){ mint p = 1; for(int h=0,s=0; s>= 1){ mint p = 1; for(int h=0,s=0; s= mod2) ? a[i].v - mod2 : a[i].v; // < 2M a[i+m].v = a[i].v + mod - v; // < 3M a[i].v += v; // < 3M } p *= magic[__builtin_ctz(++h)]; } } } rep(i,n){ a[i].v = (a[i].v >= mod2) ? a[i].v - mod2 : a[i].v; // < 2M a[i].v = (a[i].v >= mod) ? a[i].v - mod : a[i].v; // < M } // finally < mod !! } void invfft(V& a){ static constexpr uint mod = mint::mod; static constexpr uint mod2 = mod + mod; static const int H = ORDER_2_MOD_MINUS_1; static const mint root = PRIMITIVE_ROOT; static mint magic[H-1]; int n = si(a); assert(!(n & (n-1))); assert(n >= 1); assert(n <= 1<>(i+2))*3); magic[i] = w.inv(); } } int m = 1; if(m < n>>1){ mint p = 1; for(int h=0,s=0; s>1; m <<= 1){ mint p = 1; for(int h=0,s=0; s>1);i++){ ull x = a[i].v + mod2 - a[i+m].v; // < 4M a[i].v += a[i+m].v; // < 4M a[i].v = (a[i].v >= mod2) ? a[i].v - mod2 : a[i].v; // < 2M a[i+m].v = (p.v * x) % mod; // < M } for(int i=s+(m>>1); i 70ms // verify https://judge.yosupo.jp/submission/44937 V multiply(V a, V b) { int A = si(a), B = si(b); if (!A || !B) return {}; int n = A+B-1; int s = 1; while(s 2 a.resize(s); fft(a); rep(i,s) a[i] *= a[i]; }else{ a.resize(s); fft(a); b.resize(s); fft(b); rep(i,s) a[i] *= b[i]; } invfft(a); a.resize(n); return a; } /* 係数アクセス f[i] でいいが、 配列外参照する可能性があるなら at/set */ template struct Poly: public V{ template Poly(Args...args) : V(args...){} Poly(initializer_list li) : V(li){} int size() const { return V::size(); } mint at(int i) const { return i=size() && !x) return; while(i>=size()) this->pb(0); (*this)[i] = x; return; } mint operator()(mint x) const { // eval mint res = 0; int n = size(); mint a = 1; rep(i,n){ res += a * (*this)[i]; a *= x; } return res; } Poly low(int n) const { // ignore x^n (take first n), but not empty return Poly(this->begin(), this->begin()+min(max(n,1),size())); } Poly rev() const { return Poly(this->rbegin(), this->rend()); } friend ostream& operator<<(ostream &o,const Poly& f){ o << "["; rep(i,f.size()){ o << f[i]; if(i != f.size()-1) o << ","; } o << "]"; return o; } Poly operator-() const { Poly res = *this; for(auto& v: res) v = -v; return res; } template Poly& operator+=(T c){ (*this)[0] += c; return *this; } template Poly& operator-=(T c){ (*this)[0] -= c; return *this; } template Poly& operator*=(T c){ for(auto& v: *this) v *= c; return *this; } template Poly& operator/=(T c){ return *this *= mint(1)/mint(c); } Poly& operator+=(const Poly& r){ if(size() < r.size()) this->resize(r.size(),0); rep(i,r.size()) (*this)[i] += r[i]; return *this; } Poly& operator-=(const Poly& r){ if(size() < r.size()) this->resize(r.size(),0); rep(i,r.size()) (*this)[i] -= r[i]; return *this; } Poly& operator*=(const Poly& r){ return *this = multiply(*this,r); } // 何回も同じrで割り算するなら毎回rinvを計算するのは無駄なので、呼び出し側で一回計算した後直接こっちを呼ぶと良い // 取るべきinvの長さに注意 // 例えば mod r で色々計算したい時は、基本的に deg(r) * 2 長さの多項式を r で割ることになる // とはいえいったん rinv を長く計算したらより短い場合はprefix見るだけだし、 rinv としてムダに長いものを渡しても問題ないので // 割られる多項式として最大の次数を取ればよい Poly quotient(const Poly& r, const Poly& rinv){ int m = r.size(); assert(r[m-1].v); int n = size(); int s = n-m+1; if(s <= 0) return {0}; return (rev().low(s)*rinv.low(s)).low(s).rev(); } Poly& operator/=(const Poly& r){ return *this = quotient(r,r.rev().inv(max(size()-r.size(),0)+1)); } Poly& operator%=(const Poly& r){ *this -= *this/r * r; return *this = low(r.size()-1); } template Poly operator+(T c) const {return Poly(*this) += c; } template Poly operator-(T c) const {return Poly(*this) -= c; } template Poly operator*(T c) const {return Poly(*this) *= c; } template Poly operator/(T c) const {return Poly(*this) /= c; } 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; } Poly operator/(const Poly& r) const {return Poly(*this) /= r; } Poly operator%(const Poly& r) const {return Poly(*this) %= r; } Poly diff() const { Poly g(max(size()-1,0)); rep(i,g.size()) g[i] = (*this)[i+1] * (i+1); return g; } Poly intg() const { assert(si(invs) > size()); Poly g(size()+1); rep(i,size()) g[i+1] = (*this)[i] * invs[i+1]; return g; } Poly square() const { return multiply(*this,*this); } // 1/f(x) mod x^s // N = s = 500000 -> 90ms // inv は 5 回 fft(2n) を呼んでいるので、multiply が 3 回 fft(2n) を呼ぶのと比べると // だいたい multiply の 5/3 倍の時間がかかる // 導出: Newton // fg = 1 mod x^m // (fg-1)^2 = 0 mod x^2m // f(2g-fg^2) = 1 mod x^2m // verify: https://judge.yosupo.jp/submission/44938 Poly inv(int s) const { Poly r(s); r[0] = mint(1)/at(0); for(int n=1;n f = low(2*n); f.resize(2*n); fft(f); V g = r.low(2*n); g.resize(2*n); fft(g); rep(i,2*n) f[i] *= g[i]; invfft(f); rep(i,n) f[i] = 0; fft(f); rep(i,2*n) f[i] *= g[i]; invfft(f); for(int i=n;i f{1},g{1},z{1}; for(int n=1;n d(2*n); rep(i,n) d[n+i] = z[i] - at(i) - at(n+i); fft(d); V g2(2*n); rep(i,n) g2[i] = g[i]; fft(g2); rep(i,n*2) d[i] *= g2[i]; invfft(d); f.resize(n*2); for(int i=n;i= s) break; z = f; fft(z); V eps = g2; rep(i,n*2) eps[i] *= z[i]; invfft(eps); rep(i,n) eps[i] = 0; fft(eps); rep(i,n*2) eps[i] *= g2[i]; invfft(eps); g.resize(n*2); for(int i=n;i T divAt(Poly f, Poly g, ll p){ assert(g.at(0)); while(p){ auto gm = g; for(int i=1;i T linearRecurrenceAt(V a, V c, ll k){ assert(!c.empty() && c[0]); int d = si(c) - 1; assert(si(a) >= d); return divAt((Poly(a.begin(),a.begin()+d) * Poly(c)).low(d), Poly(c), k); } template Poly berlekamp_massey(const vector &u){ int N = u.size(); vector b = {D(-1)}, c = {D(-1)}; D y = D(1); rep1(n,N){ int L = c.size(), M = b.size(); D x = 0; rep(i,L) x += c[i]*u[n-L+i]; b.pb(0); M++; if(!x) continue; D coef = x/y; if(L(c); } mint dp[101][101]; mint nx[101][101]; int main(){ cin.tie(0); ios::sync_with_stdio(false); //DON'T USE scanf/printf/puts !! cout << fixed << setprecision(20); int N,M; cin >> N >> M; auto brute = [&](int N){ V f(N); rep(i,N){ if(i == 0){ for(int m=2;m<=M;m++) nx[m][1] = 1; }else{ mint off = 0; for(int m=2;m<=M;m++) rep(r,m) if(dp[m][r]){ nx[m][(r+1)%m] += dp[m][r]; if(r){ // for(int mm=2;mm<=M;mm++) if(mm != m) nx[mm][1] += dp[m][r]; nx[m][1] -= dp[m][r]; off += dp[m][r]; } } for(int m=2;m<=M;m++) nx[m][1] += off; } rep1(m,100) rep(r,m){ dp[m][r] = nx[m][r], nx[m][r] = 0; if(r) f[i] += dp[m][r]; } } return f; }; V f = brute(M*(M+1)); show(f); auto c = berlekamp_massey(f); reverse(all(c)); show(c); cout << linearRecurrenceAt(f,c,N-1) << endl; show(brute(N)[N-1]); }