FormalPowerSeries - 形式的冪級数
(lib/32-polynomial/FormalPowerSeries.cpp)

• View this file on GitHub
• Last update: 2023-06-05 03:49:56+09:00
• Link: View error logs on GitHub Actions

FormalPowerSeries

• 形式的冪級数、多項式ライブラリ
• 具体的には下記の係数$a=(a_0,…a_{n-1})$を保持するクラス
  • $f(x)=a_0 + a_1x^1 + a_2x^2 + … + a_{n-1}x^{n-1}$
• 以下FormalPowerSeriesの型のことをFpsと称す。
• 0-indexed
  • Fps f(N)で初期化した場合はN-1項までが確保される
• ModIntあるいはRuntimeModIntをテンプレートで渡す必要がある
  • RuntimeModIntのときの使い方はこの提出を参考
• 実装はvectorを継承しているので、[]演算子などが使える。
• ntt-friendlyでないmod(1000’000’007など)でも対応。使い方は同じ。
  • 内部でgarnerを使うかどうかを制御している。
• 実行時modにも対応。使い方は同じ。

コンストラクタ

• 基本的にはvectorに同じ
• FormalPowerSeries(vector< Mint > v)
  • ModInt型のvectorでも初期化できる

FormalPowerSeries メソッド

• int,long long,ModInt型に対して、Fpsとの四則演算をオーバーロードしている。
  • 演算は、fpsの各項に対して行われる。
  • +,-,*: $O(N)$
  • /: $O(N+\log mod)$
• Fps prefix(size_t n)
  • Fpsのprefix
  • $O(n)$
  • Fpsの0項目からn-1項目までのコピーが返却される
• Fps inv(size_t n)
  • Fpsの逆元$1/(f(x))$を返却する
  • $O(N\log N)$
  • nは打ち切りたい項数。指定しなかった場合は今のFpsと同じ大きさで返却される
• Fps diff(void)
  • 微分$( f’(x) )$を返却する
  • $O(N)$
• Fps intg(void)
  • 積分$( \int f(x)dx )$を返却する
  • $O(N)$
  • 初回時だけ、$O(N+\log mod)$
• Fps log(size_t n)
  • 対数$\log f(x)$を返却する
  • $O(NlogN)$
  • nは打ち切りたい項数。指定しなかった場合は今のFpsと同じ大きさで返却される
• Fps exp(size_t n)
  • 指数$\exp f(x)$を返却する
  • $O(N\log N)$
  • nは打ち切りたい項数。指定しなかった場合は今のFpsと同じ大きさで返却される
• Fps pow(long long k,size_t n)
  • 累乗$ f(x)^k$を返却する
  • $O(N\log N)$
  • nは打ち切りたい項数。指定しなかった場合は今のFps.size()*kの大きさで返却される
• Mint eval(Mint x)
  • 一点評価
  • f(x)を返す
  • $O(N)$
• Fps add(const Fps& lhs,const Fps& rhs)
  • 2個のfpsの各要素の和を取ったfpsを返却する。
  • $O(N)$
  • サイズは、lhs,rhsの大きい方に合わせる。
• Fps sub(const Fps& lhs,const Fps& rhs)
  • 2個のfpsの各要素の差を取ったfpsを返却する。
  • lhs[i]-rhs[i]が返ると思って良い
  • $O(N)$
  • サイズは、lhs,rhsの大きい方に合わせる。
• Fps mul(const Fps& lhs,const Fps& rhs)
  • 2個のfpsを畳み込んだfpsを返却する。
  • 内部的にはnttを使用
  • $O(N\log N)$
  • サイズは、lhs.size()+rhs.size()-1になる。
• Fps div(const Fps& lhs,const Fps& rhs)
  • lhsをrhsで割ったfpsを返却する。
  • lhs * rhs.inv() とは違うことに注意。単純な多項式除算である。
  • $O(N\log N)$
  • サイズは、lhs.size()-rhs.size()+1になる。
  • 提出
• Fps mod(const Fps& lhs,const Fps& rhs)
  • lhsをrhsで割ったあまりのfpsを返却する。
  • 単純な多項式除算から発生する余りであることに注意。
  • $O(N\log N)$
  • サイズは、rhs.size()-1になる
• Fps fold_all(vector vfps, size_t n=0)
  • 複数個のfpsをマージテクを使って積を取る
  • nを指定すると、そのprefixで打ち切りながら積を取る
  • 提出
  • 提出
• vector multipoint_evaluation(vector x)
  • 多点評価
  • f(x_0),f(x_1), … , f(x_N) が返る
  • $O(N\log N)$ 定数倍重め。
• inline static Fps interpolation(const vector& x,const vector& y)
  • 多項式補間
  • f(x_0)=y_0,f(x_1)=y_1, … , f(x_N)=y_N を満たすfpsが返る
  • $O(N\log N)$ 定数倍重め。
• Mint nth_term(long long n,const Fps& numerator,const Fps& denominator)
  • f(x)のn項目が返却される
  • $O(N\log N \log n)$
  • nは求めたい項数
  • numeratorは分子のFps (numerator=1かつ求めたい項が大きくないなら、invを使う方が速い)
  • denominatorは分母のFps
  • 提出
  • 提出
• Fps negative_binomial_theorem(const Mint r, const size_t n, const size_t m, const CombinationMod& cm)
  • $(1-rx)^{-n} = Σ {i+n-1}C{n-1} (rx)^{i}$ を計算して返す
  • nは分母にあった時の乗数
  • mは計算後求めたい項数

    SparseFormalPowerSeries メソッド

  • 疎なFps。発展途上なので随時メソッド追加予定。
  • 提出

参考資料

• 【競技プログラミング】形式的冪級数の応用テクニック(前編)
• 形式的冪級数(Formal-Power-Series) / Luzhiled’s memo
• 形式的べき級数解説 / maspyのHP
• (形式的)べき級数と数え上げの写像12相との関係性　前編 - Senの競技プログラミング備忘録
• 今年中に理解する！多項式、母関数、形式的べき級数の競プロでの実践的使い方 - はまやんはまやんはまやん
• 線形漸化式を満たす数列の N 項目を計算するアルゴリズム
• 形式的冪級数（FPS）の inv，log，exp，pow の定数倍の軽いアルゴリズム
• 多項式を使うテクニックたち
• 非再帰で多項式の多点評価をするよ
• 非再帰で補間多項式を求めるよ
• 桁dpもできる
  • [多項式・形式的べき級数]（２）式変形による解法の導出
  • 提出
• $(1-rx)^{-n} = Σ {i+n-1}C{n-1} (rx)^{i}$
  • [AtCoder 参加感想] 2020/06/21:ABC 171
  • 提出

Verified with

• test/polynomial/FormalPowerSeries-exp.test.cpp
• test/polynomial/FormalPowerSeries-interpolation.test.cpp
• test/polynomial/FormalPowerSeries-inv.test.cpp
• test/polynomial/FormalPowerSeries-log.test.cpp
• test/polynomial/FormalPowerSeries-multi-eval.test.cpp
• test/polynomial/FormalPowerSeries-nth.test.cpp
• test/polynomial/FormalPowerSeries-pow.test.cpp

Code

/*
 * @title FormalPowerSeries - 形式的冪級数
 * @docs md/polynomial/FormalPowerSeries.md
 */
template<long long prime, class T = ModInt<prime>> struct FormalPowerSeries;
template<long long prime, class T = ModInt<prime>> class SparseFormalPowerSeries {
    using Mint = T;
    using Sfps = SparseFormalPowerSeries<prime>;
    unordered_map<size_t, T> mp;
    size_t sz;
public:
    friend class FormalPowerSeries<prime>;
    SparseFormalPowerSeries():sz(0){}
    void update(const size_t idx, const Mint val) {
        if(val.x == 0) mp.erase(idx);
        else mp[idx]=val;
        size_t tmp_sz=0;
        for(auto& p:mp) tmp_sz=max(tmp_sz,p.first);
        sz=tmp_sz;
    }
    inline static Sfps mul(const Sfps& lhs, const Sfps& rhs) {
        unordered_map<size_t, T> ret;
        for(auto& [i,a]: lhs.mp) {
            for(auto& [j,b]: rhs.mp) {
                ret[i+j]+=a*b;
            }
        }
        return ret;
    }
    //mapのサイズじゃなくて、最大のindexを返すことに注意
    size_t size() const {return sz+1;}
};
template<long long prime, class T> struct FormalPowerSeries : public vector<T> {
    using vector<T>::vector;
    using Mint  = T;
    using Fps   = FormalPowerSeries<prime>;
    using Sfps   = SparseFormalPowerSeries<prime>;
    inline static constexpr int N_MAX = 1000000;
    Fps even(void) const {Fps ret;for(int i = 0; i < this->size(); i+=2) ret.push_back((*this)[i]);return ret;}
    Fps odd(void)  const {Fps ret;for(int i = 1; i < this->size(); i+=2) ret.push_back((*this)[i]);return ret;}
    Fps symmetry(void) const {Fps ret(this->size());for(int i = 0; i < ret.size(); ++i) ret[i] = (*this)[i]*(i&1?-1:1);return ret;}
public:
    //a0 + a_1*x^1 + a_2*x^2 + ... + a_(n-1)*x^(n-1)
    FormalPowerSeries(const vector<Mint>& v){*this=FormalPowerSeries(v.size());for(int i=0;i<v.size();++i) (*this)[i]=v[i];}
    //TODO constexprにする
    inline static vector<Mint> invs;
    static void invs_build() {
        if(invs.size()) return;
        vector<Mint> fac(N_MAX+1,1),finv(N_MAX+1,1);
        invs.resize(N_MAX+1);
        for(int i=1;i<=N_MAX;i++) fac[i]=fac[i-1]*i;
        finv[N_MAX]=fac[N_MAX].inv();
        for(int i=N_MAX;i>=1;i--) finv[i-1]=finv[i]*i;
        for(int i=1;i<=N_MAX;i++) invs[i]=finv[i]*fac[i-1];
    }
    Fps operator*(const int r) const {return Fps(*this) *= r; }
    Fps &operator*=(const int r) {for(int i=0;i< this->size(); ++i) (*this)[i] *= r; return *this; }
    Fps operator*(const long long int r) const {return Fps(*this) *= r; }
    Fps &operator*=(const long long int r) {for(int i=0;i< this->size(); ++i) (*this)[i] *= r; return *this; }
    Fps operator*(const Mint r) const {return Fps(*this) *= r; }
    Fps &operator*=(const Mint r) {for(int i=0;i< this->size(); ++i) (*this)[i] *= r; return *this; }
    Fps operator/(const int r) const {return Fps(*this) /= r; }
    Fps &operator/=(const int r) {return (*this) *= Mint(r).inv(); }
    Fps operator+(const int r) const {return Fps(*this) += r; }
    Fps &operator+=(const int r) {for(int i=0;i< this->size(); ++i) (*this)[i] += r; return *this; }
    Fps operator-(void) const {return Fps(*this) *= (-1);}
    Fps operator-(const int r) const {return Fps(*this) -= r; }
    Fps &operator-=(const int r) {for(int i=0;i< this->size(); ++i) (*this)[i] -= r; return *this; }
    Fps prefix(size_t n) const {
        return Fps(this->begin(),this->begin()+min(n,this->size()));
    }
    Fps inv(size_t n) const {
        Fps ret({Mint(1)/(*this)[0]});
        for(size_t i=2,m=(n<<1);i < m; i<<=1) {
            Fps h = mul(mul(ret,ret),(this->prefix(i)));
            ret.resize(i);
            for(int j=i>>1;j<i;++j) ret[j] -= h[j];
        }
        return ret.prefix(n);
    }
    Fps inv(void) const {return inv(this->size());}
    Fps diff(void) const {
        Fps ret(max(0,int(this->size())-1));
        for(int i=0;i<ret.size(); ++i) ret[i]=(*this)[i+1]*(i+1);
        return ret;
    }
    Fps intg(size_t n) const {
        invs_build();
        Fps ret(min(this->size()+1,n));
        for(int i=1;i<ret.size(); ++i) ret[i]=(*this)[i-1]*invs[i];
        return ret;
    }
    Fps log(size_t n) const {
        return mul(this->diff(),this->inv(n)).intg(n);
    }
    Fps log(void) const {return log(this->size());}
    Fps exp(size_t n) const {
        Fps ret(1,1);
        for(size_t i=2,m=(n<<1);i<m;i<<=1) {
            Fps h = mul(ret,(sub(this->prefix(i),ret.log(i))));
            ret.resize(i);
            for(int j=i>>1;j<i;++j) ret[j] += h[j];
        }
        return ret.prefix(n);
    }
    Fps exp(void) const {return exp(this->size());}
    Fps pow(long long k,size_t n) const {
        Fps ret(n,0);
        for(size_t i=0,m = min(n,this->size()); i < m && i*k < n; ++i) {
            if((*this)[i].x == 0) continue;
            Mint t0=(*this)[i], t0_inv=t0.inv();
            Fps tmp(n-i);for(int j=i;j<m; ++j) tmp[j-i]=(*this)[j]*t0_inv;
            tmp = (tmp.log(n)*k).exp(n)*(t0.pow(k));
            for(int j=0;j+i*k<n;++j) ret[j+i*k] = tmp[j];
            break;
        }
        return ret;
    }
    Fps pow(long long k) const {return pow(k,this->size());}
    Mint eval(Mint x) const {
        Mint base = 1,ret = 0;
        for(size_t i=0;i<this->size();++i) {
            ret += (*this)[i]*base;
            base *= x;
        }
        return ret;
    }
    inline static Fps add(const Fps& lhs,const Fps& rhs) {
        size_t n = lhs.size(), m = rhs.size();
        Fps res(max(n,m),0);
        for(int i=0;i<n;++i) res[i] += lhs[i];
        for(int i=0;i<m;++i) res[i] += rhs[i];
        return res;
    }
    inline static Fps sub(const Fps& lhs,const Fps& rhs) {
        size_t n = lhs.size(), m = rhs.size();
        Fps res(max(n,m),0);
        for(int i=0;i<n;++i) res[i] += lhs[i];
        for(int i=0;i<m;++i) res[i] -= rhs[i];
        return res;
    }
    inline static Fps mul(const Fps& lhs, const Fps& rhs) {
        return NumberTheoreticalTransform<prime>::convolution(lhs,rhs);
    }
    inline static Fps mul(const Fps& lhs, const Sfps& rhs) {
        size_t n = lhs.size() + rhs.size();
        Fps res(n,0);
        for(auto& [i,a]:rhs.mp) {
            for (int j = 0; j < lhs.size(); ++j) res[i+j]+=a*lhs[j];
        }
        return res;
    }
    inline static Fps div(Fps lhs, Fps rhs) {
        while(lhs.size() && lhs.back().x == 0) lhs.pop_back();
        while(rhs.size() && rhs.back().x == 0) rhs.pop_back();
        int n = lhs.size(), m = rhs.size();
        if(n < m) return Fps(1,0);
        reverse(lhs.begin(),lhs.end());
        reverse(rhs.begin(),rhs.end());
        auto f = mul(lhs,rhs.inv(n-m+1)).prefix(n-m+1);
        reverse(f.begin(),f.end());
        return f;
    }
    inline static Fps mod(const Fps& lhs, const Fps& rhs) {
        int m = rhs.size();
        auto f = sub(lhs,mul(div(lhs,rhs).prefix(m),rhs)).prefix(m);
        while(f.size() && f.back().x==0) f.pop_back();
        return f;
    }
    inline static Fps fold_all(vector<Fps> vfps, size_t n=0) {
        if(vfps.empty()) return {};
        priority_queue<pair<size_t,size_t>, vector<pair<size_t,size_t>>, greater<>> pq;
        for(size_t i=0;i<vfps.size(); ++i) pq.emplace(vfps[i].size(), i);
        while(pq.size()>1) {
            auto l=pq.top().second; pq.pop();
            auto r=pq.top().second; pq.pop();
            vfps[l]=mul(vfps[l],vfps[r]);
            if(n) vfps[l].resize(n);
            vfps[r]={};
            pq.emplace(vfps[l].size(), l);
        }
        auto ret=pq.top().second; pq.pop();
        return vfps[ret];
    }
    vector<Mint> multipoint_evaluation(vector<Mint> x) {
        int n = x.size(),m;
        for(m=1;m<n;m<<=1);
        vector<Fps> f(2*m,Fps(1,1));
        for(int i=0;i<n;++i) f[i+m] = Fps({-x[i],1});
        for(int i=m-1;i;--i) f[i] = mul(f[(i<<1)|0],f[(i<<1)|1]);
        f[1] = mod(*this,f[1]);
        for(int i=2;i<m+n;++i) f[i] = mod(f[i>>1],f[i]);
        for(int i=0;i<n;++i)   x[i] = f[i+m][0];
        return x;
    }
    inline static Fps interpolation(const vector<Mint>& x,const vector<Mint>& y) {
        int n = x.size(),m;
        for(m=1;m<n;m<<=1);
        vector<Fps> f(2*m,Fps(1,1)),g(2*m);
        for(int i=0;i<n;++i) f[i+m] = Fps({-x[i],1});
        for(int i=m-1;i;--i) f[i] = mul(f[(i<<1)|0],f[(i<<1)|1]);
        g[1] = mod(f[1].diff(), f[1]);
        for(int i=2;i<m+n;++i) g[i] = mod(g[i>>1],f[i]);
        for(int i=0;i<n;++i) g[i+m] = Fps(1, y[i] / g[i+m][0]);
        for(int i=m-1;i;--i) g[i] = add(mul(g[(i<<1)|0],f[(i<<1)|1]),mul(f[(i<<1)|0],g[(i<<1)|1]));
        return g[1];
    }
    inline static Mint nth_term(long long n, Fps numerator,Fps denominator) {
        while(n) {
            numerator    = mul(numerator,denominator.symmetry());
            numerator    = ((n&1)?numerator.odd():numerator.even());
            denominator  = mul(denominator,denominator.symmetry());
            denominator  = denominator.even();
            n >>= 1;
        }
        return numerator[0];
    }
    //(1-rx)^(-n) = Σ (i+n-1)_C_(n-1) (rx)^i をm項まで計算して返してくれる
    inline static Fps negative_binomial_theorem(const Mint r, const size_t n, const size_t m, const CombinationMod<prime>& cm) {
        assert(n-1+m-1 <= cm.size());
        Fps res(m,0);
        for(int i=0;i<m;++i) {
            res[i]= Mint(cm.binom(i+n-1,n-1));
        }
        return res;
    }
    friend ostream &operator<<(ostream &os, const Fps& fps) {os << "{" << fps[0];for(int i=1;i<fps.size();++i) os << ", " << fps[i];return os << "}";}
};

#line 1 "lib/32-polynomial/FormalPowerSeries.cpp"
/*
 * @title FormalPowerSeries - 形式的冪級数
 * @docs md/polynomial/FormalPowerSeries.md
 */
template<long long prime, class T = ModInt<prime>> struct FormalPowerSeries;
template<long long prime, class T = ModInt<prime>> class SparseFormalPowerSeries {
    using Mint = T;
    using Sfps = SparseFormalPowerSeries<prime>;
    unordered_map<size_t, T> mp;
    size_t sz;
public:
    friend class FormalPowerSeries<prime>;
    SparseFormalPowerSeries():sz(0){}
    void update(const size_t idx, const Mint val) {
        if(val.x == 0) mp.erase(idx);
        else mp[idx]=val;
        size_t tmp_sz=0;
        for(auto& p:mp) tmp_sz=max(tmp_sz,p.first);
        sz=tmp_sz;
    }
    inline static Sfps mul(const Sfps& lhs, const Sfps& rhs) {
        unordered_map<size_t, T> ret;
        for(auto& [i,a]: lhs.mp) {
            for(auto& [j,b]: rhs.mp) {
                ret[i+j]+=a*b;
            }
        }
        return ret;
    }
    //mapのサイズじゃなくて、最大のindexを返すことに注意
    size_t size() const {return sz+1;}
};
template<long long prime, class T> struct FormalPowerSeries : public vector<T> {
    using vector<T>::vector;
    using Mint  = T;
    using Fps   = FormalPowerSeries<prime>;
    using Sfps   = SparseFormalPowerSeries<prime>;
    inline static constexpr int N_MAX = 1000000;
    Fps even(void) const {Fps ret;for(int i = 0; i < this->size(); i+=2) ret.push_back((*this)[i]);return ret;}
    Fps odd(void)  const {Fps ret;for(int i = 1; i < this->size(); i+=2) ret.push_back((*this)[i]);return ret;}
    Fps symmetry(void) const {Fps ret(this->size());for(int i = 0; i < ret.size(); ++i) ret[i] = (*this)[i]*(i&1?-1:1);return ret;}
public:
    //a0 + a_1*x^1 + a_2*x^2 + ... + a_(n-1)*x^(n-1)
    FormalPowerSeries(const vector<Mint>& v){*this=FormalPowerSeries(v.size());for(int i=0;i<v.size();++i) (*this)[i]=v[i];}
    //TODO constexprにする
    inline static vector<Mint> invs;
    static void invs_build() {
        if(invs.size()) return;
        vector<Mint> fac(N_MAX+1,1),finv(N_MAX+1,1);
        invs.resize(N_MAX+1);
        for(int i=1;i<=N_MAX;i++) fac[i]=fac[i-1]*i;
        finv[N_MAX]=fac[N_MAX].inv();
        for(int i=N_MAX;i>=1;i--) finv[i-1]=finv[i]*i;
        for(int i=1;i<=N_MAX;i++) invs[i]=finv[i]*fac[i-1];
    }
    Fps operator*(const int r) const {return Fps(*this) *= r; }
    Fps &operator*=(const int r) {for(int i=0;i< this->size(); ++i) (*this)[i] *= r; return *this; }
    Fps operator*(const long long int r) const {return Fps(*this) *= r; }
    Fps &operator*=(const long long int r) {for(int i=0;i< this->size(); ++i) (*this)[i] *= r; return *this; }
    Fps operator*(const Mint r) const {return Fps(*this) *= r; }
    Fps &operator*=(const Mint r) {for(int i=0;i< this->size(); ++i) (*this)[i] *= r; return *this; }
    Fps operator/(const int r) const {return Fps(*this) /= r; }
    Fps &operator/=(const int r) {return (*this) *= Mint(r).inv(); }
    Fps operator+(const int r) const {return Fps(*this) += r; }
    Fps &operator+=(const int r) {for(int i=0;i< this->size(); ++i) (*this)[i] += r; return *this; }
    Fps operator-(void) const {return Fps(*this) *= (-1);}
    Fps operator-(const int r) const {return Fps(*this) -= r; }
    Fps &operator-=(const int r) {for(int i=0;i< this->size(); ++i) (*this)[i] -= r; return *this; }
    Fps prefix(size_t n) const {
        return Fps(this->begin(),this->begin()+min(n,this->size()));
    }
    Fps inv(size_t n) const {
        Fps ret({Mint(1)/(*this)[0]});
        for(size_t i=2,m=(n<<1);i < m; i<<=1) {
            Fps h = mul(mul(ret,ret),(this->prefix(i)));
            ret.resize(i);
            for(int j=i>>1;j<i;++j) ret[j] -= h[j];
        }
        return ret.prefix(n);
    }
    Fps inv(void) const {return inv(this->size());}
    Fps diff(void) const {
        Fps ret(max(0,int(this->size())-1));
        for(int i=0;i<ret.size(); ++i) ret[i]=(*this)[i+1]*(i+1);
        return ret;
    }
    Fps intg(size_t n) const {
        invs_build();
        Fps ret(min(this->size()+1,n));
        for(int i=1;i<ret.size(); ++i) ret[i]=(*this)[i-1]*invs[i];
        return ret;
    }
    Fps log(size_t n) const {
        return mul(this->diff(),this->inv(n)).intg(n);
    }
    Fps log(void) const {return log(this->size());}
    Fps exp(size_t n) const {
        Fps ret(1,1);
        for(size_t i=2,m=(n<<1);i<m;i<<=1) {
            Fps h = mul(ret,(sub(this->prefix(i),ret.log(i))));
            ret.resize(i);
            for(int j=i>>1;j<i;++j) ret[j] += h[j];
        }
        return ret.prefix(n);
    }
    Fps exp(void) const {return exp(this->size());}
    Fps pow(long long k,size_t n) const {
        Fps ret(n,0);
        for(size_t i=0,m = min(n,this->size()); i < m && i*k < n; ++i) {
            if((*this)[i].x == 0) continue;
            Mint t0=(*this)[i], t0_inv=t0.inv();
            Fps tmp(n-i);for(int j=i;j<m; ++j) tmp[j-i]=(*this)[j]*t0_inv;
            tmp = (tmp.log(n)*k).exp(n)*(t0.pow(k));
            for(int j=0;j+i*k<n;++j) ret[j+i*k] = tmp[j];
            break;
        }
        return ret;
    }
    Fps pow(long long k) const {return pow(k,this->size());}
    Mint eval(Mint x) const {
        Mint base = 1,ret = 0;
        for(size_t i=0;i<this->size();++i) {
            ret += (*this)[i]*base;
            base *= x;
        }
        return ret;
    }
    inline static Fps add(const Fps& lhs,const Fps& rhs) {
        size_t n = lhs.size(), m = rhs.size();
        Fps res(max(n,m),0);
        for(int i=0;i<n;++i) res[i] += lhs[i];
        for(int i=0;i<m;++i) res[i] += rhs[i];
        return res;
    }
    inline static Fps sub(const Fps& lhs,const Fps& rhs) {
        size_t n = lhs.size(), m = rhs.size();
        Fps res(max(n,m),0);
        for(int i=0;i<n;++i) res[i] += lhs[i];
        for(int i=0;i<m;++i) res[i] -= rhs[i];
        return res;
    }
    inline static Fps mul(const Fps& lhs, const Fps& rhs) {
        return NumberTheoreticalTransform<prime>::convolution(lhs,rhs);
    }
    inline static Fps mul(const Fps& lhs, const Sfps& rhs) {
        size_t n = lhs.size() + rhs.size();
        Fps res(n,0);
        for(auto& [i,a]:rhs.mp) {
            for (int j = 0; j < lhs.size(); ++j) res[i+j]+=a*lhs[j];
        }
        return res;
    }
    inline static Fps div(Fps lhs, Fps rhs) {
        while(lhs.size() && lhs.back().x == 0) lhs.pop_back();
        while(rhs.size() && rhs.back().x == 0) rhs.pop_back();
        int n = lhs.size(), m = rhs.size();
        if(n < m) return Fps(1,0);
        reverse(lhs.begin(),lhs.end());
        reverse(rhs.begin(),rhs.end());
        auto f = mul(lhs,rhs.inv(n-m+1)).prefix(n-m+1);
        reverse(f.begin(),f.end());
        return f;
    }
    inline static Fps mod(const Fps& lhs, const Fps& rhs) {
        int m = rhs.size();
        auto f = sub(lhs,mul(div(lhs,rhs).prefix(m),rhs)).prefix(m);
        while(f.size() && f.back().x==0) f.pop_back();
        return f;
    }
    inline static Fps fold_all(vector<Fps> vfps, size_t n=0) {
        if(vfps.empty()) return {};
        priority_queue<pair<size_t,size_t>, vector<pair<size_t,size_t>>, greater<>> pq;
        for(size_t i=0;i<vfps.size(); ++i) pq.emplace(vfps[i].size(), i);
        while(pq.size()>1) {
            auto l=pq.top().second; pq.pop();
            auto r=pq.top().second; pq.pop();
            vfps[l]=mul(vfps[l],vfps[r]);
            if(n) vfps[l].resize(n);
            vfps[r]={};
            pq.emplace(vfps[l].size(), l);
        }
        auto ret=pq.top().second; pq.pop();
        return vfps[ret];
    }
    vector<Mint> multipoint_evaluation(vector<Mint> x) {
        int n = x.size(),m;
        for(m=1;m<n;m<<=1);
        vector<Fps> f(2*m,Fps(1,1));
        for(int i=0;i<n;++i) f[i+m] = Fps({-x[i],1});
        for(int i=m-1;i;--i) f[i] = mul(f[(i<<1)|0],f[(i<<1)|1]);
        f[1] = mod(*this,f[1]);
        for(int i=2;i<m+n;++i) f[i] = mod(f[i>>1],f[i]);
        for(int i=0;i<n;++i)   x[i] = f[i+m][0];
        return x;
    }
    inline static Fps interpolation(const vector<Mint>& x,const vector<Mint>& y) {
        int n = x.size(),m;
        for(m=1;m<n;m<<=1);
        vector<Fps> f(2*m,Fps(1,1)),g(2*m);
        for(int i=0;i<n;++i) f[i+m] = Fps({-x[i],1});
        for(int i=m-1;i;--i) f[i] = mul(f[(i<<1)|0],f[(i<<1)|1]);
        g[1] = mod(f[1].diff(), f[1]);
        for(int i=2;i<m+n;++i) g[i] = mod(g[i>>1],f[i]);
        for(int i=0;i<n;++i) g[i+m] = Fps(1, y[i] / g[i+m][0]);
        for(int i=m-1;i;--i) g[i] = add(mul(g[(i<<1)|0],f[(i<<1)|1]),mul(f[(i<<1)|0],g[(i<<1)|1]));
        return g[1];
    }
    inline static Mint nth_term(long long n, Fps numerator,Fps denominator) {
        while(n) {
            numerator    = mul(numerator,denominator.symmetry());
            numerator    = ((n&1)?numerator.odd():numerator.even());
            denominator  = mul(denominator,denominator.symmetry());
            denominator  = denominator.even();
            n >>= 1;
        }
        return numerator[0];
    }
    //(1-rx)^(-n) = Σ (i+n-1)_C_(n-1) (rx)^i をm項まで計算して返してくれる
    inline static Fps negative_binomial_theorem(const Mint r, const size_t n, const size_t m, const CombinationMod<prime>& cm) {
        assert(n-1+m-1 <= cm.size());
        Fps res(m,0);
        for(int i=0;i<m;++i) {
            res[i]= Mint(cm.binom(i+n-1,n-1));
        }
        return res;
    }
    friend ostream &operator<<(ostream &os, const Fps& fps) {os << "{" << fps[0];for(int i=1;i<fps.size();++i) os << ", " << fps[i];return os << "}";}
};

Back to top page