#pragma once
#include "rational-binomial.hpp"
#include "rational.hpp"
template <typename R = Rational>
struct FormalPowerSeries_rational : vector<R> {
using vector<R>::vector;
using fps = FormalPowerSeries_rational;
fps &operator+=(const fps &r) {
if (r.size() > this->size()) this->resize(r.size());
for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i];
return *this;
}
fps &operator+=(const R &r) {
if (this->empty()) this->resize(1);
(*this)[0] += r;
return *this;
}
fps &operator-=(const fps &r) {
if (r.size() > this->size()) this->resize(r.size());
for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i];
return *this;
}
fps &operator-=(const R &r) {
if (this->empty()) this->resize(1);
(*this)[0] -= r;
return *this;
}
fps &operator*=(const fps &r) {
int n = this->size() + r.size() - 1;
fps f(n);
for (int i = 0; i < (int)this->size(); i++) {
for (int j = 0; j < (int)r.size(); j++) {
f[i + j] += (*this)[i] * r[j];
}
}
return *this = f;
}
fps &operator*=(const R &v) {
for (int k = 0; k < (int)this->size(); k++) (*this)[k] *= v;
return *this;
}
fps &operator/=(const fps &r) {
if (this->size() < r.size()) {
this->clear();
return *this;
}
int n = this->size() - r.size() + 1;
fps f(*this), g(r);
g.shrink();
R coeff = g.back().inverse();
for (auto &x : g) x *= coeff;
int deg = (int)f.size() - (int)g.size() + 1;
int gs = g.size();
fps quo(deg);
for (int i = deg - 1; i >= 0; i--) {
quo[i] = f[i + gs - 1];
for (int j = 0; j < gs; j++) f[i + j] -= quo[i] * g[j];
}
*this = quo * coeff;
this->resize(n, R(0));
return *this;
}
fps &operator%=(const fps &r) {
*this -= *this / r * r;
shrink();
return *this;
}
fps operator+(const fps &r) const { return fps(*this) += r; }
fps operator+(const R &v) const { return fps(*this) += v; }
fps operator-(const fps &r) const { return fps(*this) -= r; }
fps operator-(const R &v) const { return fps(*this) -= v; }
fps operator*(const fps &r) const { return fps(*this) *= r; }
fps operator*(const R &v) const { return fps(*this) *= v; }
fps operator/(const fps &r) const { return fps(*this) /= r; }
fps operator%(const fps &r) const { return fps(*this) %= r; }
fps operator-() const {
fps ret(this->size());
for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i];
return ret;
}
void shrink() {
while (this->size() && this->back() == R(0)) this->pop_back();
}
fps rev() const {
fps ret(*this);
reverse(begin(ret), end(ret));
return ret;
}
fps dot(fps r) const {
fps ret(min(this->size(), r.size()));
for (int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i];
return ret;
}
// 前 sz 項を取ってくる。sz に足りない項は 0 埋めする
fps pre(int sz) const {
fps ret(begin(*this), begin(*this) + min((int)this->size(), sz));
if ((int)ret.size() < sz) ret.resize(sz);
return ret;
}
fps operator>>(int sz) const {
if ((int)this->size() <= sz) return {};
fps ret(*this);
ret.erase(ret.begin(), ret.begin() + sz);
return ret;
}
fps operator<<(int sz) const {
fps ret(*this);
ret.insert(ret.begin(), sz, R(0));
return ret;
}
fps diff() const {
const int n = (int)this->size();
fps ret(max(0, n - 1));
R one(1), coeff(1);
for (int i = 1; i < n; i++) {
ret[i - 1] = (*this)[i] * coeff;
coeff += one;
}
return ret;
}
fps integral() const {
const int n = (int)this->size();
fps ret(n + 1);
for (int i = 0; i < n; i++) ret[i + 1] = (*this)[i] / (i + 1);
return ret;
}
R eval(R x) const {
R r = 0, w = 1;
for (auto &v : *this) r += w * v, w *= x;
return r;
}
fps inv(int deg = -1) const {
assert((*this)[0] != R(0));
if (deg == -1) deg = (*this).size();
fps ret{R(1) / (*this)[0]};
for (int i = 1; i < deg; i <<= 1) {
ret = (ret + ret - ret * ret * (*this).pre(i << 1)).pre(i << 1);
}
return ret.pre(deg);
}
fps log(int deg = -1) const {
assert(!(*this).empty() && (*this)[0] == R(1));
if (deg == -1) deg = (int)this->size();
return (this->diff() * this->inv(deg)).pre(deg - 1).integral();
}
fps exp(int deg = -1) const {
assert((*this).size() == 0 || (*this)[0] == R(0));
if (deg == -1) deg = (int)this->size();
fps ret{R(1)};
for (int i = 1; i < deg; i <<= 1) {
ret = (ret * (pre(i << 1) + R(1) - ret.log(i << 1))).pre(i << 1);
}
return ret.pre(deg);
}
fps pow(int64_t k, int deg = -1) const {
const int n = (int)this->size();
if (deg == -1) deg = n;
if (k == 0) {
fps ret(deg);
if (deg) ret[0] = 1;
return ret;
}
for (int i = 0; i < n; i++) {
if ((*this)[i] != R(0)) {
R rev = R(1) / (*this)[i];
fps ret = (((*this * rev) >> i).log(deg) * k).exp(deg);
ret *= (*this)[i].pow(k);
ret = (ret << (i * k)).pre(deg);
if ((int)ret.size() < deg) ret.resize(deg, R(0));
return ret;
}
if (__int128_t(i + 1) * k >= deg) return fps(deg, R(0));
}
return fps(deg, R(0));
}
};
#line 2 "math/rational-fps.hpp"
#line 2 "math/rational-binomial.hpp"
#line 2 "math/rational.hpp"
#include <cassert>
#include <numeric>
#include <vector>
using namespace std;
#line 2 "internal/internal-type-traits.hpp"
#include <type_traits>
using namespace std;
namespace internal {
template <typename T>
using is_broadly_integral =
typename conditional_t<is_integral_v<T> || is_same_v<T, __int128_t> ||
is_same_v<T, __uint128_t>,
true_type, false_type>::type;
template <typename T>
using is_broadly_signed =
typename conditional_t<is_signed_v<T> || is_same_v<T, __int128_t>,
true_type, false_type>::type;
template <typename T>
using is_broadly_unsigned =
typename conditional_t<is_unsigned_v<T> || is_same_v<T, __uint128_t>,
true_type, false_type>::type;
#define ENABLE_VALUE(x) \
template <typename T> \
constexpr bool x##_v = x<T>::value;
ENABLE_VALUE(is_broadly_integral);
ENABLE_VALUE(is_broadly_signed);
ENABLE_VALUE(is_broadly_unsigned);
#undef ENABLE_VALUE
#define ENABLE_HAS_TYPE(var) \
template <class, class = void> \
struct has_##var : false_type {}; \
template <class T> \
struct has_##var<T, void_t<typename T::var>> : true_type {}; \
template <class T> \
constexpr auto has_##var##_v = has_##var<T>::value;
#define ENABLE_HAS_VAR(var) \
template <class, class = void> \
struct has_##var : false_type {}; \
template <class T> \
struct has_##var<T, void_t<decltype(T::var)>> : true_type {}; \
template <class T> \
constexpr auto has_##var##_v = has_##var<T>::value;
} // namespace internal
#line 2 "math-fast/gcd.hpp"
#include <algorithm>
using namespace std;
namespace BinaryGCDImpl {
using u64 = unsigned long long;
using i8 = char;
u64 binary_gcd(u64 a, u64 b) {
if (a == 0 || b == 0) return a + b;
i8 n = __builtin_ctzll(a);
i8 m = __builtin_ctzll(b);
a >>= n;
b >>= m;
n = min(n, m);
while (a != b) {
u64 d = a - b;
i8 s = __builtin_ctzll(d);
bool f = a > b;
b = f ? b : a;
a = (f ? d : -d) >> s;
}
return a << n;
}
using u128 = __uint128_t;
// a > 0
int ctz128(u128 a) {
u64 lo = a & u64(-1);
return lo ? __builtin_ctzll(lo) : 64 + __builtin_ctzll(a >> 64);
}
u128 binary_gcd128(u128 a, u128 b) {
if (a == 0 || b == 0) return a + b;
i8 n = ctz128(a);
i8 m = ctz128(b);
a >>= n;
b >>= m;
n = min(n, m);
while (a != b) {
u128 d = a - b;
i8 s = ctz128(d);
bool f = a > b;
b = f ? b : a;
a = (f ? d : -d) >> s;
}
return a << n;
}
} // namespace BinaryGCDImpl
long long binary_gcd(long long a, long long b) {
return BinaryGCDImpl::binary_gcd(abs(a), abs(b));
}
__int128_t binary_gcd128(__int128_t a, __int128_t b) {
if (a < 0) a = -a;
if (b < 0) b = -b;
return BinaryGCDImpl::binary_gcd128(a, b);
}
/**
* @brief binary GCD
*/
#line 10 "math/rational.hpp"
// T : 値, U : 比較用
template <typename T, typename U>
struct RationalBase {
using R = RationalBase;
using Key = T;
T x, y;
RationalBase() : x(0), y(1) {}
template <typename T1>
RationalBase(const T1& _x) : RationalBase<T, U>(_x, T1{1}) {}
template <typename T1, typename T2>
RationalBase(const pair<T1, T2>& _p)
: RationalBase<T, U>(_p.first, _p.second) {}
template <typename T1, typename T2>
RationalBase(const T1& _x, const T2& _y) : x(_x), y(_y) {
assert(y != 0);
if (y == -1) x = -x, y = -y;
if (y != 1) {
T g;
if constexpr (internal::is_broadly_integral_v<T>) {
if constexpr (sizeof(T) == 16) {
g = binary_gcd128(x, y);
} else {
g = binary_gcd(x, y);
}
} else {
g = gcd(x, y);
}
if (g != 0) x /= g, y /= g;
if (y < 0) x = -x, y = -y;
}
}
// y = 0 の代入も認める
static R raw(T _x, T _y) {
R r;
r.x = _x, r.y = _y;
return r;
}
friend R operator+(const R& l, const R& r) {
if (l.y == r.y) return R{l.x + r.x, l.y};
return R{l.x * r.y + l.y * r.x, l.y * r.y};
}
friend R operator-(const R& l, const R& r) {
if (l.y == r.y) return R{l.x - r.x, l.y};
return R{l.x * r.y - l.y * r.x, l.y * r.y};
}
friend R operator*(const R& l, const R& r) { return R{l.x * r.x, l.y * r.y}; }
friend R operator/(const R& l, const R& r) { return R{l.x * r.y, l.y * r.x}; }
R& operator+=(const R& r) { return (*this) = (*this) + r; }
R& operator-=(const R& r) { return (*this) = (*this) - r; }
R& operator*=(const R& r) { return (*this) = (*this) * r; }
R& operator/=(const R& r) { return (*this) = (*this) / r; }
R operator-() const { return raw(-x, y); }
R inverse() const {
assert(x != 0);
R r = raw(y, x);
if (r.y < 0) r.x = -r.x, r.y = -r.y;
return r;
}
R pow(long long p) const {
R res{1}, base{*this};
while (p) {
if (p & 1) res *= base;
base *= base;
p >>= 1;
}
return res;
}
friend bool operator==(const R& l, const R& r) {
return l.x == r.x && l.y == r.y;
};
friend bool operator!=(const R& l, const R& r) {
return l.x != r.x || l.y != r.y;
};
friend bool operator<(const R& l, const R& r) {
return U{l.x} * r.y < U{l.y} * r.x;
};
friend bool operator<=(const R& l, const R& r) { return l < r || l == r; }
friend bool operator>(const R& l, const R& r) {
return U{l.x} * r.y > U{l.y} * r.x;
};
friend bool operator>=(const R& l, const R& r) { return l > r || l == r; }
friend ostream& operator<<(ostream& os, const R& r) {
os << r.x;
if (r.x != 0 && r.y != 1) os << "/" << r.y;
return os;
}
// T にキャストされるので T が bigint の場合は to_ll も要る
T to_mint(T mod) const {
assert(mod != 0);
T a = y, b = mod, u = 1, v = 0, t;
while (b > 0) {
t = a / b;
swap(a -= t * b, b);
swap(u -= t * v, v);
}
return U((u % mod + mod) % mod) * x % mod;
}
};
using Rational = RationalBase<long long, __int128_t>;
#line 4 "math/rational-binomial.hpp"
template <typename R = Rational>
struct Binomial_rational {
vector<R> fc;
Binomial_rational(int = 0) { fc.emplace_back(1); }
void extend() {
int n = fc.size();
R nxt = fc.back() * n;
fc.push_back(nxt);
}
R fac(int n) {
if (n < 0) return 0;
while ((int)fc.size() <= n) extend();
return fc[n];
}
R finv(int n) {
if (n < 0) return 0;
return fac(n).inverse();
}
R inv(int n) {
if (n < 0) return -inv(-n);
return R{1, max(n, 1)};
}
R C(int n, int r) {
if (n < 0 or r < 0 or n < r) return R{0};
return fac(n) * finv(n - r) * finv(r);
}
R operator()(int n, int r) { return C(n, r); }
template <typename I>
R multinomial(const vector<I>& r) {
static_assert(is_integral<I>::value == true);
int n = 0;
for (auto& x : r) {
if (x < 0) return R{0};
n += x;
}
R res = fac(n);
for (auto& x : r) res *= finv(x);
return res;
}
template <typename I>
R operator()(const vector<I>& r) {
return multinomial(r);
}
R P(int n, int r) {
if (n < 0 || n < r || r < 0) return R(0);
return fac(n) * finv(n - r);
}
// [x^r] 1 / (1-x)^n
R H(int n, int r) {
if (n < 0 || r < 0) return R(0);
return r == 0 ? 1 : C(n + r - 1, r);
}
};
#line 5 "math/rational-fps.hpp"
template <typename R = Rational>
struct FormalPowerSeries_rational : vector<R> {
using vector<R>::vector;
using fps = FormalPowerSeries_rational;
fps &operator+=(const fps &r) {
if (r.size() > this->size()) this->resize(r.size());
for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i];
return *this;
}
fps &operator+=(const R &r) {
if (this->empty()) this->resize(1);
(*this)[0] += r;
return *this;
}
fps &operator-=(const fps &r) {
if (r.size() > this->size()) this->resize(r.size());
for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i];
return *this;
}
fps &operator-=(const R &r) {
if (this->empty()) this->resize(1);
(*this)[0] -= r;
return *this;
}
fps &operator*=(const fps &r) {
int n = this->size() + r.size() - 1;
fps f(n);
for (int i = 0; i < (int)this->size(); i++) {
for (int j = 0; j < (int)r.size(); j++) {
f[i + j] += (*this)[i] * r[j];
}
}
return *this = f;
}
fps &operator*=(const R &v) {
for (int k = 0; k < (int)this->size(); k++) (*this)[k] *= v;
return *this;
}
fps &operator/=(const fps &r) {
if (this->size() < r.size()) {
this->clear();
return *this;
}
int n = this->size() - r.size() + 1;
fps f(*this), g(r);
g.shrink();
R coeff = g.back().inverse();
for (auto &x : g) x *= coeff;
int deg = (int)f.size() - (int)g.size() + 1;
int gs = g.size();
fps quo(deg);
for (int i = deg - 1; i >= 0; i--) {
quo[i] = f[i + gs - 1];
for (int j = 0; j < gs; j++) f[i + j] -= quo[i] * g[j];
}
*this = quo * coeff;
this->resize(n, R(0));
return *this;
}
fps &operator%=(const fps &r) {
*this -= *this / r * r;
shrink();
return *this;
}
fps operator+(const fps &r) const { return fps(*this) += r; }
fps operator+(const R &v) const { return fps(*this) += v; }
fps operator-(const fps &r) const { return fps(*this) -= r; }
fps operator-(const R &v) const { return fps(*this) -= v; }
fps operator*(const fps &r) const { return fps(*this) *= r; }
fps operator*(const R &v) const { return fps(*this) *= v; }
fps operator/(const fps &r) const { return fps(*this) /= r; }
fps operator%(const fps &r) const { return fps(*this) %= r; }
fps operator-() const {
fps ret(this->size());
for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i];
return ret;
}
void shrink() {
while (this->size() && this->back() == R(0)) this->pop_back();
}
fps rev() const {
fps ret(*this);
reverse(begin(ret), end(ret));
return ret;
}
fps dot(fps r) const {
fps ret(min(this->size(), r.size()));
for (int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i];
return ret;
}
// 前 sz 項を取ってくる。sz に足りない項は 0 埋めする
fps pre(int sz) const {
fps ret(begin(*this), begin(*this) + min((int)this->size(), sz));
if ((int)ret.size() < sz) ret.resize(sz);
return ret;
}
fps operator>>(int sz) const {
if ((int)this->size() <= sz) return {};
fps ret(*this);
ret.erase(ret.begin(), ret.begin() + sz);
return ret;
}
fps operator<<(int sz) const {
fps ret(*this);
ret.insert(ret.begin(), sz, R(0));
return ret;
}
fps diff() const {
const int n = (int)this->size();
fps ret(max(0, n - 1));
R one(1), coeff(1);
for (int i = 1; i < n; i++) {
ret[i - 1] = (*this)[i] * coeff;
coeff += one;
}
return ret;
}
fps integral() const {
const int n = (int)this->size();
fps ret(n + 1);
for (int i = 0; i < n; i++) ret[i + 1] = (*this)[i] / (i + 1);
return ret;
}
R eval(R x) const {
R r = 0, w = 1;
for (auto &v : *this) r += w * v, w *= x;
return r;
}
fps inv(int deg = -1) const {
assert((*this)[0] != R(0));
if (deg == -1) deg = (*this).size();
fps ret{R(1) / (*this)[0]};
for (int i = 1; i < deg; i <<= 1) {
ret = (ret + ret - ret * ret * (*this).pre(i << 1)).pre(i << 1);
}
return ret.pre(deg);
}
fps log(int deg = -1) const {
assert(!(*this).empty() && (*this)[0] == R(1));
if (deg == -1) deg = (int)this->size();
return (this->diff() * this->inv(deg)).pre(deg - 1).integral();
}
fps exp(int deg = -1) const {
assert((*this).size() == 0 || (*this)[0] == R(0));
if (deg == -1) deg = (int)this->size();
fps ret{R(1)};
for (int i = 1; i < deg; i <<= 1) {
ret = (ret * (pre(i << 1) + R(1) - ret.log(i << 1))).pre(i << 1);
}
return ret.pre(deg);
}
fps pow(int64_t k, int deg = -1) const {
const int n = (int)this->size();
if (deg == -1) deg = n;
if (k == 0) {
fps ret(deg);
if (deg) ret[0] = 1;
return ret;
}
for (int i = 0; i < n; i++) {
if ((*this)[i] != R(0)) {
R rev = R(1) / (*this)[i];
fps ret = (((*this * rev) >> i).log(deg) * k).exp(deg);
ret *= (*this)[i].pow(k);
ret = (ret << (i * k)).pre(deg);
if ((int)ret.size() < deg) ret.resize(deg, R(0));
return ret;
}
if (__int128_t(i + 1) * k >= deg) return fps(deg, R(0));
}
return fps(deg, R(0));
}
};
math/rational-fps.hpp
math/bigint-all.hpp