library/formalpowerseries/functions/log.hpp
Depends on
Required by
Verified with
Code
#pragma once
#include "../Base.hpp"
#include "./differential.hpp"
#include "./integral.hpp"
namespace fps {
template <typename T, int MX>
FormalPowerSeries<T, MX> log(const FormalPowerSeries<T, MX> &f) {
assert(f.size() and f[0] == 1);
return integral(differential(f) / f);
}
} // namespace fps
#line 1 "library/util/Valarray.hpp"
#include <functional>
#include <ranges>
#include <vector>
template <typename T> struct Valarray : std::vector<T> {
using std::vector<T>::vector; // コンストラクタ継承
Valarray(const std::vector<T> &v) : std::vector<T>(v.begin(), v.end()) {}
private:
template <typename Op>
Valarray &apply_inplace(const Valarray &other, Op op) {
if (this->size() < other.size())
this->resize(other.size(), T(0));
for (auto [a, b] : std::views::zip(*this, other))
a = op(a, b);
return *this;
}
public:
Valarray &operator+=(const Valarray &other) {
return apply_inplace(other, std::plus<>());
}
Valarray &operator-=(const Valarray &other) {
return apply_inplace(other, std::minus<>());
}
Valarray &operator*=(const Valarray &other) {
return apply_inplace(other, std::multiplies<>());
}
Valarray &operator/=(const Valarray &other) {
return apply_inplace(other, std::divides<>());
}
friend Valarray operator+(Valarray a, const Valarray &b) { return a += b; }
friend Valarray operator-(Valarray a, const Valarray &b) { return a -= b; }
friend Valarray operator*(Valarray a, const Valarray &b) { return a *= b; }
friend Valarray operator/(Valarray a, const Valarray &b) { return a /= b; }
Valarray operator-() const {
Valarray g = *this;
for (T &a : g)
a = -a;
return g;
}
};
#line 3 "library/formalpowerseries/Base.hpp"
template <typename T, int MX> struct FormalPowerSeries : Valarray<T> {
using FPS = FormalPowerSeries;
static constexpr int max_size = MX;
using Valarray<T>::Valarray;
using Valarray<T>::size;
using Valarray<T>::resize;
using Valarray<T>::at;
using Valarray<T>::begin;
using Valarray<T>::end;
using Valarray<T>::back;
using Valarray<T>::pop_back;
using value_type = T;
void strict(int n) {
if (size() > n)
resize(n);
shrink();
}
void shrink() {
while (size() and back() == 0)
pop_back();
}
FormalPowerSeries() = default;
FormalPowerSeries(const std::vector<T> &f) : Valarray<T>(f) {
strict(MX);
shrink();
}
static FPS unit() { return {1}; }
static FPS x() { return {0, 1}; }
#pragma region operator
FPS operator-() const { return FPS(Valarray<T>::operator-()); }
FPS &operator+=(const FPS &g) {
Valarray<T>::operator+=(g);
shrink();
return *this;
}
FPS operator+(const FPS &g) const { return FPS(*this) += g; }
FPS &operator-=(const FPS &g) {
Valarray<T>::operator-=(g);
shrink();
return *this;
}
FPS operator-(const FPS &g) const { return FPS(*this) -= g; }
FPS &operator+=(const T &a) {
if (!size())
resize(1);
at(0) += a;
return *this;
}
FPS operator+(const T &a) const { return FPS(*this) += a; }
friend FPS operator+(const T &a, const FPS &f) { return f + a; }
FPS &operator-=(const T &a) {
if (!size())
resize(1);
at(0) -= a;
return *this;
}
FPS operator-(const T &a) { return FPS(*this) -= a; }
friend FPS operator-(const T &a, const FPS &f) { return a + (-f); }
FPS operator*(const FPS &g) const { return FPS(convolution(*this, g)); }
FPS &operator*=(const FPS &g) { return (*this) = (*this) * g; }
FPS &operator*=(const T &a) {
for (size_t i = 0; i < size(); i++)
at(i) *= a;
return *this;
}
FPS operator*(const T &a) const { return FPS(*this) *= a; }
friend FPS operator*(const T &a, const FPS &f) { return f * a; }
FPS operator/(const FPS &g) const { return (*this) * g.inv(); }
FPS &operator/=(const FPS &g) { return (*this) = (*this) / g; }
FPS &operator/=(const T &a) { return *this *= a.inv(); }
FPS operator/(const T &a) { return FPS(*this) /= a; }
FPS &operator<<=(const int d) {
if (d >= MX)
return *this = FPS(0);
resize(std::min(MX, int(size()) + d));
for (int i = int(size()) - 1 - d; i >= 0; i--)
at(i + d) = at(i);
for (int i = d - 1; i >= 0; i--)
at(i) = 0;
return *this;
}
FPS operator<<(const int d) const { return FPS(*this) <<= d; }
FPS &operator>>=(const int d) {
if (d >= size())
return *this = FPS(0);
for (size_t i = d; i < size(); i++)
at(i - d) = at(i);
strict(int(size()) - d);
return *this;
}
FPS operator>>(const int d) const { return FPS(*this) >>= d; }
#pragma endregion operator
FPS pre(int n) const {
if (size() <= n)
return *this;
return FPS(Valarray<T>(this->begin(), this->begin() + n));
}
// 最小の非ゼロ次数(すべて 0 のときは size())を返す
int order() const {
for (int i = 0; i < int(size()); i++) {
if (at(i) != 0)
return i;
}
return int(size());
}
FPS inv(int SZ = MX) const {
assert(size() and at(0) != 0);
FPS res = {at(0).inv()};
for (int n = 1; n < SZ; n <<= 1) {
res *= (2 - this->pre(n << 1) * res);
res.strict(n << 1);
}
res.strict(SZ);
return res;
}
// *this = f_1 + f_2 x^n ⇒ [*this←f_1, return f_2]
FPS separate(int n) {
if (size() <= n)
return FPS(0);
FPS f_2(size() - n);
for (size_t i = n; i < size(); i++)
f_2[i - n] = at(i);
strict(n);
return f_2;
}
T operator()(T a) const {
T res = 0, b = 1;
for (size_t i = 0; i < size(); i++, b *= a)
res += at(i) * b;
return res;
}
};
#line 3 "library/formalpowerseries/functions/differential.hpp"
namespace fps {
template <typename T, int MX>
FormalPowerSeries<T, MX> differential(FormalPowerSeries<T, MX> f) {
if (f.size() <= 1) {
return FormalPowerSeries<T, MX>{};
}
for (std::size_t i = 0; i < f.size() - 1; i++) {
f[i] = (i + 1) * f[i + 1];
}
f.pop_back();
return f;
}
} // namespace fps
#line 3 "library/formalpowerseries/functions/integral.hpp"
namespace fps {
template <typename T, int MX>
FormalPowerSeries<T, MX> integral(FormalPowerSeries<T, MX> f) {
if (f.size() < MX) {
f.resize(f.size() + 1);
}
for (int i = f.size() - 1; i > 0; i--) {
f[i] = f[i - 1] / i;
}
f[0] = 0;
return f;
}
} // namespace fps
#line 5 "library/formalpowerseries/functions/log.hpp"
namespace fps {
template <typename T, int MX>
FormalPowerSeries<T, MX> log(const FormalPowerSeries<T, MX> &f) {
assert(f.size() and f[0] == 1);
return integral(differential(f) / f);
}
} // namespace fps
Back to top page