polynomial/fps_div_at.hpp

• View this file on GitHub
• Last update: 2022-07-17 14:27:07+09:00
• Include: #include "polynomial/fps_div_at.hpp"

Depends on

• polynomial/polynomial.hpp

Verified with

• polynomial/test/kth_term_of_linearly_recurrent_sequence.test.cpp

Code

#pragma once

#include "polynomial.hpp"

// Bostan-Mori algorithm
template <typename T, typename Mul>
T fps_div_at(Polynomial<T, Mul> f, Polynomial<T, Mul> g, unsigned long long n) {
    const auto even = [](const Polynomial<T, Mul> &f) -> Polynomial<T, Mul> {
        Polynomial<T, Mul> g((f.size() + 1) / 2);
        for (int i = 0; i < g.size(); ++i) {
            g[i] = f[2 * i];
        }
        return g;
    };
    const auto odd = [](const Polynomial<T, Mul> &f) -> Polynomial<T, Mul> {
        Polynomial<T, Mul> g(f.size() / 2);
        for (int i = 0; i < g.size(); ++i) {
            g[i] = f[2 * i + 1];
        }
        return g;
    };
    assert(g.size() > 0 && g[0] != T(0));
    while (n > 0) {
        Polynomial<T, Mul> g_ = g;
        for (int i = 1; i < g_.size(); i += 2) {
            g_[i] = -g_[i];
        }
        g = even(g * g_);
        Polynomial<T, Mul> tmp_f = f * g_;
        if (n % 2 == 1) {
            f = odd(tmp_f);
        } else {
            f = even(tmp_f);
        }
        n /= 2;
    }
    if (f.size() == 0) {
        return T(0);
    } else {
        return f[0] / g[0];
    }
}

#line 2 "polynomial/fps_div_at.hpp"

#line 2 "polynomial/polynomial.hpp"

#include <vector>
#include <utility>
#include <cassert>
#include <algorithm>

template <typename T, typename Mul>
class Polynomial {
    std::vector<T> coeff;
    
public:
    using This = Polynomial<T, Mul>;
    
    Polynomial() : coeff() {}
    Polynomial(int n) : coeff(n, T(0)) {}
    Polynomial(std::vector<T> c) : coeff(std::move(c)) {}
    
    const std::vector<T> &vec() const {
        return coeff;
    }
    
    int size() const {
        return (int) coeff.size();
    }
    
    const T &operator[](int idx) const {
        return coeff[idx];
    }
    T &operator[](int idx) {
        return coeff[idx];
    }
    
    T at(int idx) const {
        if (idx < size()) {
            return coeff[idx];
        } else {
            return T(0);
        }
    }
    
    void pre_(int n) {
        assert(n >= 0);
        coeff.resize(n, T(0));
    }
    This pre(int n) const {
        This tmp(*this);
        tmp.pre_(n);
        return tmp;
    }
    
    T operator()(const T &x) const {
        T p(1), sum(0);
        for (const T &ele : coeff) {
            sum += p * ele;
            p *= x;
        }
        return sum;
    }
    
    This &operator+=(const This &rhs) {
        if (coeff.size() < rhs.coeff.size()) {
            coeff.resize(rhs.coeff.size(), T(0));
        }
        for (int i = 0; i < (int) rhs.coeff.size(); ++i) {
            coeff[i] += rhs.coeff[i];
        }
        return *this;
    }
    friend This operator+(This lhs, const This &rhs) {
        lhs += rhs;
        return lhs;
    }
    This &operator-=(const This &rhs) {
        if (coeff.size() < rhs.coeff.size()) {
            coeff.resize(rhs.coeff.size(), T(0));
        }
        for (int i = 0; i < (int) rhs.coeff.size(); ++i) {
            coeff[i] -= rhs.coeff[i];
        }
        return *this;
    }
    friend This operator-(This lhs, const This &rhs) {
        lhs -= rhs;
        return lhs;
    }
    
    This &operator*=(This rhs) {
        coeff = Mul::mul(std::move(coeff), std::move(rhs.coeff));
        return *this;
    }
    friend This operator*(This lhs, This rhs) {
        return This(Mul::mul(std::move(lhs.coeff), std::move(rhs.coeff)));
    }
    
    This diff() const {
        if (coeff.empty()) {
            return This();
        }
        std::vector<T> c(coeff.size() - 1);
        for (int i = 0; i < (int) c.size(); ++i) {
            c[i] = T(i + 1) * coeff[i + 1];
        }
        return This(c);
    }
    This integ() const {
        std::vector<T> c(coeff.size() + 1, T(0));
        for (int i = 0; i < (int) coeff.size(); ++i) {
            c[i + 1] = coeff[i] / T(i + 1);
        }
        return This(c);
    }
};
#line 4 "polynomial/fps_div_at.hpp"

// Bostan-Mori algorithm
template <typename T, typename Mul>
T fps_div_at(Polynomial<T, Mul> f, Polynomial<T, Mul> g, unsigned long long n) {
    const auto even = [](const Polynomial<T, Mul> &f) -> Polynomial<T, Mul> {
        Polynomial<T, Mul> g((f.size() + 1) / 2);
        for (int i = 0; i < g.size(); ++i) {
            g[i] = f[2 * i];
        }
        return g;
    };
    const auto odd = [](const Polynomial<T, Mul> &f) -> Polynomial<T, Mul> {
        Polynomial<T, Mul> g(f.size() / 2);
        for (int i = 0; i < g.size(); ++i) {
            g[i] = f[2 * i + 1];
        }
        return g;
    };
    assert(g.size() > 0 && g[0] != T(0));
    while (n > 0) {
        Polynomial<T, Mul> g_ = g;
        for (int i = 1; i < g_.size(); i += 2) {
            g_[i] = -g_[i];
        }
        g = even(g * g_);
        Polynomial<T, Mul> tmp_f = f * g_;
        if (n % 2 == 1) {
            f = odd(tmp_f);
        } else {
            f = even(tmp_f);
        }
        n /= 2;
    }
    if (f.size() == 0) {
        return T(0);
    } else {
        return f[0] / g[0];
    }
}

Back to top page