#include "fft/series.hpp"

View this file on GitHub · Last update: 2026-08-17 04:22:06+05:00

References:

• https://arxiv.org/pdf/2404.05177
• https://mathexp.eu/bostan/publications/BoSc09a.pdf

Depends on

• fft/common.hpp
• fft/engine.hpp
• fft/multiply.hpp
• fft/series_core.hpp

Required by

• fft/all.hpp
• 🧪 fft/poly.test.cpp
• 🧪 fft/series.test.cpp

Verified with

• verify/composition_of_formal_power_series.test.cpp
• verify/composition_of_formal_power_series_large.test.cpp
• verify/exp_of_formal_power_series.test.cpp
• verify/kth_term_of_linearly_recurrent_sequence.test.cpp
• verify/log_of_formal_power_series.test.cpp
• verify/pow_of_formal_power_series.test.cpp

Code

	Coverage	Exec / Excl / Total
Lines	84.9%	203 / 0 / 239
Functions	100.0%	23 / 0 / 23
Branches	74.9%	555 / 0 / 741
Full report

#pragma once

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

#include "fft/series_core.hpp"

// ==== analytic ops ====
// Free functions over series-like operands; each borrows the operand's span
// and writes a fresh result.
// TODO: reuse/populate the operands' whole/prefix transform caches

namespace ecnerwala::series {

template <like S>
vec<typename S::engine_t, S::exact_v> stretch(const S& a_, int n) {
	using E = typename S::engine_t;
	span<E, S::exact_v> a = a_;
	vec<E, S::exact_v> r(size_t(a.len()));
	for (int i = 0; i*n < a.len(); i++) {
		r[i*n] = a[i];
	}
	return r;
}
template <like S>
vec<typename S::engine_t, S::exact_v> deriv_shift(const S& a_) {
	using E = typename S::engine_t;
	span<E, S::exact_v> a = a_;
	vec<E, S::exact_v> r(a.begin(), a.end());
	for (int i = 0; i < r.len(); i++) {
		r[i] *= i;
	}
	return r;
}
template <like S>
vec<typename S::engine_t, S::exact_v> integ_shift(const S& a_) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, S::exact_v> a = a_;
	assert(a[0] == 0);
	vec<E, S::exact_v> r(a.begin(), a.end());
	T f = 1;
	for (int i = 1; i < r.len(); i++) {
		r[i] *= f;
		f *= i;
	}
	f = inv(f);
	for (int i = r.len() - 1; i > 0; i--) {
		r[i] *= f;
		f *= i;
	}
	return r;
}
template <like S>
vec<typename S::engine_t, S::exact_v> integ_shift_offset(const S& a_, int offset) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, S::exact_v> a = a_;
	vec<E, S::exact_v> r(a.begin(), a.end());
	T f = 1;
	for (int i = 0; i < r.len(); i++) {
		r[i] *= f;
		f *= i + offset;
	}
	assert(f != 0);
	f = inv(f);
	for (int i = r.len() - 1; i >= 0; i--) {
		r[i] *= f;
		f *= i + offset;
	}
	return r;
}
template <trunc_like S>
trunc<typename S::engine_t> deriv_shift_log(const S& a) {
	return deriv_shift(a) * ps_inv(a);
}
template <trunc_like S>
trunc<typename S::engine_t> ps_log(const S& a) {
	assert(a[0] == 1);
	return integ_shift(deriv_shift_log(a));
}
template <trunc_like S>
trunc<typename S::engine_t> ps_exp(const S& a_) {
	// See https://mathexp.eu/bostan/publications/BoSc09a.pdf for details
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, false> a = a_;
	assert(a.len() >= 1);
	assert(a[0] == 0);
	trunc<E> r(1, T(1)); r.reserve(size_t(a.len()));
	trunc<E> invR(1, T(1)); invR.reserve(size_t(a.len()));
	while (r.len() < a.len()) {
		int o_sz = r.len();
		int n_sz = std::min(o_sz * 2, a.len());
		trunc<E> t = deriv_shift(trunc<E>(a.begin(), a.begin() + o_sz));
		fft::multiply_circular<E>(std::span<const T>(t), std::span<const T>(r).first(o_sz), std::span<T>(t), o_sz);
		t = deriv_shift(r) - t;
		t *= invR;
		t.resize(size_t(n_sz - o_sz));
		trunc<E> v(a.begin() + o_sz, a.begin() + n_sz);
		v -= integ_shift_offset(t, o_sz);
		v *= r;
		r.resize(size_t(n_sz));
		std::copy(v.begin(), v.end(), r.begin() + o_sz);
		if (r.len() < a.len()) {
			// double invR via a Newton step
			assert(r.len() == 2 * invR.len());
			int n = invR.len();
			int nn = r.len();
			trunc<E> tmp(size_t(4) * n);
			fft::square<E>(std::span<const T>(invR).first(n), std::span<T>(tmp));
			fft::multiply<E>(std::span<const T>(tmp).first(nn), std::span<const T>(r).first(nn), std::span<T>(tmp));
			invR.resize(size_t(nn));
			for (int i = n; i < nn; i++) invR[i] = -tmp[i];
		}
	}
	return r;
}
template <trunc_like S>
trunc<typename S::engine_t> ps_pow_monic(const S& a_, typename S::engine_t::value_type k) {
	using E = typename S::engine_t;
	span<E, false> a = a_;
	if (a.len() == 0) return {};
	assert(a[0] == 1);
	trunc<E> l = ps_log(a_);
	l *= k;
	return ps_exp(l);
}
template <trunc_like S>
trunc<typename S::engine_t> ps_pow(const S& a_, int64_t k) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, false> a = a_;
	assert(k >= 0);
	if (k == 0) {
		trunc<E> r(size_t(a.len()), T(0));
		if (r.len() > 0) r[0] = T(1);
		return r;
	}

	int st = 0;
	while (st < a.len() && a[st] == 0) st++;

	if (st > 0 && k > (a.len() - 1) / st) {
		return trunc<E>(size_t(a.len()), T(0));
	}

	trunc<E> r(a.begin() + st, a.end() - (st * (k-1)));
	T leading_coeff = r[0];
	r *= inv(leading_coeff);
	r = ps_pow_monic(r, T(k));
	r *= power(leading_coeff, k);
	r.insert(r.begin(), size_t(st * k), T(0));
	assert(r.len() == a.len());
	return r;
}

template <trunc_like S>
trunc<typename S::engine_t> to_newton_sums(const S& a, int deg) {
	auto r = deriv_shift_log(a);
	r[0] = deg;
	for (int i = 1; i < r.len(); i++) r[i] = -r[i];
	return r;
}
template <trunc_like S>
trunc<typename S::engine_t> from_newton_sums(const S& s_, int deg) {
	using E = typename S::engine_t;
	span<E, false> s = s_;
	assert(s[0] == deg);
	trunc<E> r(s.begin(), s.end());
	r[0] = 0;
	for (int i = 1; i < r.len(); i++) r[i] = -r[i];
	return ps_exp(integ_shift(std::move(r)));
}

// Calculates prod 1/(1-x^i)^{a[i]}
template <trunc_like S>
trunc<typename S::engine_t> euler_transform(const S& a) {
	using E = typename S::engine_t;
	trunc<E> r = deriv_shift(a);
	std::vector<bool> is_prime(size_t(r.len()), true);
	for (int p = 2; p < r.len(); p++) {
		if (!is_prime[p]) continue;
		for (int i = 1; i*p < r.len(); i++) {
			r[i*p] += r[i];
			is_prime[i*p] = false;
		}
	}
	return ps_exp(integ_shift(r));
}
template <trunc_like S>
trunc<typename S::engine_t> inverse_euler_transform(const S& a) {
	using E = typename S::engine_t;
	trunc<E> r = deriv_shift(ps_log(a));
	std::vector<bool> is_prime(size_t(r.len()), true);
	for (int p = 2; p < r.len(); p++) {
		if (!is_prime[p]) continue;
		for (int i = (r.len()-1)/p; i >= 1; i--) {
			r[i*p] -= r[i];
			is_prime[i*p] = false;
		}
	}
	return integ_shift(r);
}

// Helper packed bivariate buffer for Kinoshita-Li composition (arXiv:2404.05177).
//
// The motivation is performing Bostan-Mori (Graeffe root-squaring) to compute
// something like [x^n] P / Q_0(x, y) with deg_y(Q_0) = 1 and deg_x(Q_0) = n.
//
// In each step, we want to compute Q_{i+1}(x^2, y) = Q_i(x, y) * Q_i(-x, y).
// This doubles the degree of y and also lets us truncate x at half the previous
// degree, leaving the total size invariant.
//
// We will store Q as a packed buffer with x as the inner dimension to facilitate easy Q(-x) substitution.
// The inner span will be 2*deg(x), and the outer span will be 2*deg(y).
// As we advance, we will also return the cached transform of Q_i(-x, y) for the caller to use in the numerator.
template <fft::engine E> struct packed_bivariate {
	using T = typename E::value_type;
	int L, l;
	std::vector<T> c;

	// Q_0 = 1 - y g(x), deg g < n <= 2^L
	packed_bivariate(int L_, std::span<const T> g) : L(L_), l(0), c(size_t(4) << L) {
		c[0] = T(1);
		for (int i = 0; i < sz(g); i++) c[(2 << L) + i] = -g[i];
	}

	fft::transformed<E> advance() {
		int B = 4 << L;
		auto tq = E::transform(std::span<const T>(c), B);
		auto tn = E::negate_arg(tq, B);
		E::finish(
			E::downsample(E::mul(tq, tn, B), B/2, false),
			std::span<T>(c).first(B/2)
		);
		l++;
		// undo the circular wraparound using monicity in y
		for (int i = 0; i < (2 << (L - l)); i++) {
			c[(2 << L) + i] = c[i];
			c[i] = T(0);
		}
		c[2 << L] -= T(1);
		c[0] = T(1);
		// zero x coefficients beyond the level's truncation mod x^(2^(L-l))
		std::fill(c.begin() + (2 << L) + (1 << (L - l)), c.end(), T(0));
		for (int i = 0; i < (2 << L); i += 2 << (L - l)) {
			for (int j = 0; j < (1 << (L - l)); j++) {
				c[i + (1 << (L - l)) + j] = T(0);
			}
		}
		return tn;
	}
};

// Calculates f(g(x)) mod x^n where deg(g) == n
template <trunc_like SF, trunc_like SG> requires fft::same_engine<SF, SG>
trunc<typename SF::engine_t> ps_compose(const SF& f_, const SG& g_) {
	using E = typename SF::engine_t;
	using T = typename E::value_type;
	span<E, false> f = f_;
	span<E, false> g = g_;
	if (g.len() == 0) return {};

	int m = f.len();
	int n = g.len();

	// https://arxiv.org/pdf/2404.05177
	// Consider P(y) = f(1/y) has terms from y^{-(m-1)}...y^0 (Laurent series)
	// We want [y^0] P(y) / (1 - y g(x))
	// Let Q_0 = 1 - yg(x)
	// Q_{i+1}(x^2, y) = Q_i(x, y) * Q_i(-x, y) mod x^{ceil(n / 2^i)}
	// deg_y(Q_i) = 2^i, deg_x(Q_i) = ceil(n / 2^i) - 1
	//
	// [y^0] P(y) / Q_l(x^2^l, y) * Q_{l-1}(-x^2^{l-1}, y) * Q_{l-2}(-x^2^{l-2}, y) * ... * Q_0(-x, y)
	// The total y deg of Q_{k-1} ... Q_0 is 2^k-1
	int L = __builtin_ctz(unsigned(nextPow2(n)));
	int B = 4 << L;
	packed_bivariate<E> Q(L, g.coeffs());
	// tneg[l] is the transform of Q_l(-x, y), reused by the pushdown pass below
	std::vector<fft::transformed<E>> tneg;
	tneg.reserve(L);
	for (int l = 1; l <= L; l++) tneg.push_back(Q.advance());
	trunc<E> P;
	{
		P = trunc<E>(f.begin(), f.end());
		std::reverse(P.begin(), P.end());
		trunc<E> QL((1 << L) + 1);
		for (int i = 0; i <= (1 << L); i++) {
			QL[i] = Q.c[2 * i];
		}
		QL.resize(size_t(m), T(0));
		P *= ps_inv(QL);
		std::reverse(P.begin(), P.end());
		P.resize(size_t(1) << L, T(0));
		std::reverse(P.begin(), P.end());
		P.resize(size_t(B), T(0));
		for (int i = (1 << L) - 1; i > 0; i--) {
			P[2*i] = P[i];
			P[i] = T(0);
		}
	}
	for (int l = L-1; l >= 0; l--) {
		// Spread it out, clear the high terms
		for (int i = (2 << L) - 1; i > 0; i--) {
			T v = P[i];
			P[2*i] = ((2*i) & (1 << (L-l))) ? T(0) : v;
			P[i] = T(0);
		}
		auto tp = E::transform(std::span<const T>(P), B);
		E::finish(E::mul(tneg[l], tp, B), std::span<T>(P));
		for (int i = 0; i < (2 << L); i++) {
			P[i] = P[(2 << L) + i];
			P[(2 << L) + i] = T(0);
		}
	}
	return trunc<E>(P.begin(), P.begin() + n);
}

// [x^k] p(x)/q(x) (Bostan-Mori) for an exact rational function.
template <exact_like P, exact_like Q> requires fft::same_engine<P, Q>
P::engine_t::value_type kth_term_of_rational_function(
	const P& p,
	const Q& q,
	uint64_t k
) {
	using E = P::engine_t;
	using T = E::value_type;

	assert(q.len() > 0 && q[0] != T(0));
	// Check this here so we avoid accessing p[0]
	if (p.len() == 0) return T(0);

	// Size up in a pretty conservative way
	int d = std::max(p.len() + 1, q.len());
	assert(d >= 2);

	int n = nextPow2((d-1) + d - 1); // >= d

	// Seed the loop transforms from any whole caches; the buffers below hold the
	// current p, q (zero-padded, which extend_to tolerates).
	fft::transformed<E> tq, tp;
	if (auto cq = detail::cache_of(q)) { E::extend_to(cq->get(), n, q); tq = cq->get(); }
	if (auto cp = detail::cache_of(p)) { E::extend_to(cp->get(), n, p); tp = cp->get(); }

	std::vector<T> p_buf(d-1, T(0));
	std::ranges::copy(std::span<const T>(p), p_buf.begin());
	std::vector<T> q_buf(d, T(0));
	std::ranges::copy(std::span<const T>(q), q_buf.begin());

	while (k > 0) {
		E::extend_to(tq, n, q_buf);
		auto tnq = E::negate_arg(tq, n);
		E::extend_to(tp, n, p_buf);

		// P <- downsample(P(x) * Q(-x))
		auto ntp = E::downsample(E::mul(tp, tnq, n), n/2, bool(k & 1));
		assert(ntp.size() == n/2);
		if constexpr (std::same_as<typename E::product, typename E::transformed>) {
			tp = ntp;
		} else {
			tp = {};
		}
		E::finish(std::move(ntp), std::span(p_buf));
		k >>= 1;

		// Save the last iteration if we're done
		if (!k) {
			// HACK: fix the constant coefficient of q only
			q_buf[0] *= q_buf[0];
			break;
		}

		// Q <- downsample(Q(x) * Q(-x))
		auto ntq = E::downsample(E::mul(tq, tnq, n), n/2, false);
		assert(ntq.size() == n/2);
		if constexpr (std::same_as<typename E::product, typename E::transformed>) {
			tq = ntq;
		} else {
			tq = {};
		}
		if (n/2 == d-1) {
			// Fix the wraparound
			T v0 = q_buf[0] * q_buf[0];
			E::finish(std::move(ntq), std::span(q_buf).first(d-1));
			q_buf[d-1] = std::exchange(q_buf[0], v0) - v0;
		} else {
			E::finish(std::move(ntq), std::span(q_buf));
		}
	}
	return p_buf[0] * inv(q_buf[0]);
}

// Find the kth term of linearly recurrent sequence S with char poly Q and len(S) >= len(Q)-1
template <trunc_like S, exact_like Q> requires fft::same_engine<S, Q>
S::engine_t::value_type kth_term_of_linear_recurrence(
	const S& s,
	const Q& q,
	uint64_t k
) {
	using E = S::engine_t;
	using T = E::value_type;

	assert(q.len() > 0 && q[0] != T(0));
	assert(s.len() >= q.len()-1);

	// Don't even bother with P so we don't have to do truncation checks
	// TODO: Could use generic multiply for this whole part?
	fft::transformed<E> tq;
	auto q_cached = detail::as_cached_span(q, tq);

	// Compute the prefix and then hard-cast it to exact
	span<E, false> sv = s;
	auto p = exact<E>(sv.first(q.len()-1) * q_cached);
	return kth_term_of_rational_function(p, q_cached, k);
}

/* namespace ecnerwala::series */ }

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <span>
#include <utility>
#include <vector>
#include <concepts>
#include <cstddef>
#include <functional>
#include <optional>
#include <type_traits>
#include <iterator>
#line 2 "src/fft/series.hpp"

#line 9 "src/fft/series.hpp"

#line 2 "src/fft/series_core.hpp"

#line 12 "src/fft/series_core.hpp"

#line 2 "src/fft/multiply.hpp"

#line 10 "src/fft/multiply.hpp"

#line 2 "src/fft/engine.hpp"

#line 7 "src/fft/engine.hpp"

#line 2 "src/fft/common.hpp"

#line 8 "src/fft/common.hpp"

/**
 * Author: Andrew He
 * Source: http://neerc.ifmo.ru/trains/toulouse/2017/fft2.pdf
 * Papers about accuracy: http://www.daemonology.net/papers/fft.pdf, http://www.cs.berkeley.edu/~fateman/papers/fftvsothers.pdf
 * For integers rounding works if $(|a| + |b|)\max(a, b) < \mathtt{\sim} 10^9$, or in theory maybe $10^6$.
 *
 * Abstraction layers:
 *   fft_core<num>   FFT itself and other ops on rings with 2^k-th roots of unity. We use bit-reversed indexing in the frequency domain.
 *
 *   engines         Engines for packing/unpacking arbitrary rings for convolution (the `engine` concept).
 *                   Still expose (opaque) transform-domain objects for caching/fusion.
 *
 *   multiply layer  Wrappers for convolving bounded sequences: track length/truncation.
 *
 *   value types     series::vec<E, exact> - R[[x]]
 *                   series::exact<E> - exact (finite-support) power series
 *                   series::trunc<E> - truncated prefix of an (infinite) power series
 *
 *                   polynomials - R[x]. Under x -> 1/x a polynomial becomes a Laurent polynomial in 1/x;
 *                                 shifting by x^{deg P} (reversal) lands it in R[[x]], and we store that exact series.
 *                   poly::vec<E> - polynomial type, supporting natural indexing
 *                   poly::form<E> - finite-support linear forms, via the pairing <P, S> = [x^0] P(1/x) S(x)
 *                                    a linear form is one side of this pairing, applied to the other
 *
 *                   online_multiplier<E> - online (relaxed) multiplication of 2 sequences in n log^2 n time
 *                   ap_sampled_poly<E> - a polynomial stored as its evaluations on an arithmetic progression
 */

namespace ecnerwala {

template<class T> int sz(T&& arg) { using std::size; return int(size(std::forward<T>(arg))); }
inline int nextPow2(int s) { return 1 << (s > 1 ? 32 - __builtin_clz(s-1) : 0); }

namespace fft {

using std::swap;
using std::vector;
using std::min;
using std::max;

// Reusable scratch buffers. Not thread-safe by default: this is deliberately plain
// static storage so single-threaded programs pay no TLS indirection; define
// ECNERWALA_FFT_POOL_STORAGE to `thread_local` for multithreaded use.
#ifndef ECNERWALA_FFT_POOL_STORAGE
#define ECNERWALA_FFT_POOL_STORAGE
#endif
template <typename T> struct buffer_pool {
	static inline ECNERWALA_FFT_POOL_STORAGE std::vector<std::vector<T>> free_list;
	struct handle {
		std::vector<T> v;
		explicit handle(int n) {
			if (!free_list.empty()) {
				v = std::move(free_list.back());
				free_list.pop_back();
			}
			v.assign(n, T());
		}
		handle(const handle&) = delete;
		handle& operator=(const handle&) = delete;
		handle(handle&& o) noexcept : v(std::move(o.v)) {}
		~handle() { if (v.capacity()) free_list.push_back(std::move(v)); }
		T& operator[](int i) { return v[i]; }
		operator std::span<T>() { return std::span<T>(v); }
		std::span<T> span() { return std::span<T>(v); }
	};
	static handle get(int n) { return handle(n); }
};

/* namespace fft */ }

/* namespace ecnerwala */ }
#line 9 "src/fft/engine.hpp"

namespace ecnerwala::fft {

// ==== engine concept ====

// Output operations for the finish step to express arbitrary fusion into the output buffer.
struct assign_op { template <typename T> void operator()(T& d, T v) const { d = v; } };
struct add_op { template <typename T> void operator()(T& d, T v) const { d += v; } };
struct sub_op { template <typename T> void operator()(T& d, T v) const { d -= v; } };
struct add_twice_op { template <typename T> void operator()(T& d, T v) const { d += v + v; } };

// `engine` contract
//   engine represents a way of packing/unpacking sequences over an arbitrary ring into FFT-style transforms.
//   We expect transforms/products of transforms to be linear but potentially lossy/imprecise, so we'll track precision
//   at compile-time as a template parameter.
//
//   E::value_type   The ring we operate over
//   E::unit_scale   0 or 1 depending on whether there's error that can accumulate
//   E::commutative  A marker for whether the ring is commutative
//
//   transformed_t<A>  The transform of a sequence. This object owns its data buffer.
//   product_t<A>      The product of 2 transforms. May equal transformed_t, particularly when unit_scale = 0.
//
//   transformed       alias for transformed_t<unit_scale>
//   product           alias for product_t<unit_scale>
//
//   The basic multiplication API is
//      transform(span<const value_type> in, int n) -> transformed_t<unit_scale>
//      mul(transformed_t<A>, transformed_t<B>, int n) -> product_t<A*B>
//      mul2(a1, b1, a2, b2, int n) -> product_t<A1*B1 + A2*B2>, computing a1*b1 + a2*b2 in one pass
//      finish(product_t<A>, span<value_type>& out, Op) -> void
//
//   Input span can be length up to 2n.
//   Output spans can be length up to n; only the prefix that exists is filled.
//   finish applies Op exactly once per out element, in index order, to
//   value_type targets (so ops may be stateful).
//   Transforms can be longer than necessary, and only the relevant prefix is used.
//
//   For non-exact engines, there's some subtlety in whether we wrap before or after packing.
//   We will choose to wrap *after* packing, which hurts error bounds but makes the prefix condition more uniform.
//
//   Additionally, we have APIs to take advantage of linearity in both transformed and product space:
//      add(transformed_t<A>, transformed_t<B>) -> transformed_t<A+B>
//      add(product_t<K1>, product_t<K2>) -> product_t<K1+K2>
//
//   Finally, we expose some additional fast-transform optimization paths.
//   extend_to only operates on transformed_t<unit_scale>; the others are scale-generic.
//   downsample is also defined on product_t (halving before finish saves inverse-transform work).
//      extend_to       build (if empty) or grow a transform to size m by repeated doubling; feed the same coefficients (sz <= 2m) every time
//                      (each doubling step reads only the coefficients that fit, so zero-padded buffers are fine:
//                       coefficients past twice the existing transform's size must be zero, or it couldn't be a prefix)
//      downsample      compute the half-sized transform/product of just the even (odd = false) or odd terms of the input
//      negate_arg      size n transform of A(-x)
template <typename E>
concept engine = requires(
	std::span<const typename E::value_type> in,
	std::span<typename E::value_type> out,
	typename E::transformed& t,
	const typename E::transformed& ct,
	typename E::product& p,
	const typename E::product& cp,
	int n
) {
	typename E::value_type;
	{ E::transform(in, n) } -> std::same_as<typename E::transformed>;
	{ ct.size() } -> std::same_as<int>;
	E::extend_to(t, n, in);
	{ E::downsample(ct, n, false) } -> std::same_as<typename E::transformed>;
	{ E::downsample(cp, n, false) } -> std::same_as<typename E::product>;
	{ E::negate_arg(ct, n) } -> std::same_as<typename E::transformed>;
	{ E::mul(ct, ct, n) } -> std::same_as<typename E::product>;
	{ E::sq(ct, n) } -> std::same_as<typename E::product>;
	{ E::mul2(ct, ct, ct, ct, n) } -> std::same_as<typename E::template product_t<2 * E::unit_scale>>;
	E::finish(std::move(p), out);
	E::finish(std::move(p), out, add_op{});
	E::finish(E::add(std::move(p), std::move(p)), out);
	{ E::add(E::transform(in, n), ct) } -> std::same_as<typename E::template transformed_t<2 * E::unit_scale>>;
	{ E::add(std::move(p), std::move(p)) } -> std::same_as<typename E::template product_t<2 * E::unit_scale>>;
	requires std::same_as<std::remove_cvref_t<decltype(E::commutative)>, bool>;
	requires std::same_as<std::remove_cvref_t<decltype(E::unit_scale)>, int>;
};

// Constrains two engine-parameterized value types to share the same engine.
template <typename A, typename B>
concept same_engine = std::same_as<typename A::engine_t, typename B::engine_t>;

// short spelling for E::transformed at use sites
template <engine E> using transformed = typename E::transformed;

/* namespace ecnerwala::fft */ }
#line 12 "src/fft/multiply.hpp"

namespace ecnerwala::fft {

// ==== multiply layer ====
// These are free functions to convolve spans.
//
// The interfaces will typically take input spans, an output span, and an Op representing how to fold the result into the output.
// Output spans may alias one of the input spans.
// Output spans may be shorter than expected; the output will just be truncated.
//
// Some functions may also take E::transformed& objects associated with the input
// spans. These will be lazily filled (see E::extend_to) and used if available.

// Circular convolution mod n (power of 2)
template <engine E, typename Op = assign_op>
void multiply_circular(std::span<const typename E::value_type> a, std::span<const typename E::value_type> b,
		std::span<typename E::value_type> out, int n, Op op = {}) {
	assert(!(n & (n-1)));
	auto ta = E::transform(a, n);
	auto tb = E::transform(b, n);
	E::finish(E::mul(ta, tb, n), out, op);
}

template <engine E, typename Op = assign_op>
void square_circular(std::span<const typename E::value_type> a, std::span<typename E::value_type> out, int n, Op op = {}) {
	assert(!(n & (n-1)));
	auto ta = E::transform(a, n);
	E::finish(E::sq(ta, n), out, op);
}

namespace detail {
// Arrays of length 2^k + 1 are somewhat common, so we will optimize them by
// multiplying mod 2^k, and fixing up the leading coefficient.

// Helpers to detect and perform this optimization.
struct conv_size { int n; bool cut; };
inline conv_size conv_size_for(int s) {
	int n = nextPow2(s);
	bool cut = (n == 2 * (s - 1));
	return {cut ? n / 2 : n, cut};
}

// Call op while lazily applying the correction if necessary
template <typename T, typename Op>
void emit_linear(std::span<T> buf, int n, int s, bool cut, T c0, std::span<T> out, Op op) {
	T cn{};
	if (cut) {
		cn = buf[0] - c0;
		buf[0] = c0;
	}
	int lim = min(sz(out), min(s, n));
	for (int i = 0; i < lim; i++) op(out[i], buf[i]);
	if (cut && sz(out) >= s) op(out[s-1], cn);
}

// Applies op, diverting the wrapped leading coefficient of a cut product:
// out[0] receives c0 and the wraparound term is captured into cn for the
// caller to emit at out[s-1].
template <typename T, typename Op>
struct cut_op {
	Op op;
	T* out0;
	T c0;
	T& cn;
	void operator()(T& x, T v) const {
		if (&x == out0) { cn = v - c0; v = c0; }
		op(x, v);
	}
};

// finish + emit_linear fused: write the finished product directly into out,
// applying the cut correction in place.
template <engine E, typename P, typename Op = assign_op>
void finish_linear(
	P&& p, int n, int s, bool cut,
	typename E::value_type c0, std::span<typename E::value_type> out, Op op = {}
) {
	using T = typename E::value_type;
	if (sz(out) == 0) return;
	int lim = min(sz(out), min(s, n));
	if (!cut) {
		E::finish(std::move(p), out.subspan(0, lim), op);
	} else {
		T cn{};
		E::finish(std::move(p), out.subspan(0, lim), cut_op<T, Op>{op, &out[0], c0, cn});
		if (sz(out) >= s) op(out[s-1], cn);
	}
}

}

template <engine E, typename Op = assign_op>
void multiply(std::span<const typename E::value_type> a, std::span<const typename E::value_type> b,
		std::span<typename E::value_type> out, Op op = {}) {
	using T = typename E::value_type;
	if (sz(a) == 0 || sz(b) == 0) return;
	int s = sz(a) + sz(b) - 1;
	auto [n, cut] = detail::conv_size_for(s);
	T c0 = a[0] * b[0];
	auto buf = buffer_pool<T>::get(n);
	multiply_circular<E>(a, b, buf.span(), n);
	detail::emit_linear<T>(buf.span(), n, s, cut, c0, out, op);
}

template <engine E, typename Op = assign_op>
void multiply(std::span<const typename E::value_type> a, transformed<E>& ta,
		std::span<const typename E::value_type> b, transformed<E>& tb,
		std::span<typename E::value_type> out, Op op = {}) {
	using T = typename E::value_type;
	if (sz(a) == 0 || sz(b) == 0) return;
	int s = sz(a) + sz(b) - 1;
	auto [n, cut] = detail::conv_size_for(s);
	T c0 = a[0] * b[0];
	E::extend_to(ta, n, a);
	E::extend_to(tb, n, b);
	detail::finish_linear<E>(E::mul(ta, tb, n), n, s, cut, c0, out, op);
}

template <engine E, typename Op = assign_op>
void multiply_add2(std::span<const typename E::value_type> a1, transformed<E>& ta1,
		std::span<const typename E::value_type> b1, transformed<E>& tb1,
		std::span<const typename E::value_type> a2, transformed<E>& ta2,
		std::span<const typename E::value_type> b2, transformed<E>& tb2,
		std::span<typename E::value_type> out, Op op = {}) {
	using T = typename E::value_type;
	assert(sz(a1) > 0 && sz(b1) > 0 && sz(a2) > 0 && sz(b2) > 0);
	int s = sz(a1) + sz(b1) - 1;
	assert(sz(a2) + sz(b2) - 1 == s);
	auto [n, cut] = detail::conv_size_for(s);
	T c0 = a1[0] * b1[0] + a2[0] * b2[0];
	E::extend_to(ta1, n, a1); E::extend_to(tb1, n, b1);
	E::extend_to(ta2, n, a2); E::extend_to(tb2, n, b2);
	detail::finish_linear<E>(E::mul2(ta1, tb1, ta2, tb2, n), n, s, cut, c0, out, op);
}

// As multiply_add2, but also outputs the summed pointwise product as a reusable
// transform of the (full-length) result, like multiply_cached.
template <engine E>
void multiply_add2_cached(
		std::span<const typename E::value_type> a1, transformed<E>& ta1,
		std::span<const typename E::value_type> b1, transformed<E>& tb1,
		std::span<const typename E::value_type> a2, transformed<E>& ta2,
		std::span<const typename E::value_type> b2, transformed<E>& tb2,
		std::vector<typename E::value_type>& coeffs, transformed<E>& t) {
	using T = typename E::value_type;
	assert(sz(a1) > 0 && sz(b1) > 0 && sz(a2) > 0 && sz(b2) > 0);
	int s = sz(a1) + sz(b1) - 1;
	assert(sz(a2) + sz(b2) - 1 == s);
	coeffs.assign(size_t(s), T{});
	t = transformed<E>{};
	if constexpr (std::same_as<typename E::product, transformed<E>>) {
		auto [n, cut] = detail::conv_size_for(s);
		T c0 = a1[0] * b1[0] + a2[0] * b2[0];
		E::extend_to(ta1, n, a1); E::extend_to(tb1, n, b1);
		E::extend_to(ta2, n, a2); E::extend_to(tb2, n, b2);
		auto p = E::mul2(ta1, tb1, ta2, tb2, n);
		auto tp = p;
		detail::finish_linear<E>(std::move(p), n, s, cut, c0, std::span<T>(coeffs));
		t = std::move(tp);
	} else {
		multiply_add2<E>(a1, ta1, b1, tb1, a2, ta2, b2, tb2, std::span<T>(coeffs));
	}
}

// This helper also accepts an output transform which will be populated if it is cheap to do so
template <engine E>
void multiply_cached(std::span<const typename E::value_type> a, transformed<E>& ta,
		std::span<const typename E::value_type> b, transformed<E>& tb,
		std::vector<typename E::value_type>& coeffs, transformed<E>& t) {
	using T = typename E::value_type;
	coeffs.assign(size_t(sz(a) && sz(b) ? sz(a) + sz(b) - 1 : 0), T{});
	t = transformed<E>{};
	if (coeffs.empty()) return;
	int s = sz(coeffs);
	if constexpr (std::same_as<typename E::product, transformed<E>>) {
		auto [n, cut] = detail::conv_size_for(s);
		T c0 = a[0] * b[0];
		E::extend_to(ta, n, a);
		E::extend_to(tb, n, b);
		auto p = E::mul(ta, tb, n);
		auto tp = p;
		detail::finish_linear<E>(std::move(p), n, s, cut, c0, std::span<T>(coeffs));
		t = std::move(tp);
	} else {
		multiply<E>(a, ta, b, tb, std::span<T>(coeffs));
	}
}

template <engine E, typename Op = assign_op>
void square(std::span<const typename E::value_type> a, std::span<typename E::value_type> out, Op op = {}) {
	using T = typename E::value_type;
	if (sz(a) == 0) return;
	int s = 2 * sz(a) - 1;
	auto [n, cut] = detail::conv_size_for(s);
	T c0 = a[0] * a[0];
	auto buf = buffer_pool<T>::get(n);
	square_circular<E>(a, buf.span(), n);
	detail::emit_linear<T>(buf.span(), n, s, cut, c0, out, op);
}

template <engine E, typename Op = assign_op>
void square(std::span<const typename E::value_type> a, transformed<E>& ta,
		std::span<typename E::value_type> out, Op op = {}) {
	using T = typename E::value_type;
	if (sz(a) == 0) return;
	int s = 2 * sz(a) - 1;
	auto [n, cut] = detail::conv_size_for(s);
	T c0 = a[0] * a[0];
	E::extend_to(ta, n, a);
	detail::finish_linear<E>(E::sq(ta, n), n, s, cut, c0, out, op);
}

// As square, but also outputs the pointwise product as a reusable transform of
// the result (empty when the engine's product isn't a transform).
template <engine E>
void square_cached(std::span<const typename E::value_type> a, transformed<E>& ta,
		std::vector<typename E::value_type>& coeffs, transformed<E>& t) {
	using T = typename E::value_type;
	coeffs.assign(size_t(sz(a) ? 2 * sz(a) - 1 : 0), T{});
	t = transformed<E>{};
	if (coeffs.empty()) return;
	int s = sz(coeffs);
	if constexpr (std::same_as<typename E::product, transformed<E>>) {
		auto [n, cut] = detail::conv_size_for(s);
		T c0 = a[0] * a[0];
		E::extend_to(ta, n, a);
		auto p = E::sq(ta, n);
		auto tp = p;
		detail::finish_linear<E>(std::move(p), n, s, cut, c0, std::span<T>(coeffs));
		t = std::move(tp);
	} else {
		square<E>(a, ta, std::span<T>(coeffs));
	}
}

template <engine E> vector<typename E::value_type> multiply(
		const vector<typename E::value_type>& a, const vector<typename E::value_type>& b) {
	using T = typename E::value_type;
	if (sz(a) == 0 || sz(b) == 0) return {};
	vector<T> r(sz(a) + sz(b) - 1);
	multiply<E>(std::span<const T>(a), std::span<const T>(b), std::span<T>(r));
	return r;
}

template <engine E> vector<typename E::value_type> square(const vector<typename E::value_type>& a) {
	using T = typename E::value_type;
	if (sz(a) == 0) return {};
	vector<T> r(2 * sz(a) - 1);
	square<E>(std::span<const T>(a), std::span<T>(r));
	return r;
}

namespace detail {
// emit_linear but for middle_product
template <typename T, typename Op>
void emit_middle(std::span<T> buf, bool cut, int la, int lb, T c0, T ctop, std::span<T> out, Op op) {
	int m = la - lb + 1;
	T cn{};
	if (cut) {
		cn = buf[0] - c0; // for lb == 1 these coincide: slot 0 = c_0 + c_n and ctop = c_n
		buf[lb - 1] -= ctop;
	}
	int lim = min(sz(out), cut ? m - 1 : m);
	for (int t = 0; t < lim; t++) op(out[t], buf[lb - 1 + t]);
	if (cut && sz(out) >= m) op(out[m-1], cn);
}
}

// Middle product (the transposed multiplication): takes only coefficients of a * b which include terms from all of b.
// Must have len(a) >= len(b)
template <engine E, typename Op = assign_op>
void middle_product(std::span<const typename E::value_type> a, std::span<const typename E::value_type> b,
		std::span<typename E::value_type> out, Op op = {}) {
	using T = typename E::value_type;
	if (sz(a) == 0 || sz(b) == 0) return;
	assert(sz(a) >= sz(b));
	if (sz(a) == sz(b)) {
		T r{};
		for (int i = 0; i < sz(a); i++) {
			r += a[i] * b[sz(b) - 1 - i];
		}
		if (sz(out) > 0) op(out[0], r);
		return;
	}
	auto [n, cut] = detail::conv_size_for(sz(a));
	auto buf = buffer_pool<T>::get(n);
	multiply_circular<E>(a, b, buf.span(), n);
	detail::emit_middle<T>(buf.span(), cut, sz(a), sz(b),
			a[0] * b[0], a[sz(a) - 1] * b[sz(b) - 1], out, op);
}

// TODO: Let's decide whether to keep vector<> returning forms or not; this
// largely depends on whether we think these functions are a public interface or
// merely convenience for value type implementors.
template <engine E> vector<typename E::value_type> middle_product(
		std::span<const typename E::value_type> a, std::span<const typename E::value_type> b) {
	using T = typename E::value_type;
	if (sz(a) == 0 || sz(b) == 0) return {};
	assert(sz(a) >= sz(b));
	vector<T> r(size_t(sz(a) - sz(b) + 1));
	middle_product<E>(a, b, std::span<T>(r));
	return r;
}

template <engine E, typename Op = assign_op>
void middle_product(std::span<const typename E::value_type> a, transformed<E>& ta,
		std::span<const typename E::value_type> b, transformed<E>& tb,
		std::span<typename E::value_type> out, Op op = {}) {
	using T = typename E::value_type;
	if (sz(a) == 0 || sz(b) == 0) return;
	assert(sz(a) >= sz(b));
	if (sz(a) == sz(b)) {
		T r{};
		for (int i = 0; i < sz(a); i++) {
			r += a[i] * b[sz(b) - 1 - i];
		}
		if (sz(out) > 0) op(out[0], r);
		return;
	}
	auto [n, cut] = detail::conv_size_for(sz(a));
	E::extend_to(ta, n, a);
	E::extend_to(tb, n, b);
	auto buf = buffer_pool<T>::get(n);
	E::finish(E::mul(ta, tb, n), buf.span());
	detail::emit_middle<T>(buf.span(), cut, sz(a), sz(b),
			a[0] * b[0], a[sz(a) - 1] * b[sz(b) - 1], out, op);
}

template <engine E>
vector<typename E::value_type> middle_product(std::span<const typename E::value_type> a, transformed<E>& ta,
		std::span<const typename E::value_type> b, transformed<E>& tb) {
	using T = typename E::value_type;
	if (sz(a) == 0 || sz(b) == 0) return {};
	assert(sz(a) >= sz(b));
	vector<T> r(size_t(sz(a) - sz(b) + 1));
	middle_product<E>(a, ta, b, tb, std::span<T>(r));
	return r;
}

/* namespace ecnerwala::fft */ }
#line 14 "src/fft/series_core.hpp"

// ==== value types ====

namespace ecnerwala::series {

// A series is either exact (a finite series, R[x] sitting inside R[[x]]: the
// length is just the support bound) or trunc (a known prefix of an infinite
// series: the length is the precision, and products truncate to it).

// Non-owning view of power series coefficients: the span pattern (contiguous
// window + series semantics), borrowed from an owning series-like type.
template <fft::engine E, bool exact_>
struct span {
	using T = typename E::value_type;
	using engine_t = E;
	static constexpr bool exact_v = exact_;

	span() = default;
	explicit span(std::span<const T> s_) : s(s_) {}
	// exact -> trunc is implicit, trunc -> exact is explicit
	template <bool oe> requires (oe != exact_)
	explicit(oe < exact_) span(span<E, oe> o) : s(o.coeffs()) {}

	int len() const { return sz(s); }
	const T& operator[](int i) const { return s[size_t(i)]; }
	auto begin() const { return s.begin(); }
	auto end() const { return s.end(); }
	// engine primitives borrow through std::span's range constructor
	std::span<const T> coeffs() const { return s; }
	// the first n coefficients; requires n <= len().
	// Widening past len() is explicit: see with_len.
	span first(int n) const {
		assert(n <= len());
		return span(s.first(size_t(n)));
	}

private:
	std::span<const T> s;
};

// `vec` represents both exact (finite) power series (R[x]) and prefixes of infinite power series (R[[x]]), depending on the flag.
// `exact` and `trunc` are aliases.
//
// Operators here are typically permissive: they will accept combinations of unequal types and lengths.
template <fft::engine E, bool exact_>
struct vec : public std::vector<typename E::value_type> {
	using T = typename E::value_type;
	using engine_t = E;
	static constexpr bool exact_v = exact_;
	using std::vector<T>::vector;

	// a free const borrow of the coefficients: implicit
	operator span<E, exact_>() const {
		return span<E, exact_>(std::span<const T>(*this));
	}

	// exact -> trunc is implicit, trunc -> exact is explicit
	template <bool oe> requires (oe != exact_)
	explicit(oe < exact_) vec(const vec<E, oe>& p) : std::vector<T>(p) {}
	template <bool oe> requires (oe != exact_)
	explicit(oe < exact_) vec(vec<E, oe>&& p) : std::vector<T>(std::move(p)) {}

	// adopt a plain coefficient vector
	explicit vec(std::vector<T> v) : std::vector<T>(std::move(v)) {}
	// materialize an owned copy of any borrowed series, of either exactness
	explicit vec(span<E, exact_> s) : std::vector<T>(s.begin(), s.end()) {}
	explicit vec(span<E, !exact_> s) : std::vector<T>(s.begin(), s.end()) {}

	span<E, exact_> first(int n) const { return span<E, exact_>(*this).first(n); }

	int len() const {
		return int(this->size());
	}
	int degree() const requires (exact_) {
		return len() - 1;
	}
	void extend(int sz) {
		assert(sz >= len());
		this->resize(sz);
	}
	void shrink(int sz) {
		assert(sz <= len());
		this->resize(sz);
	}
	// multiply by x^n within the fixed precision window
	void shift_trunc(int n = 1) requires (!exact_) {
		assert(n >= 0 && n <= len());
		std::rotate(this->begin(), this->end()-n, this->end());
		std::fill(this->begin(), this->begin()+n, T(0));
	}
	// divide by x^n and 0-pad within the fixed precision window
	void unshift_trunc(int n = 1) requires (!exact_) {
		assert(n >= 0 && n <= len());
		std::fill(this->begin(), this->begin()+n, T(0));
		std::rotate(this->begin(), this->begin()+n, this->end());
	}

	// in-place forms require that the result's exactness/length must equal this operand's
	template <bool oe> requires (exact_ <= oe)
	vec& operator += (const vec<E, oe>& o) {
		if constexpr (exact_) { if (o.len() > len()) this->resize(o.len()); }
		else if constexpr (!oe) { if (o.len() < len()) this->resize(o.len()); }
		for (int i = 0; i < std::min(len(), o.len()); i++) {
			(*this)[i] += o[i];
		}
		return *this;
	}
	template <bool oe> requires (exact_ <= oe)
	vec& operator -= (const vec<E, oe>& o) {
		if constexpr (exact_) { if (o.len() > len()) this->resize(o.len()); }
		else if constexpr (!oe) { if (o.len() < len()) this->resize(o.len()); }
		for (int i = 0; i < std::min(len(), o.len()); i++) {
			(*this)[i] -= o[i];
		}
		return *this;
	}

	vec& operator *= (const T& n) {
		for (auto& v : *this) v *= n;
		return *this;
	}
	friend vec operator * (const vec& a, const T& n) {
		vec r(a.size());
		for (int i = 0; i < a.len(); i++) {
			r[i] = a[i] * n;
		}
		return r;
	}
	friend vec operator * (const T& n, const vec& a) {
		vec r(a.size());
		for (int i = 0; i < a.len(); i++) {
			r[i] = n * a[i];
		}
		return r;
	}

	vec& operator *= (const vec& o) {
		return *this = (*this) * o;
	}

};

template <fft::engine E> using exact = vec<E, true>;
template <fft::engine E> using trunc = vec<E, false>;

// Series-like concepts: the binary operators below are written once as constrained
// templates and dispatch on which memoized transforms an operand carries.
// A series-like type exposes its engine/exactness and its coefficients as a
// span borrow of exactly len() coefficients; cached wrappers additionally
// expose their transform caches (filling them is logically const).
template <typename S>
concept like = fft::engine<typename S::engine_t> && requires(const S& s, int i) {
	{ S::exact_v } -> std::convertible_to<bool>;
	{ s.len() } -> std::same_as<int>;
	{ s[i] } -> std::convertible_to<const typename S::engine_t::value_type&>;
	// borrows straight into the engine primitives
	{ std::span<const typename S::engine_t::value_type>(s) };
	// and into the series layer's own span, keeping the exactness tag
	requires std::convertible_to<const S&, span<typename S::engine_t, S::exact_v>>;
	// first(n): the first n coefficients, borrowed; requires n <= len().
	// The result is itself like (concepts can't self-reference).
	// Cached types keep a cache in the result only when it still serves the whole borrow.
	{ s.first(i) } -> std::convertible_to<span<typename S::engine_t, S::exact_v>>;
};
template <typename S>
concept exact_like = like<S> && S::exact_v;
template <typename S>
concept trunc_like = like<S> && !S::exact_v;

// carries one extendable transform of the whole coefficient sequence
template <typename S>
concept has_cache = like<S> && requires(const S& s) {
	{ s.cache() } -> std::same_as<fft::transformed<typename S::engine_t>&>;
};

template <fft::engine E, bool exact_>
struct maybe_cached;

// A borrowed series paired with the transform serving it: the
// normalized operand form fed to the cached fft:: entry points. Models has_cache.
template <fft::engine E, bool exact_>
struct cached_span {
	using engine_t = E;
	static constexpr bool exact_v = exact_;

	span<E, exact_> s;
	std::reference_wrapper<fft::transformed<E>> f;

	cached_span(span<E, exact_> s_, fft::transformed<E>& f_) : s(s_), f(f_) {}
	// exact -> trunc is implicit, trunc -> exact is explicit
	template <bool oe> requires (oe != exact_)
	explicit(oe < exact_) cached_span(cached_span<E, oe> o) : s(span<E, exact_>(o.s)), f(o.f) {}

	int len() const { return s.len(); }
	const typename E::value_type& operator[](int i) const { return s[i]; }
	operator std::span<const typename E::value_type>() const { return s.coeffs(); }
	operator span<E, exact_>() const { return s; }
	maybe_cached<E, exact_> first(int n) const;
	fft::transformed<E>& cache() const { return f; }
};

// carries a whole-sequence cache only sometimes, decided at runtime
template <typename S>
concept has_cache_opt = like<S> && requires(const S& s) {
	{ s.cache_opt() } -> std::same_as<std::optional<std::reference_wrapper<fft::transformed<typename S::engine_t>>>>;
};

namespace detail {
// the operand's whole cache, if it carries one
template <like S>
std::optional<std::reference_wrapper<fft::transformed<typename S::engine_t>>> cache_of(const S& s) {
	if constexpr (has_cache<S>) return s.cache();
	else if constexpr (has_cache_opt<S>) return s.cache_opt();
	else return std::nullopt;
}
/* namespace detail */ }

// A borrowed series which may carry the transform serving it: the runtime
// counterpart of cached_span in the borrow hierarchy
// prefix_cached/cached -> maybe_cached/cached_span -> span.
template <fft::engine E, bool exact_>
struct maybe_cached {
	using T = typename E::value_type;
	using engine_t = E;
	static constexpr bool exact_v = exact_;

	span<E, exact_> s;
	std::optional<std::reference_wrapper<fft::transformed<E>>> f;

	explicit maybe_cached(span<E, exact_> s_) : s(s_) {}
	maybe_cached(span<E, exact_> s_, fft::transformed<E>& f_) : s(s_), f(f_) {}
	maybe_cached(cached_span<E, exact_> c) : s(c.s), f(c.f) {}
	// borrow any like operand whole, taking along whatever cache it carries
	template <like S> requires std::same_as<typename S::engine_t, E> && (S::exact_v == exact_)
	maybe_cached(const S& o) : s(o), f(detail::cache_of(o)) {}

	int len() const { return s.len(); }
	const T& operator[](int i) const { return s[i]; }
	operator std::span<const T>() const { return s.coeffs(); }
	operator span<E, exact_>() const { return s; }
	maybe_cached first(int n) const {
		return n == len() ? *this : maybe_cached(s.first(n));
	}
	std::optional<std::reference_wrapper<fft::transformed<E>>> cache_opt() const { return f; }
};

template <fft::engine E, bool exact_>
maybe_cached<E, exact_> cached_span<E, exact_>::first(int n) const {
	return maybe_cached<E, exact_>(*this).first(n);
}

// An owned series copy at an adjusted logical length: the result of with_len.
// Carries a reference to a source cache when it still serves the copied
// coefficients (a zero tail doesn't change the transform), so the source
// must outlive the result.
template <fft::engine E>
struct resized {
	using T = typename E::value_type;
	using engine_t = E;
	static constexpr bool exact_v = false;

	trunc<E> s;
	std::optional<std::reference_wrapper<fft::transformed<E>>> f;

	int len() const { return s.len(); }
	const T& operator[](int i) const { return s[size_t(i)]; }
	operator std::span<const T>() const { return std::span<const T>(s); }
	operator span<E, false>() const { return s; }
	maybe_cached<E, false> first(int n) const {
		span<E, false> v = s;
		if (n == len() && f) return {v, f->get()};
		return maybe_cached<E, false>(v.first(n));
	}
	std::optional<std::reference_wrapper<fft::transformed<E>>> cache_opt() const { return f; }
};

// copy any operand to logical length n: extending zero-fills, shrinking truncates
template <like S>
resized<typename S::engine_t> with_len(const S& s, int n) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	auto p = s.first(std::min(n, s.len()));
	resized<E> r;
	r.s.assign(size_t(n), T{});
	std::span<const T> pc(p);
	std::copy(pc.begin(), pc.end(), r.s.begin());
	r.f = detail::cache_of(p);
	return r;
}

// carries memoized transforms of power-of-two prefixes (see prefix_cached):
// product operands truncate to a covered scale to reuse them.
// Trunc-only: an exact operand participates whole, so has_cache covers it.
template <typename S>
concept has_prefix_cache = like<S> && !S::exact_v && requires(const S& s, int n) {
	{ s.prefix_cache(n) } -> std::same_as<fft::transformed<typename S::engine_t>&>;
};

// Wrapper around vec which caches the transform of the whole series.
// Ops exploit the cache whenever the whole span participates; a trunc series'
// whole-sequence transform is still useful for middle products and repeated
// full-precision use.
template <fft::engine E, bool exact_>
struct cached {
	using T = typename E::value_type;
	using engine_t = E;
	static constexpr bool exact_v = exact_;

	cached() = default;
	// moving coefficients in or out is free: implicit on rvalues, explicit copy otherwise
	cached(vec<E, exact_>&& s_) : s(std::move(s_)) {}
	explicit cached(const vec<E, exact_>& s_) : s(s_) {}
	operator vec<E, exact_>() && { return std::move(s); }

	int len() const { return s.len(); }
	// unwrap to the owned coefficients
	const vec<E, exact_>& uncached() const { return s; }
	const T& operator[](int i) const { return s[size_t(i)]; }
	auto begin() const { return s.cbegin(); }
	auto end() const { return s.cend(); }
	operator span<E, exact_>() const { return s; }
	maybe_cached<E, exact_> first(int n) const {
		return n == len() ? maybe_cached<E, exact_>(s, f) : maybe_cached<E, exact_>(s.first(n));
	}
	// the transform of the coefficients, fed to the cached fft:: entry points alongside them
	fft::transformed<E>& cache() const { return f; }

	template <like S>
	friend bool operator==(const cached& a, const S& b) {
		span<E, S::exact_v> bs = b;
		return a.len() == bs.len() && std::equal(a.s.begin(), a.s.end(), bs.begin());
	}

private:
	vec<E, exact_> s;
	mutable fft::transformed<E> f; // memoized transform: filling it is logically const
};

template <fft::engine E> using cached_exact = cached<E, true>;
template <fft::engine E> using cached_trunc = cached<E, false>;

namespace detail {
// Normalize a whole-span operand to a cached_span: the coefficients borrowed
// together with the cache serving them (the operand's own, or tmp otherwise).
// The whole-span multiply/square/middle_product paths run entirely on this form.
template <like S>
cached_span<typename S::engine_t, S::exact_v> as_cached_span(const S& s, fft::transformed<typename S::engine_t>& tmp) {
	auto co = cache_of(s);
	return {s, co ? co->get() : tmp};
}
/* namespace detail */ }

// Newton inversion: 1/a mod x^a.len(). Generic over any engine; per doubling step
// n -> m = 2n this is 5 transforms of size m, reusing b's transform for both circular
// products; in each product the wraparound only contaminates coefficients [0, n)
// which are already known.
//
// This is correct for non-commutative rings.
// TODO: reuse/populate the operand's whole/prefix transform caches
template <trunc_like S>
trunc<typename S::engine_t> ps_inv(const S& a) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	int N = a.len();
	trunc<E> r(size_t(N), T{});
	if (N == 0) return r;
	int s = nextPow2(N);
	std::vector<T> b(size_t(s), T{});
	b[0] = inv(a[0]);
	for (int n = 1; n < N; n *= 2) {
		int m = 2 * n;
		auto ta = E::transform(a.first(std::min(N, m)), m);
		auto tb = E::transform(std::span<const T>(b).first(n), m);
		// e = a*b mod x^m; only e[n..m) is needed (and is wraparound-free).
		auto e = fft::buffer_pool<T>::get(m);
		E::finish(E::mul(ta, tb, m), e.span());
		for (int i = 0; i < n; i++) e[i] = T{};
		auto te = E::transform(std::span<const T>(e.span()), m);
		auto c = fft::buffer_pool<T>::get(m);
		// b' = 2b - b*(a*b): keep b on the left of e = a*b
		E::finish(E::mul(tb, te, m), c.span());
		for (int i = n; i < std::min(m, N); i++) b[i] = -c[i];
	}
	std::copy(b.begin(), b.begin() + N, r.begin());
	return r;
}
// TODO: operator / can be done slightly faster than ps_inv:
// we only need the n/2 terms of ps_inv(), and can do the last Newton step directly on the quotient

// Both consume whole-sequence transforms by nature (the full span always
// participates), so only whole caches apply, never prefix caches.
template <like A>
auto square(const A& a) {
	using E = typename A::engine_t;
	using T = typename E::value_type;
	fft::transformed<E> ta_;
	auto av = detail::as_cached_span(a, ta_);
	if constexpr (A::exact_v) {
		// like operator*, an exact square returns has_cache, adopting the
		// pointwise product as the result's transform when the engine supports it
		std::vector<T> coeffs;
		fft::transformed<E> f;
		fft::square_cached<E>(av, av.cache(), coeffs, f);
		cached<E, true> w(exact<E>(std::move(coeffs)));
		w.cache() = std::move(f);
		return w;
	} else {
		trunc<E> r(size_t(a.len()), T{});
		fft::square<E>(av, av.cache(), std::span<T>(r));
		return r;
	}
}

// a*b + c*d, all exact; returns has_cache, adopting the summed pointwise
// product as the result's transform when the engine supports it. Reuses each
// operand's whole cache. Requires a*b and c*d to have equal length.
template <like A, like B, like C, like D>
	requires fft::same_engine<A, B> && fft::same_engine<A, C> && fft::same_engine<A, D>
		&& A::exact_v && B::exact_v && C::exact_v && D::exact_v
cached<typename A::engine_t, true> multiply_add2(
		const A& a, const B& b, const C& c, const D& d) {
	using E = typename A::engine_t;
	using T = typename E::value_type;
	fft::transformed<E> ta_, tb_, tc_, td_;
	auto av = detail::as_cached_span(a, ta_), bv = detail::as_cached_span(b, tb_);
	auto cv = detail::as_cached_span(c, tc_), dv = detail::as_cached_span(d, td_);
	std::vector<T> coeffs;
	fft::transformed<E> f;
	fft::multiply_add2_cached<E>(
		av, av.cache(),
		bv, bv.cache(),
		cv, cv.cache(),
		dv, dv.cache(),
		coeffs, f
	);
	cached<E, true> w(exact<E>(std::move(coeffs)));
	w.cache() = std::move(f);
	return w;
}

// coefficients [b.len()-1, a.len()) of a*b; requires a.len() >= b.len() > 0.
// The kernel b participates whole, so it must be exact; the result mirrors a's kind.
template <like A, exact_like B> requires fft::same_engine<A, B>
vec<typename A::engine_t, A::exact_v> middle_product(const A& a, const B& b) {
	using E = typename A::engine_t;
	fft::transformed<E> ta_, tb_;
	auto av = detail::as_cached_span(a, ta_);
	auto bv = detail::as_cached_span(b, tb_);
	return vec<E, A::exact_v>(fft::middle_product<E>(
		av, av.cache(),
		bv, bv.cache()
	));
}

namespace detail {
template <bool ea, bool eb> int product_prec(int la, int lb) {
	if constexpr (ea && eb) return la > 0 && lb > 0 ? la + lb - 1 : 0;
	else return ea ? lb : eb ? la : std::min(la, lb);
}

// Normalize a product operand at the given precision to a borrowed series + the
// whole cache serving it: a prefix cache at scale nextPow2(prec), or the whole
// span with the operand's own cache, or a truncated span with the caller's
// throwaway cache.
// A whole cache of an over-length operand (len > prec, which pins the other,
// necessarily trunc, operand's span at exactly prec) is only worth using when
// the untruncated span doesn't grow the transform size: a 2x'd inverse
// transform costs more than the saved forward transform.
template <like S>
auto product_operand(const S& s, int prec, fft::transformed<typename S::engine_t>& tmp) {
	using E = typename S::engine_t;
	if constexpr (has_prefix_cache<S>) {
		return s.first(std::min(s.len(), nextPow2(prec)));
	} else {
		span<E, S::exact_v> v = s;
		int used = std::min(v.len(), prec);
		if (auto co = cache_of(s)) {
			if (s.len() <= prec || fft::detail::conv_size_for(s.len() + prec - 1).n
					== fft::detail::conv_size_for(2 * prec - 1).n) {
				return cached_span<E, S::exact_v>{v, co->get()};
			}
		}
		return cached_span<E, S::exact_v>{v.first(used), tmp};
	}
}
/* namespace detail */ }

template <like A, like B> requires fft::same_engine<A, B>
vec<typename A::engine_t, A::exact_v && B::exact_v> operator + (const A& a, const B& b) {
	using T = typename A::engine_t::value_type;
	int n = (A::exact_v && B::exact_v) ? std::max(a.len(), b.len())
		: A::exact_v ? b.len() : B::exact_v ? a.len() : std::min(a.len(), b.len());
	vec<typename A::engine_t, A::exact_v && B::exact_v> r(size_t(n), T(0));
	for (int i = 0; i < n; i++) {
		r[i] = (i < a.len() ? a[i] : T(0)) + (i < b.len() ? b[i] : T(0));
	}
	return r;
}
template <like A, like B> requires fft::same_engine<A, B>
vec<typename A::engine_t, A::exact_v && B::exact_v> operator - (const A& a, const B& b) {
	using T = typename A::engine_t::value_type;
	int n = (A::exact_v && B::exact_v) ? std::max(a.len(), b.len())
		: A::exact_v ? b.len() : B::exact_v ? a.len() : std::min(a.len(), b.len());
	vec<typename A::engine_t, A::exact_v && B::exact_v> r(size_t(n), T(0));
	for (int i = 0; i < n; i++) {
		r[i] = (i < a.len() ? a[i] : T(0)) - (i < b.len() ? b[i] : T(0));
	}
	return r;
}

// The single multiplication operator: each operand is normalized to a borrowed
// series + whole cache (see detail::product_operand), then multiplied once.
// An exact x exact product returns a has_cache result, going through
// fft::multiply_cached so the pointwise product is adopted as the result's
// transform whenever the engine supports it.
template <like A, like B> requires fft::same_engine<A, B>
auto operator * (const A& a, const B& b) {
	using E = typename A::engine_t;
	using T = typename E::value_type;
	constexpr bool ea = A::exact_v, eb = B::exact_v;
	int prec = detail::product_prec<ea, eb>(a.len(), b.len());
	if (prec == 0 || a.len() == 0 || b.len() == 0) {
		if constexpr (ea && eb) return cached<E, true>{};
		else return trunc<E>(size_t(prec), T(0));
	}
	fft::transformed<E> ta_, tb_;
	auto va = detail::product_operand(a, prec, ta_);
	auto vb = detail::product_operand(b, prec, tb_);
	if constexpr (ea && eb) {
		std::vector<T> coeffs;
		fft::transformed<E> f;
		auto ca = detail::as_cached_span(va, ta_), cb = detail::as_cached_span(vb, tb_);
		fft::multiply_cached<E>(
			ca, ca.cache(),
			cb, cb.cache(),
			coeffs, f
		);
		cached<E, true> w(exact<E>(std::move(coeffs)));
		w.cache() = std::move(f);
		return w;
	} else {
		trunc<E> r(size_t(prec), T(0));
		auto ca = detail::as_cached_span(va, ta_), cb = detail::as_cached_span(vb, tb_);
		fft::multiply<E>(
			ca, ca.cache(),
			cb, cb.cache(),
			std::span<T>(r)
		);
		return r;
	}
}

// Wrapper around trunc which caches transform(s[:2^k]) for all k,
// matching the doubling shape of ps_inv/exp so they can populate the caches.
// TODO: make ps_inv/exp populate these
template <fft::engine E>
struct prefix_cached {
	using T = typename E::value_type;

	using engine_t = E;
	static constexpr bool exact_v = false;

	prefix_cached() = default;
	// moving coefficients in or out is free: implicit on rvalues, explicit copy otherwise
	prefix_cached(trunc<E>&& s_) : s(std::move(s_)) {}
	explicit prefix_cached(const trunc<E>& s_) : s(s_) {}
	operator trunc<E>() && { return std::move(s); }

	int len() const { return s.len(); }
	// unwrap to the owned coefficients
	const trunc<E>& uncached() const { return s; }
	const T& operator[](int i) const { return s[size_t(i)]; }
	auto begin() const { return s.cbegin(); }
	auto end() const { return s.cend(); }
	operator span<E, false>() const { return s; }

	// extend precision: appends coefficients, keeping all covering caches valid
	void append(std::span<const T> tail) {
		s.insert(s.end(), tail.begin(), tail.end());
	}

	// The first k coefficients, with the covering prefix cache when one lines
	// up with k (k a power of two, or k == len()); uncached otherwise.
	maybe_cached<E, false> first(int k) const {
		assert(k <= len());
		span<E, false> v = s.first(k);
		int n = nextPow2(k);
		if (std::min(n, len()) == k) return {v, prefix_cache(n)};
		return maybe_cached<E, false>(v);
	}
	// the whole-sequence transform: the prefix cache covering all of len()
	fft::transformed<E>& cache() const { return prefix_cache(nextPow2(len())); }
	// cache over the prefix of length min(n, len()); n a power of two
	fft::transformed<E>& prefix_cache(int n) const {
		assert(n > 0 && !(n & (n-1)));
		int k = __builtin_ctz(unsigned(n));
		if (k >= sz(caches)) caches.resize(size_t(k) + 1);
		auto& c = caches[k];
		int e = std::min(n, len());
		if (c.len != e) {
			c.t = E::transform(s.first(e), 2 * n);
			c.len = e;
		}
		return c.t;
	}

private:
	trunc<E> s;
	// memoized transforms: logically const; len tracks how much of s each covers
	struct entry { fft::transformed<E> t; int len = 0; };
	mutable std::vector<entry> caches;
};

/* namespace ecnerwala::series */ }
#line 11 "src/fft/series.hpp"

// ==== analytic ops ====
// Free functions over series-like operands; each borrows the operand's span
// and writes a fresh result.
// TODO: reuse/populate the operands' whole/prefix transform caches

namespace ecnerwala::series {

template <like S>
vec<typename S::engine_t, S::exact_v> stretch(const S& a_, int n) {
	using E = typename S::engine_t;
	span<E, S::exact_v> a = a_;
	vec<E, S::exact_v> r(size_t(a.len()));
	for (int i = 0; i*n < a.len(); i++) {
		r[i*n] = a[i];
	}
	return r;
}
template <like S>
vec<typename S::engine_t, S::exact_v> deriv_shift(const S& a_) {
	using E = typename S::engine_t;
	span<E, S::exact_v> a = a_;
	vec<E, S::exact_v> r(a.begin(), a.end());
	for (int i = 0; i < r.len(); i++) {
		r[i] *= i;
	}
	return r;
}
template <like S>
vec<typename S::engine_t, S::exact_v> integ_shift(const S& a_) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, S::exact_v> a = a_;
	assert(a[0] == 0);
	vec<E, S::exact_v> r(a.begin(), a.end());
	T f = 1;
	for (int i = 1; i < r.len(); i++) {
		r[i] *= f;
		f *= i;
	}
	f = inv(f);
	for (int i = r.len() - 1; i > 0; i--) {
		r[i] *= f;
		f *= i;
	}
	return r;
}
template <like S>
vec<typename S::engine_t, S::exact_v> integ_shift_offset(const S& a_, int offset) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, S::exact_v> a = a_;
	vec<E, S::exact_v> r(a.begin(), a.end());
	T f = 1;
	for (int i = 0; i < r.len(); i++) {
		r[i] *= f;
		f *= i + offset;
	}
	assert(f != 0);
	f = inv(f);
	for (int i = r.len() - 1; i >= 0; i--) {
		r[i] *= f;
		f *= i + offset;
	}
	return r;
}
template <trunc_like S>
trunc<typename S::engine_t> deriv_shift_log(const S& a) {
	return deriv_shift(a) * ps_inv(a);
}
template <trunc_like S>
trunc<typename S::engine_t> ps_log(const S& a) {
	assert(a[0] == 1);
	return integ_shift(deriv_shift_log(a));
}
template <trunc_like S>
trunc<typename S::engine_t> ps_exp(const S& a_) {
	// See https://mathexp.eu/bostan/publications/BoSc09a.pdf for details
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, false> a = a_;
	assert(a.len() >= 1);
	assert(a[0] == 0);
	trunc<E> r(1, T(1)); r.reserve(size_t(a.len()));
	trunc<E> invR(1, T(1)); invR.reserve(size_t(a.len()));
	while (r.len() < a.len()) {
		int o_sz = r.len();
		int n_sz = std::min(o_sz * 2, a.len());
		trunc<E> t = deriv_shift(trunc<E>(a.begin(), a.begin() + o_sz));
		fft::multiply_circular<E>(std::span<const T>(t), std::span<const T>(r).first(o_sz), std::span<T>(t), o_sz);
		t = deriv_shift(r) - t;
		t *= invR;
		t.resize(size_t(n_sz - o_sz));
		trunc<E> v(a.begin() + o_sz, a.begin() + n_sz);
		v -= integ_shift_offset(t, o_sz);
		v *= r;
		r.resize(size_t(n_sz));
		std::copy(v.begin(), v.end(), r.begin() + o_sz);
		if (r.len() < a.len()) {
			// double invR via a Newton step
			assert(r.len() == 2 * invR.len());
			int n = invR.len();
			int nn = r.len();
			trunc<E> tmp(size_t(4) * n);
			fft::square<E>(std::span<const T>(invR).first(n), std::span<T>(tmp));
			fft::multiply<E>(std::span<const T>(tmp).first(nn), std::span<const T>(r).first(nn), std::span<T>(tmp));
			invR.resize(size_t(nn));
			for (int i = n; i < nn; i++) invR[i] = -tmp[i];
		}
	}
	return r;
}
template <trunc_like S>
trunc<typename S::engine_t> ps_pow_monic(const S& a_, typename S::engine_t::value_type k) {
	using E = typename S::engine_t;
	span<E, false> a = a_;
	if (a.len() == 0) return {};
	assert(a[0] == 1);
	trunc<E> l = ps_log(a_);
	l *= k;
	return ps_exp(l);
}
template <trunc_like S>
trunc<typename S::engine_t> ps_pow(const S& a_, int64_t k) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, false> a = a_;
	assert(k >= 0);
	if (k == 0) {
		trunc<E> r(size_t(a.len()), T(0));
		if (r.len() > 0) r[0] = T(1);
		return r;
	}

	int st = 0;
	while (st < a.len() && a[st] == 0) st++;

	if (st > 0 && k > (a.len() - 1) / st) {
		return trunc<E>(size_t(a.len()), T(0));
	}

	trunc<E> r(a.begin() + st, a.end() - (st * (k-1)));
	T leading_coeff = r[0];
	r *= inv(leading_coeff);
	r = ps_pow_monic(r, T(k));
	r *= power(leading_coeff, k);
	r.insert(r.begin(), size_t(st * k), T(0));
	assert(r.len() == a.len());
	return r;
}

template <trunc_like S>
trunc<typename S::engine_t> to_newton_sums(const S& a, int deg) {
	auto r = deriv_shift_log(a);
	r[0] = deg;
	for (int i = 1; i < r.len(); i++) r[i] = -r[i];
	return r;
}
template <trunc_like S>
trunc<typename S::engine_t> from_newton_sums(const S& s_, int deg) {
	using E = typename S::engine_t;
	span<E, false> s = s_;
	assert(s[0] == deg);
	trunc<E> r(s.begin(), s.end());
	r[0] = 0;
	for (int i = 1; i < r.len(); i++) r[i] = -r[i];
	return ps_exp(integ_shift(std::move(r)));
}

// Calculates prod 1/(1-x^i)^{a[i]}
template <trunc_like S>
trunc<typename S::engine_t> euler_transform(const S& a) {
	using E = typename S::engine_t;
	trunc<E> r = deriv_shift(a);
	std::vector<bool> is_prime(size_t(r.len()), true);
	for (int p = 2; p < r.len(); p++) {
		if (!is_prime[p]) continue;
		for (int i = 1; i*p < r.len(); i++) {
			r[i*p] += r[i];
			is_prime[i*p] = false;
		}
	}
	return ps_exp(integ_shift(r));
}
template <trunc_like S>
trunc<typename S::engine_t> inverse_euler_transform(const S& a) {
	using E = typename S::engine_t;
	trunc<E> r = deriv_shift(ps_log(a));
	std::vector<bool> is_prime(size_t(r.len()), true);
	for (int p = 2; p < r.len(); p++) {
		if (!is_prime[p]) continue;
		for (int i = (r.len()-1)/p; i >= 1; i--) {
			r[i*p] -= r[i];
			is_prime[i*p] = false;
		}
	}
	return integ_shift(r);
}

// Helper packed bivariate buffer for Kinoshita-Li composition (arXiv:2404.05177).
//
// The motivation is performing Bostan-Mori (Graeffe root-squaring) to compute
// something like [x^n] P / Q_0(x, y) with deg_y(Q_0) = 1 and deg_x(Q_0) = n.
//
// In each step, we want to compute Q_{i+1}(x^2, y) = Q_i(x, y) * Q_i(-x, y).
// This doubles the degree of y and also lets us truncate x at half the previous
// degree, leaving the total size invariant.
//
// We will store Q as a packed buffer with x as the inner dimension to facilitate easy Q(-x) substitution.
// The inner span will be 2*deg(x), and the outer span will be 2*deg(y).
// As we advance, we will also return the cached transform of Q_i(-x, y) for the caller to use in the numerator.
template <fft::engine E> struct packed_bivariate {
	using T = typename E::value_type;
	int L, l;
	std::vector<T> c;

	// Q_0 = 1 - y g(x), deg g < n <= 2^L
	packed_bivariate(int L_, std::span<const T> g) : L(L_), l(0), c(size_t(4) << L) {
		c[0] = T(1);
		for (int i = 0; i < sz(g); i++) c[(2 << L) + i] = -g[i];
	}

	fft::transformed<E> advance() {
		int B = 4 << L;
		auto tq = E::transform(std::span<const T>(c), B);
		auto tn = E::negate_arg(tq, B);
		E::finish(
			E::downsample(E::mul(tq, tn, B), B/2, false),
			std::span<T>(c).first(B/2)
		);
		l++;
		// undo the circular wraparound using monicity in y
		for (int i = 0; i < (2 << (L - l)); i++) {
			c[(2 << L) + i] = c[i];
			c[i] = T(0);
		}
		c[2 << L] -= T(1);
		c[0] = T(1);
		// zero x coefficients beyond the level's truncation mod x^(2^(L-l))
		std::fill(c.begin() + (2 << L) + (1 << (L - l)), c.end(), T(0));
		for (int i = 0; i < (2 << L); i += 2 << (L - l)) {
			for (int j = 0; j < (1 << (L - l)); j++) {
				c[i + (1 << (L - l)) + j] = T(0);
			}
		}
		return tn;
	}
};

// Calculates f(g(x)) mod x^n where deg(g) == n
template <trunc_like SF, trunc_like SG> requires fft::same_engine<SF, SG>
trunc<typename SF::engine_t> ps_compose(const SF& f_, const SG& g_) {
	using E = typename SF::engine_t;
	using T = typename E::value_type;
	span<E, false> f = f_;
	span<E, false> g = g_;
	if (g.len() == 0) return {};

	int m = f.len();
	int n = g.len();

	// https://arxiv.org/pdf/2404.05177
	// Consider P(y) = f(1/y) has terms from y^{-(m-1)}...y^0 (Laurent series)
	// We want [y^0] P(y) / (1 - y g(x))
	// Let Q_0 = 1 - yg(x)
	// Q_{i+1}(x^2, y) = Q_i(x, y) * Q_i(-x, y) mod x^{ceil(n / 2^i)}
	// deg_y(Q_i) = 2^i, deg_x(Q_i) = ceil(n / 2^i) - 1
	//
	// [y^0] P(y) / Q_l(x^2^l, y) * Q_{l-1}(-x^2^{l-1}, y) * Q_{l-2}(-x^2^{l-2}, y) * ... * Q_0(-x, y)
	// The total y deg of Q_{k-1} ... Q_0 is 2^k-1
	int L = __builtin_ctz(unsigned(nextPow2(n)));
	int B = 4 << L;
	packed_bivariate<E> Q(L, g.coeffs());
	// tneg[l] is the transform of Q_l(-x, y), reused by the pushdown pass below
	std::vector<fft::transformed<E>> tneg;
	tneg.reserve(L);
	for (int l = 1; l <= L; l++) tneg.push_back(Q.advance());
	trunc<E> P;
	{
		P = trunc<E>(f.begin(), f.end());
		std::reverse(P.begin(), P.end());
		trunc<E> QL((1 << L) + 1);
		for (int i = 0; i <= (1 << L); i++) {
			QL[i] = Q.c[2 * i];
		}
		QL.resize(size_t(m), T(0));
		P *= ps_inv(QL);
		std::reverse(P.begin(), P.end());
		P.resize(size_t(1) << L, T(0));
		std::reverse(P.begin(), P.end());
		P.resize(size_t(B), T(0));
		for (int i = (1 << L) - 1; i > 0; i--) {
			P[2*i] = P[i];
			P[i] = T(0);
		}
	}
	for (int l = L-1; l >= 0; l--) {
		// Spread it out, clear the high terms
		for (int i = (2 << L) - 1; i > 0; i--) {
			T v = P[i];
			P[2*i] = ((2*i) & (1 << (L-l))) ? T(0) : v;
			P[i] = T(0);
		}
		auto tp = E::transform(std::span<const T>(P), B);
		E::finish(E::mul(tneg[l], tp, B), std::span<T>(P));
		for (int i = 0; i < (2 << L); i++) {
			P[i] = P[(2 << L) + i];
			P[(2 << L) + i] = T(0);
		}
	}
	return trunc<E>(P.begin(), P.begin() + n);
}

// [x^k] p(x)/q(x) (Bostan-Mori) for an exact rational function.
template <exact_like P, exact_like Q> requires fft::same_engine<P, Q>
P::engine_t::value_type kth_term_of_rational_function(
	const P& p,
	const Q& q,
	uint64_t k
) {
	using E = P::engine_t;
	using T = E::value_type;

	assert(q.len() > 0 && q[0] != T(0));
	// Check this here so we avoid accessing p[0]
	if (p.len() == 0) return T(0);

	// Size up in a pretty conservative way
	int d = std::max(p.len() + 1, q.len());
	assert(d >= 2);

	int n = nextPow2((d-1) + d - 1); // >= d

	// Seed the loop transforms from any whole caches; the buffers below hold the
	// current p, q (zero-padded, which extend_to tolerates).
	fft::transformed<E> tq, tp;
	if (auto cq = detail::cache_of(q)) { E::extend_to(cq->get(), n, q); tq = cq->get(); }
	if (auto cp = detail::cache_of(p)) { E::extend_to(cp->get(), n, p); tp = cp->get(); }

	std::vector<T> p_buf(d-1, T(0));
	std::ranges::copy(std::span<const T>(p), p_buf.begin());
	std::vector<T> q_buf(d, T(0));
	std::ranges::copy(std::span<const T>(q), q_buf.begin());

	while (k > 0) {
		E::extend_to(tq, n, q_buf);
		auto tnq = E::negate_arg(tq, n);
		E::extend_to(tp, n, p_buf);

		// P <- downsample(P(x) * Q(-x))
		auto ntp = E::downsample(E::mul(tp, tnq, n), n/2, bool(k & 1));
		assert(ntp.size() == n/2);
		if constexpr (std::same_as<typename E::product, typename E::transformed>) {
			tp = ntp;
		} else {
			tp = {};
		}
		E::finish(std::move(ntp), std::span(p_buf));
		k >>= 1;

		// Save the last iteration if we're done
		if (!k) {
			// HACK: fix the constant coefficient of q only
			q_buf[0] *= q_buf[0];
			break;
		}

		// Q <- downsample(Q(x) * Q(-x))
		auto ntq = E::downsample(E::mul(tq, tnq, n), n/2, false);
		assert(ntq.size() == n/2);
		if constexpr (std::same_as<typename E::product, typename E::transformed>) {
			tq = ntq;
		} else {
			tq = {};
		}
		if (n/2 == d-1) {
			// Fix the wraparound
			T v0 = q_buf[0] * q_buf[0];
			E::finish(std::move(ntq), std::span(q_buf).first(d-1));
			q_buf[d-1] = std::exchange(q_buf[0], v0) - v0;
		} else {
			E::finish(std::move(ntq), std::span(q_buf));
		}
	}
	return p_buf[0] * inv(q_buf[0]);
}

// Find the kth term of linearly recurrent sequence S with char poly Q and len(S) >= len(Q)-1
template <trunc_like S, exact_like Q> requires fft::same_engine<S, Q>
S::engine_t::value_type kth_term_of_linear_recurrence(
	const S& s,
	const Q& q,
	uint64_t k
) {
	using E = S::engine_t;
	using T = E::value_type;

	assert(q.len() > 0 && q[0] != T(0));
	assert(s.len() >= q.len()-1);

	// Don't even bother with P so we don't have to do truncation checks
	// TODO: Could use generic multiply for this whole part?
	fft::transformed<E> tq;
	auto q_cached = detail::as_cached_span(q, tq);

	// Compute the prefix and then hard-cast it to exact
	span<E, false> sv = s;
	auto p = exact<E>(sv.first(q.len()-1) * q_cached);
	return kth_term_of_rational_function(p, q_cached, k);
}

/* namespace ecnerwala::series */ }

// clang-format off
// @formatter:off
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wpragmas"
#pragma GCC diagnostic ignored "-Wunknown-warning-option"
#pragma GCC diagnostic ignored "-Wmisleading-indentation"
#pragma GCC diagnostic ignored "-Wmultistatement-macros"
#include <bits/stdc++.h>
#include <cassert>
// src/fft/common.hpp
namespace ecnerwala{
template<class T>int sz(T&&arg){using std::size;return int(size(std::forward<T>(arg)));}
inline int nextPow2(int s){return 1<<(s>1?32-__builtin_clz(s-1):0);}
namespace fft{
using std::swap;
using std::vector;
using std::min;
using std::max;
#ifndef ECNERWALA_FFT_POOL_STORAGE
#define ECNERWALA_FFT_POOL_STORAGE
#endif
template<typename T>struct buffer_pool{
static inline ECNERWALA_FFT_POOL_STORAGE std::vector<std::vector<T>>free_list;
struct handle{
std::vector<T>v;
explicit handle(int n){
if(!free_list.empty()){
v=std::move(free_list.back());
free_list.pop_back();
}
v.assign(n,T());
}
handle(const handle&)=delete;
handle&operator=(const handle&)=delete;
handle(handle&&o)noexcept:v(std::move(o.v)){}
~handle(){if(v.capacity())free_list.push_back(std::move(v));}
T&operator[](int i){return v[i];}
operator std::span<T>(){return std::span<T>(v);}
std::span<T>span(){return std::span<T>(v);}
};
static handle get(int n){return handle(n);}
};
}
}
// src/fft/engine.hpp
namespace ecnerwala::fft{
struct assign_op{template<typename T>void operator()(T&d,T v)const{d=v;}};
struct add_op{template<typename T>void operator()(T&d,T v)const{d+=v;}};
struct sub_op{template<typename T>void operator()(T&d,T v)const{d-=v;}};
struct add_twice_op{template<typename T>void operator()(T&d,T v)const{d+=v+v;}};
template<typename E>
concept engine=requires(
std::span<const typename E::value_type>in,
std::span<typename E::value_type>out,
typename E::transformed&t,
const typename E::transformed&ct,
typename E::product&p,
const typename E::product&cp,
int n
){
typename E::value_type;
{E::transform(in,n)}->std::same_as<typename E::transformed>;
{ct.size()}->std::same_as<int>;
E::extend_to(t,n,in);
{E::downsample(ct,n,false)}->std::same_as<typename E::transformed>;
{E::downsample(cp,n,false)}->std::same_as<typename E::product>;
{E::negate_arg(ct,n)}->std::same_as<typename E::transformed>;
{E::mul(ct,ct,n)}->std::same_as<typename E::product>;
{E::sq(ct,n)}->std::same_as<typename E::product>;
{E::mul2(ct,ct,ct,ct,n)}->std::same_as<typename E::template product_t<2*E::unit_scale>>;
E::finish(std::move(p),out);
E::finish(std::move(p),out,add_op{});
E::finish(E::add(std::move(p),std::move(p)),out);
{E::add(E::transform(in,n),ct)}->std::same_as<typename E::template transformed_t<2*E::unit_scale>>;
{E::add(std::move(p),std::move(p))}->std::same_as<typename E::template product_t<2*E::unit_scale>>;
requires std::same_as<std::remove_cvref_t<decltype(E::commutative)>,bool>;
requires std::same_as<std::remove_cvref_t<decltype(E::unit_scale)>,int>;
};
template<typename A,typename B>
concept same_engine=std::same_as<typename A::engine_t,typename B::engine_t>;
template<engine E>using transformed=typename E::transformed;
}
// src/fft/multiply.hpp
namespace ecnerwala::fft{
template<engine E,typename Op=assign_op>
void multiply_circular(std::span<const typename E::value_type>a,std::span<const typename E::value_type>b,
std::span<typename E::value_type>out,int n,Op op={}){
assert(!(n&(n-1)));
auto ta=E::transform(a,n);
auto tb=E::transform(b,n);
E::finish(E::mul(ta,tb,n),out,op);
}
template<engine E,typename Op=assign_op>
void square_circular(std::span<const typename E::value_type>a,std::span<typename E::value_type>out,int n,Op op={}){
assert(!(n&(n-1)));
auto ta=E::transform(a,n);
E::finish(E::sq(ta,n),out,op);
}
namespace detail{
struct conv_size{int n;bool cut;};
inline conv_size conv_size_for(int s){
int n=nextPow2(s);
bool cut=(n==2*(s-1));
return{cut?n/2:n,cut};
}
template<typename T,typename Op>
void emit_linear(std::span<T>buf,int n,int s,bool cut,T c0,std::span<T>out,Op op){
T cn{};
if(cut){
cn=buf[0]-c0;
buf[0]=c0;
}
int lim=min(sz(out),min(s,n));
for(int i=0;i<lim;i++)op(out[i],buf[i]);
if(cut&&sz(out)>=s)op(out[s-1],cn);
}
template<typename T,typename Op>
struct cut_op{
Op op;
T*out0;
T c0;
T&cn;
void operator()(T&x,T v)const{
if(&x==out0){cn=v-c0;v=c0;}
op(x,v);
}
};
template<engine E,typename P,typename Op=assign_op>
void finish_linear(
P&&p,int n,int s,bool cut,
typename E::value_type c0,std::span<typename E::value_type>out,Op op={}
){
using T=typename E::value_type;
if(sz(out)==0)return;
int lim=min(sz(out),min(s,n));
if(!cut){
E::finish(std::move(p),out.subspan(0,lim),op);
}else{
T cn{};
E::finish(std::move(p),out.subspan(0,lim),cut_op<T,Op>{op,&out[0],c0,cn});
if(sz(out)>=s)op(out[s-1],cn);
}
}
}
template<engine E,typename Op=assign_op>
void multiply(std::span<const typename E::value_type>a,std::span<const typename E::value_type>b,
std::span<typename E::value_type>out,Op op={}){
using T=typename E::value_type;
if(sz(a)==0||sz(b)==0)return;
int s=sz(a)+sz(b)-1;
auto[n,cut]=detail::conv_size_for(s);
T c0=a[0]*b[0];
auto buf=buffer_pool<T>::get(n);
multiply_circular<E>(a,b,buf.span(),n);
detail::emit_linear<T>(buf.span(),n,s,cut,c0,out,op);
}
template<engine E,typename Op=assign_op>
void multiply(std::span<const typename E::value_type>a,transformed<E>&ta,
std::span<const typename E::value_type>b,transformed<E>&tb,
std::span<typename E::value_type>out,Op op={}){
using T=typename E::value_type;
if(sz(a)==0||sz(b)==0)return;
int s=sz(a)+sz(b)-1;
auto[n,cut]=detail::conv_size_for(s);
T c0=a[0]*b[0];
E::extend_to(ta,n,a);
E::extend_to(tb,n,b);
detail::finish_linear<E>(E::mul(ta,tb,n),n,s,cut,c0,out,op);
}
template<engine E,typename Op=assign_op>
void multiply_add2(std::span<const typename E::value_type>a1,transformed<E>&ta1,
std::span<const typename E::value_type>b1,transformed<E>&tb1,
std::span<const typename E::value_type>a2,transformed<E>&ta2,
std::span<const typename E::value_type>b2,transformed<E>&tb2,
std::span<typename E::value_type>out,Op op={}){
using T=typename E::value_type;
assert(sz(a1)>0&&sz(b1)>0&&sz(a2)>0&&sz(b2)>0);
int s=sz(a1)+sz(b1)-1;
assert(sz(a2)+sz(b2)-1==s);
auto[n,cut]=detail::conv_size_for(s);
T c0=a1[0]*b1[0]+a2[0]*b2[0];
E::extend_to(ta1,n,a1);E::extend_to(tb1,n,b1);
E::extend_to(ta2,n,a2);E::extend_to(tb2,n,b2);
detail::finish_linear<E>(E::mul2(ta1,tb1,ta2,tb2,n),n,s,cut,c0,out,op);
}
template<engine E>
void multiply_add2_cached(
std::span<const typename E::value_type>a1,transformed<E>&ta1,
std::span<const typename E::value_type>b1,transformed<E>&tb1,
std::span<const typename E::value_type>a2,transformed<E>&ta2,
std::span<const typename E::value_type>b2,transformed<E>&tb2,
std::vector<typename E::value_type>&coeffs,transformed<E>&t){
using T=typename E::value_type;
assert(sz(a1)>0&&sz(b1)>0&&sz(a2)>0&&sz(b2)>0);
int s=sz(a1)+sz(b1)-1;
assert(sz(a2)+sz(b2)-1==s);
coeffs.assign(size_t(s),T{});
t=transformed<E>{};
if constexpr(std::same_as<typename E::product,transformed<E>>){
auto[n,cut]=detail::conv_size_for(s);
T c0=a1[0]*b1[0]+a2[0]*b2[0];
E::extend_to(ta1,n,a1);E::extend_to(tb1,n,b1);
E::extend_to(ta2,n,a2);E::extend_to(tb2,n,b2);
auto p=E::mul2(ta1,tb1,ta2,tb2,n);
auto tp=p;
detail::finish_linear<E>(std::move(p),n,s,cut,c0,std::span<T>(coeffs));
t=std::move(tp);
}else{
multiply_add2<E>(a1,ta1,b1,tb1,a2,ta2,b2,tb2,std::span<T>(coeffs));
}
}
template<engine E>
void multiply_cached(std::span<const typename E::value_type>a,transformed<E>&ta,
std::span<const typename E::value_type>b,transformed<E>&tb,
std::vector<typename E::value_type>&coeffs,transformed<E>&t){
using T=typename E::value_type;
coeffs.assign(size_t(sz(a)&&sz(b)?sz(a)+sz(b)-1:0),T{});
t=transformed<E>{};
if(coeffs.empty())return;
int s=sz(coeffs);
if constexpr(std::same_as<typename E::product,transformed<E>>){
auto[n,cut]=detail::conv_size_for(s);
T c0=a[0]*b[0];
E::extend_to(ta,n,a);
E::extend_to(tb,n,b);
auto p=E::mul(ta,tb,n);
auto tp=p;
detail::finish_linear<E>(std::move(p),n,s,cut,c0,std::span<T>(coeffs));
t=std::move(tp);
}else{
multiply<E>(a,ta,b,tb,std::span<T>(coeffs));
}
}
template<engine E,typename Op=assign_op>
void square(std::span<const typename E::value_type>a,std::span<typename E::value_type>out,Op op={}){
using T=typename E::value_type;
if(sz(a)==0)return;
int s=2*sz(a)-1;
auto[n,cut]=detail::conv_size_for(s);
T c0=a[0]*a[0];
auto buf=buffer_pool<T>::get(n);
square_circular<E>(a,buf.span(),n);
detail::emit_linear<T>(buf.span(),n,s,cut,c0,out,op);
}
template<engine E,typename Op=assign_op>
void square(std::span<const typename E::value_type>a,transformed<E>&ta,
std::span<typename E::value_type>out,Op op={}){
using T=typename E::value_type;
if(sz(a)==0)return;
int s=2*sz(a)-1;
auto[n,cut]=detail::conv_size_for(s);
T c0=a[0]*a[0];
E::extend_to(ta,n,a);
detail::finish_linear<E>(E::sq(ta,n),n,s,cut,c0,out,op);
}
template<engine E>
void square_cached(std::span<const typename E::value_type>a,transformed<E>&ta,
std::vector<typename E::value_type>&coeffs,transformed<E>&t){
using T=typename E::value_type;
coeffs.assign(size_t(sz(a)?2*sz(a)-1:0),T{});
t=transformed<E>{};
if(coeffs.empty())return;
int s=sz(coeffs);
if constexpr(std::same_as<typename E::product,transformed<E>>){
auto[n,cut]=detail::conv_size_for(s);
T c0=a[0]*a[0];
E::extend_to(ta,n,a);
auto p=E::sq(ta,n);
auto tp=p;
detail::finish_linear<E>(std::move(p),n,s,cut,c0,std::span<T>(coeffs));
t=std::move(tp);
}else{
square<E>(a,ta,std::span<T>(coeffs));
}
}
template<engine E>vector<typename E::value_type>multiply(
const vector<typename E::value_type>&a,const vector<typename E::value_type>&b){
using T=typename E::value_type;
if(sz(a)==0||sz(b)==0)return{};
vector<T>r(sz(a)+sz(b)-1);
multiply<E>(std::span<const T>(a),std::span<const T>(b),std::span<T>(r));
return r;
}
template<engine E>vector<typename E::value_type>square(const vector<typename E::value_type>&a){
using T=typename E::value_type;
if(sz(a)==0)return{};
vector<T>r(2*sz(a)-1);
square<E>(std::span<const T>(a),std::span<T>(r));
return r;
}
namespace detail{
template<typename T,typename Op>
void emit_middle(std::span<T>buf,bool cut,int la,int lb,T c0,T ctop,std::span<T>out,Op op){
int m=la-lb+1;
T cn{};
if(cut){
cn=buf[0]-c0;
buf[lb-1]-=ctop;
}
int lim=min(sz(out),cut?m-1:m);
for(int t=0;t<lim;t++)op(out[t],buf[lb-1+t]);
if(cut&&sz(out)>=m)op(out[m-1],cn);
}
}
template<engine E,typename Op=assign_op>
void middle_product(std::span<const typename E::value_type>a,std::span<const typename E::value_type>b,
std::span<typename E::value_type>out,Op op={}){
using T=typename E::value_type;
if(sz(a)==0||sz(b)==0)return;
assert(sz(a)>=sz(b));
if(sz(a)==sz(b)){
T r{};
for(int i=0;i<sz(a);i++){
r+=a[i]*b[sz(b)-1-i];
}
if(sz(out)>0)op(out[0],r);
return;
}
auto[n,cut]=detail::conv_size_for(sz(a));
auto buf=buffer_pool<T>::get(n);
multiply_circular<E>(a,b,buf.span(),n);
detail::emit_middle<T>(buf.span(),cut,sz(a),sz(b),
a[0]*b[0],a[sz(a)-1]*b[sz(b)-1],out,op);
}
template<engine E>vector<typename E::value_type>middle_product(
std::span<const typename E::value_type>a,std::span<const typename E::value_type>b){
using T=typename E::value_type;
if(sz(a)==0||sz(b)==0)return{};
assert(sz(a)>=sz(b));
vector<T>r(size_t(sz(a)-sz(b)+1));
middle_product<E>(a,b,std::span<T>(r));
return r;
}
template<engine E,typename Op=assign_op>
void middle_product(std::span<const typename E::value_type>a,transformed<E>&ta,
std::span<const typename E::value_type>b,transformed<E>&tb,
std::span<typename E::value_type>out,Op op={}){
using T=typename E::value_type;
if(sz(a)==0||sz(b)==0)return;
assert(sz(a)>=sz(b));
if(sz(a)==sz(b)){
T r{};
for(int i=0;i<sz(a);i++){
r+=a[i]*b[sz(b)-1-i];
}
if(sz(out)>0)op(out[0],r);
return;
}
auto[n,cut]=detail::conv_size_for(sz(a));
E::extend_to(ta,n,a);
E::extend_to(tb,n,b);
auto buf=buffer_pool<T>::get(n);
E::finish(E::mul(ta,tb,n),buf.span());
detail::emit_middle<T>(buf.span(),cut,sz(a),sz(b),
a[0]*b[0],a[sz(a)-1]*b[sz(b)-1],out,op);
}
template<engine E>
vector<typename E::value_type>middle_product(std::span<const typename E::value_type>a,transformed<E>&ta,
std::span<const typename E::value_type>b,transformed<E>&tb){
using T=typename E::value_type;
if(sz(a)==0||sz(b)==0)return{};
assert(sz(a)>=sz(b));
vector<T>r(size_t(sz(a)-sz(b)+1));
middle_product<E>(a,ta,b,tb,std::span<T>(r));
return r;
}
}
// src/fft/series_core.hpp
namespace ecnerwala::series{
template<fft::engine E,bool exact_>
struct span{
using T=typename E::value_type;
using engine_t=E;
static constexpr bool exact_v=exact_;
span()=default;
explicit span(std::span<const T>s_):s(s_){}
template<bool oe>requires(oe!=exact_)
explicit(oe<exact_)span(span<E,oe>o):s(o.coeffs()){}
int len()const{return sz(s);}
const T&operator[](int i)const{return s[size_t(i)];}
auto begin()const{return s.begin();}
auto end()const{return s.end();}
std::span<const T>coeffs()const{return s;}
span first(int n)const{
assert(n<=len());
return span(s.first(size_t(n)));
}
private:
std::span<const T>s;
};
template<fft::engine E,bool exact_>
struct vec:public std::vector<typename E::value_type>{
using T=typename E::value_type;
using engine_t=E;
static constexpr bool exact_v=exact_;
using std::vector<T>::vector;
operator span<E,exact_>()const{
return span<E,exact_>(std::span<const T>(*this));
}
template<bool oe>requires(oe!=exact_)
explicit(oe<exact_)vec(const vec<E,oe>&p):std::vector<T>(p){}
template<bool oe>requires(oe!=exact_)
explicit(oe<exact_)vec(vec<E,oe>&&p):std::vector<T>(std::move(p)){}
explicit vec(std::vector<T>v):std::vector<T>(std::move(v)){}
explicit vec(span<E,exact_>s):std::vector<T>(s.begin(),s.end()){}
explicit vec(span<E,!exact_>s):std::vector<T>(s.begin(),s.end()){}
span<E,exact_>first(int n)const{return span<E,exact_>(*this).first(n);}
int len()const{
return int(this->size());
}
int degree()const requires(exact_){
return len()-1;
}
void extend(int sz){
assert(sz>=len());
this->resize(sz);
}
void shrink(int sz){
assert(sz<=len());
this->resize(sz);
}
void shift_trunc(int n=1)requires(!exact_){
assert(n>=0&&n<=len());
std::rotate(this->begin(),this->end()-n,this->end());
std::fill(this->begin(),this->begin()+n,T(0));
}
void unshift_trunc(int n=1)requires(!exact_){
assert(n>=0&&n<=len());
std::fill(this->begin(),this->begin()+n,T(0));
std::rotate(this->begin(),this->begin()+n,this->end());
}
template<bool oe>requires(exact_<=oe)
vec&operator+=(const vec<E,oe>&o){
if constexpr(exact_){if(o.len()>len())this->resize(o.len());}
else if constexpr(!oe){if(o.len()<len())this->resize(o.len());}
for(int i=0;i<std::min(len(),o.len());i++){
(*this)[i]+=o[i];
}
return*this;
}
template<bool oe>requires(exact_<=oe)
vec&operator-=(const vec<E,oe>&o){
if constexpr(exact_){if(o.len()>len())this->resize(o.len());}
else if constexpr(!oe){if(o.len()<len())this->resize(o.len());}
for(int i=0;i<std::min(len(),o.len());i++){
(*this)[i]-=o[i];
}
return*this;
}
vec&operator*=(const T&n){
for(auto&v:*this)v*=n;
return*this;
}
friend vec operator*(const vec&a,const T&n){
vec r(a.size());
for(int i=0;i<a.len();i++){
r[i]=a[i]*n;
}
return r;
}
friend vec operator*(const T&n,const vec&a){
vec r(a.size());
for(int i=0;i<a.len();i++){
r[i]=n*a[i];
}
return r;
}
vec&operator*=(const vec&o){
return*this=(*this)*o;
}
};
template<fft::engine E>using exact=vec<E,true>;
template<fft::engine E>using trunc=vec<E,false>;
template<typename S>
concept like=fft::engine<typename S::engine_t>&&requires(const S&s,int i){
{S::exact_v}->std::convertible_to<bool>;
{s.len()}->std::same_as<int>;
{s[i]}->std::convertible_to<const typename S::engine_t::value_type&>;
{std::span<const typename S::engine_t::value_type>(s)};
requires std::convertible_to<const S&,span<typename S::engine_t,S::exact_v>>;
{s.first(i)}->std::convertible_to<span<typename S::engine_t,S::exact_v>>;
};
template<typename S>
concept exact_like=like<S>&&S::exact_v;
template<typename S>
concept trunc_like=like<S>&&!S::exact_v;
template<typename S>
concept has_cache=like<S>&&requires(const S&s){
{s.cache()}->std::same_as<fft::transformed<typename S::engine_t>&>;
};
template<fft::engine E,bool exact_>
struct maybe_cached;
template<fft::engine E,bool exact_>
struct cached_span{
using engine_t=E;
static constexpr bool exact_v=exact_;
span<E,exact_>s;
std::reference_wrapper<fft::transformed<E>>f;
cached_span(span<E,exact_>s_,fft::transformed<E>&f_):s(s_),f(f_){}
template<bool oe>requires(oe!=exact_)
explicit(oe<exact_)cached_span(cached_span<E,oe>o):s(span<E,exact_>(o.s)),f(o.f){}
int len()const{return s.len();}
const typename E::value_type&operator[](int i)const{return s[i];}
operator std::span<const typename E::value_type>()const{return s.coeffs();}
operator span<E,exact_>()const{return s;}
maybe_cached<E,exact_>first(int n)const;
fft::transformed<E>&cache()const{return f;}
};
template<typename S>
concept has_cache_opt=like<S>&&requires(const S&s){
{s.cache_opt()}->std::same_as<std::optional<std::reference_wrapper<fft::transformed<typename S::engine_t>>>>;
};
namespace detail{
template<like S>
std::optional<std::reference_wrapper<fft::transformed<typename S::engine_t>>>cache_of(const S&s){
if constexpr(has_cache<S>)return s.cache();
else if constexpr(has_cache_opt<S>)return s.cache_opt();
else return std::nullopt;
}
}
template<fft::engine E,bool exact_>
struct maybe_cached{
using T=typename E::value_type;
using engine_t=E;
static constexpr bool exact_v=exact_;
span<E,exact_>s;
std::optional<std::reference_wrapper<fft::transformed<E>>>f;
explicit maybe_cached(span<E,exact_>s_):s(s_){}
maybe_cached(span<E,exact_>s_,fft::transformed<E>&f_):s(s_),f(f_){}
maybe_cached(cached_span<E,exact_>c):s(c.s),f(c.f){}
template<like S>requires std::same_as<typename S::engine_t,E>&&(S::exact_v==exact_)
maybe_cached(const S&o):s(o),f(detail::cache_of(o)){}
int len()const{return s.len();}
const T&operator[](int i)const{return s[i];}
operator std::span<const T>()const{return s.coeffs();}
operator span<E,exact_>()const{return s;}
maybe_cached first(int n)const{
return n==len()?*this:maybe_cached(s.first(n));
}
std::optional<std::reference_wrapper<fft::transformed<E>>>cache_opt()const{return f;}
};
template<fft::engine E,bool exact_>
maybe_cached<E,exact_>cached_span<E,exact_>::first(int n)const{
return maybe_cached<E,exact_>(*this).first(n);
}
template<fft::engine E>
struct resized{
using T=typename E::value_type;
using engine_t=E;
static constexpr bool exact_v=false;
trunc<E>s;
std::optional<std::reference_wrapper<fft::transformed<E>>>f;
int len()const{return s.len();}
const T&operator[](int i)const{return s[size_t(i)];}
operator std::span<const T>()const{return std::span<const T>(s);}
operator span<E,false>()const{return s;}
maybe_cached<E,false>first(int n)const{
span<E,false>v=s;
if(n==len()&&f)return{v,f->get()};
return maybe_cached<E,false>(v.first(n));
}
std::optional<std::reference_wrapper<fft::transformed<E>>>cache_opt()const{return f;}
};
template<like S>
resized<typename S::engine_t>with_len(const S&s,int n){
using E=typename S::engine_t;
using T=typename E::value_type;
auto p=s.first(std::min(n,s.len()));
resized<E>r;
r.s.assign(size_t(n),T{});
std::span<const T>pc(p);
std::copy(pc.begin(),pc.end(),r.s.begin());
r.f=detail::cache_of(p);
return r;
}
template<typename S>
concept has_prefix_cache=like<S>&&!S::exact_v&&requires(const S&s,int n){
{s.prefix_cache(n)}->std::same_as<fft::transformed<typename S::engine_t>&>;
};
template<fft::engine E,bool exact_>
struct cached{
using T=typename E::value_type;
using engine_t=E;
static constexpr bool exact_v=exact_;
cached()=default;
cached(vec<E,exact_>&&s_):s(std::move(s_)){}
explicit cached(const vec<E,exact_>&s_):s(s_){}
operator vec<E,exact_>()&&{return std::move(s);}
int len()const{return s.len();}
const vec<E,exact_>&uncached()const{return s;}
const T&operator[](int i)const{return s[size_t(i)];}
auto begin()const{return s.cbegin();}
auto end()const{return s.cend();}
operator span<E,exact_>()const{return s;}
maybe_cached<E,exact_>first(int n)const{
return n==len()?maybe_cached<E,exact_>(s,f):maybe_cached<E,exact_>(s.first(n));
}
fft::transformed<E>&cache()const{return f;}
template<like S>
friend bool operator==(const cached&a,const S&b){
span<E,S::exact_v>bs=b;
return a.len()==bs.len()&&std::equal(a.s.begin(),a.s.end(),bs.begin());
}
private:
vec<E,exact_>s;
mutable fft::transformed<E>f;
};
template<fft::engine E>using cached_exact=cached<E,true>;
template<fft::engine E>using cached_trunc=cached<E,false>;
namespace detail{
template<like S>
cached_span<typename S::engine_t,S::exact_v>as_cached_span(const S&s,fft::transformed<typename S::engine_t>&tmp){
auto co=cache_of(s);
return{s,co?co->get():tmp};
}
}
template<trunc_like S>
trunc<typename S::engine_t>ps_inv(const S&a){
using E=typename S::engine_t;
using T=typename E::value_type;
int N=a.len();
trunc<E>r(size_t(N),T{});
if(N==0)return r;
int s=nextPow2(N);
std::vector<T>b(size_t(s),T{});
b[0]=inv(a[0]);
for(int n=1;n<N;n*=2){
int m=2*n;
auto ta=E::transform(a.first(std::min(N,m)),m);
auto tb=E::transform(std::span<const T>(b).first(n),m);
auto e=fft::buffer_pool<T>::get(m);
E::finish(E::mul(ta,tb,m),e.span());
for(int i=0;i<n;i++)e[i]=T{};
auto te=E::transform(std::span<const T>(e.span()),m);
auto c=fft::buffer_pool<T>::get(m);
E::finish(E::mul(tb,te,m),c.span());
for(int i=n;i<std::min(m,N);i++)b[i]=-c[i];
}
std::copy(b.begin(),b.begin()+N,r.begin());
return r;
}
template<like A>
auto square(const A&a){
using E=typename A::engine_t;
using T=typename E::value_type;
fft::transformed<E>ta_;
auto av=detail::as_cached_span(a,ta_);
if constexpr(A::exact_v){
std::vector<T>coeffs;
fft::transformed<E>f;
fft::square_cached<E>(av,av.cache(),coeffs,f);
cached<E,true>w(exact<E>(std::move(coeffs)));
w.cache()=std::move(f);
return w;
}else{
trunc<E>r(size_t(a.len()),T{});
fft::square<E>(av,av.cache(),std::span<T>(r));
return r;
}
}
template<like A,like B,like C,like D>
requires fft::same_engine<A,B>&&fft::same_engine<A,C>&&fft::same_engine<A,D>
&&A::exact_v&&B::exact_v&&C::exact_v&&D::exact_v
cached<typename A::engine_t,true>multiply_add2(
const A&a,const B&b,const C&c,const D&d){
using E=typename A::engine_t;
using T=typename E::value_type;
fft::transformed<E>ta_,tb_,tc_,td_;
auto av=detail::as_cached_span(a,ta_),bv=detail::as_cached_span(b,tb_);
auto cv=detail::as_cached_span(c,tc_),dv=detail::as_cached_span(d,td_);
std::vector<T>coeffs;
fft::transformed<E>f;
fft::multiply_add2_cached<E>(
av,av.cache(),
bv,bv.cache(),
cv,cv.cache(),
dv,dv.cache(),
coeffs,f
);
cached<E,true>w(exact<E>(std::move(coeffs)));
w.cache()=std::move(f);
return w;
}
template<like A,exact_like B>requires fft::same_engine<A,B>
vec<typename A::engine_t,A::exact_v>middle_product(const A&a,const B&b){
using E=typename A::engine_t;
fft::transformed<E>ta_,tb_;
auto av=detail::as_cached_span(a,ta_);
auto bv=detail::as_cached_span(b,tb_);
return vec<E,A::exact_v>(fft::middle_product<E>(
av,av.cache(),
bv,bv.cache()
));
}
namespace detail{
template<bool ea,bool eb>int product_prec(int la,int lb){
if constexpr(ea&&eb)return la>0&&lb>0?la+lb-1:0;
else return ea?lb:eb?la:std::min(la,lb);
}
template<like S>
auto product_operand(const S&s,int prec,fft::transformed<typename S::engine_t>&tmp){
using E=typename S::engine_t;
if constexpr(has_prefix_cache<S>){
return s.first(std::min(s.len(),nextPow2(prec)));
}else{
span<E,S::exact_v>v=s;
int used=std::min(v.len(),prec);
if(auto co=cache_of(s)){
if(s.len()<=prec||fft::detail::conv_size_for(s.len()+prec-1).n
==fft::detail::conv_size_for(2*prec-1).n){
return cached_span<E,S::exact_v>{v,co->get()};
}
}
return cached_span<E,S::exact_v>{v.first(used),tmp};
}
}
}
template<like A,like B>requires fft::same_engine<A,B>
vec<typename A::engine_t,A::exact_v&&B::exact_v>operator+(const A&a,const B&b){
using T=typename A::engine_t::value_type;
int n=(A::exact_v&&B::exact_v)?std::max(a.len(),b.len())
:A::exact_v?b.len():B::exact_v?a.len():std::min(a.len(),b.len());
vec<typename A::engine_t,A::exact_v&&B::exact_v>r(size_t(n),T(0));
for(int i=0;i<n;i++){
r[i]=(i<a.len()?a[i]:T(0))+(i<b.len()?b[i]:T(0));
}
return r;
}
template<like A,like B>requires fft::same_engine<A,B>
vec<typename A::engine_t,A::exact_v&&B::exact_v>operator-(const A&a,const B&b){
using T=typename A::engine_t::value_type;
int n=(A::exact_v&&B::exact_v)?std::max(a.len(),b.len())
:A::exact_v?b.len():B::exact_v?a.len():std::min(a.len(),b.len());
vec<typename A::engine_t,A::exact_v&&B::exact_v>r(size_t(n),T(0));
for(int i=0;i<n;i++){
r[i]=(i<a.len()?a[i]:T(0))-(i<b.len()?b[i]:T(0));
}
return r;
}
template<like A,like B>requires fft::same_engine<A,B>
auto operator*(const A&a,const B&b){
using E=typename A::engine_t;
using T=typename E::value_type;
constexpr bool ea=A::exact_v,eb=B::exact_v;
int prec=detail::product_prec<ea,eb>(a.len(),b.len());
if(prec==0||a.len()==0||b.len()==0){
if constexpr(ea&&eb)return cached<E,true>{};
else return trunc<E>(size_t(prec),T(0));
}
fft::transformed<E>ta_,tb_;
auto va=detail::product_operand(a,prec,ta_);
auto vb=detail::product_operand(b,prec,tb_);
if constexpr(ea&&eb){
std::vector<T>coeffs;
fft::transformed<E>f;
auto ca=detail::as_cached_span(va,ta_),cb=detail::as_cached_span(vb,tb_);
fft::multiply_cached<E>(
ca,ca.cache(),
cb,cb.cache(),
coeffs,f
);
cached<E,true>w(exact<E>(std::move(coeffs)));
w.cache()=std::move(f);
return w;
}else{
trunc<E>r(size_t(prec),T(0));
auto ca=detail::as_cached_span(va,ta_),cb=detail::as_cached_span(vb,tb_);
fft::multiply<E>(
ca,ca.cache(),
cb,cb.cache(),
std::span<T>(r)
);
return r;
}
}
template<fft::engine E>
struct prefix_cached{
using T=typename E::value_type;
using engine_t=E;
static constexpr bool exact_v=false;
prefix_cached()=default;
prefix_cached(trunc<E>&&s_):s(std::move(s_)){}
explicit prefix_cached(const trunc<E>&s_):s(s_){}
operator trunc<E>()&&{return std::move(s);}
int len()const{return s.len();}
const trunc<E>&uncached()const{return s;}
const T&operator[](int i)const{return s[size_t(i)];}
auto begin()const{return s.cbegin();}
auto end()const{return s.cend();}
operator span<E,false>()const{return s;}
void append(std::span<const T>tail){
s.insert(s.end(),tail.begin(),tail.end());
}
maybe_cached<E,false>first(int k)const{
assert(k<=len());
span<E,false>v=s.first(k);
int n=nextPow2(k);
if(std::min(n,len())==k)return{v,prefix_cache(n)};
return maybe_cached<E,false>(v);
}
fft::transformed<E>&cache()const{return prefix_cache(nextPow2(len()));}
fft::transformed<E>&prefix_cache(int n)const{
assert(n>0&&!(n&(n-1)));
int k=__builtin_ctz(unsigned(n));
if(k>=sz(caches))caches.resize(size_t(k)+1);
auto&c=caches[k];
int e=std::min(n,len());
if(c.len!=e){
c.t=E::transform(s.first(e),2*n);
c.len=e;
}
return c.t;
}
private:
trunc<E>s;
struct entry{fft::transformed<E>t;int len=0;};
mutable std::vector<entry>caches;
};
}
// src/fft/series.hpp
namespace ecnerwala::series{
template<like S>
vec<typename S::engine_t,S::exact_v>stretch(const S&a_,int n){
using E=typename S::engine_t;
span<E,S::exact_v>a=a_;
vec<E,S::exact_v>r(size_t(a.len()));
for(int i=0;i*n<a.len();i++){
r[i*n]=a[i];
}
return r;
}
template<like S>
vec<typename S::engine_t,S::exact_v>deriv_shift(const S&a_){
using E=typename S::engine_t;
span<E,S::exact_v>a=a_;
vec<E,S::exact_v>r(a.begin(),a.end());
for(int i=0;i<r.len();i++){
r[i]*=i;
}
return r;
}
template<like S>
vec<typename S::engine_t,S::exact_v>integ_shift(const S&a_){
using E=typename S::engine_t;
using T=typename E::value_type;
span<E,S::exact_v>a=a_;
assert(a[0]==0);
vec<E,S::exact_v>r(a.begin(),a.end());
T f=1;
for(int i=1;i<r.len();i++){
r[i]*=f;
f*=i;
}
f=inv(f);
for(int i=r.len()-1;i>0;i--){
r[i]*=f;
f*=i;
}
return r;
}
template<like S>
vec<typename S::engine_t,S::exact_v>integ_shift_offset(const S&a_,int offset){
using E=typename S::engine_t;
using T=typename E::value_type;
span<E,S::exact_v>a=a_;
vec<E,S::exact_v>r(a.begin(),a.end());
T f=1;
for(int i=0;i<r.len();i++){
r[i]*=f;
f*=i+offset;
}
assert(f!=0);
f=inv(f);
for(int i=r.len()-1;i>=0;i--){
r[i]*=f;
f*=i+offset;
}
return r;
}
template<trunc_like S>
trunc<typename S::engine_t>deriv_shift_log(const S&a){
return deriv_shift(a)*ps_inv(a);
}
template<trunc_like S>
trunc<typename S::engine_t>ps_log(const S&a){
assert(a[0]==1);
return integ_shift(deriv_shift_log(a));
}
template<trunc_like S>
trunc<typename S::engine_t>ps_exp(const S&a_){
using E=typename S::engine_t;
using T=typename E::value_type;
span<E,false>a=a_;
assert(a.len()>=1);
assert(a[0]==0);
trunc<E>r(1,T(1));r.reserve(size_t(a.len()));
trunc<E>invR(1,T(1));invR.reserve(size_t(a.len()));
while(r.len()<a.len()){
int o_sz=r.len();
int n_sz=std::min(o_sz*2,a.len());
trunc<E>t=deriv_shift(trunc<E>(a.begin(),a.begin()+o_sz));
fft::multiply_circular<E>(std::span<const T>(t),std::span<const T>(r).first(o_sz),std::span<T>(t),o_sz);
t=deriv_shift(r)-t;
t*=invR;
t.resize(size_t(n_sz-o_sz));
trunc<E>v(a.begin()+o_sz,a.begin()+n_sz);
v-=integ_shift_offset(t,o_sz);
v*=r;
r.resize(size_t(n_sz));
std::copy(v.begin(),v.end(),r.begin()+o_sz);
if(r.len()<a.len()){
assert(r.len()==2*invR.len());
int n=invR.len();
int nn=r.len();
trunc<E>tmp(size_t(4)*n);
fft::square<E>(std::span<const T>(invR).first(n),std::span<T>(tmp));
fft::multiply<E>(std::span<const T>(tmp).first(nn),std::span<const T>(r).first(nn),std::span<T>(tmp));
invR.resize(size_t(nn));
for(int i=n;i<nn;i++)invR[i]=-tmp[i];
}
}
return r;
}
template<trunc_like S>
trunc<typename S::engine_t>ps_pow_monic(const S&a_,typename S::engine_t::value_type k){
using E=typename S::engine_t;
span<E,false>a=a_;
if(a.len()==0)return{};
assert(a[0]==1);
trunc<E>l=ps_log(a_);
l*=k;
return ps_exp(l);
}
template<trunc_like S>
trunc<typename S::engine_t>ps_pow(const S&a_,int64_t k){
using E=typename S::engine_t;
using T=typename E::value_type;
span<E,false>a=a_;
assert(k>=0);
if(k==0){
trunc<E>r(size_t(a.len()),T(0));
if(r.len()>0)r[0]=T(1);
return r;
}
int st=0;
while(st<a.len()&&a[st]==0)st++;
if(st>0&&k>(a.len()-1)/st){
return trunc<E>(size_t(a.len()),T(0));
}
trunc<E>r(a.begin()+st,a.end()-(st*(k-1)));
T leading_coeff=r[0];
r*=inv(leading_coeff);
r=ps_pow_monic(r,T(k));
r*=power(leading_coeff,k);
r.insert(r.begin(),size_t(st*k),T(0));
assert(r.len()==a.len());
return r;
}
template<trunc_like S>
trunc<typename S::engine_t>to_newton_sums(const S&a,int deg){
auto r=deriv_shift_log(a);
r[0]=deg;
for(int i=1;i<r.len();i++)r[i]=-r[i];
return r;
}
template<trunc_like S>
trunc<typename S::engine_t>from_newton_sums(const S&s_,int deg){
using E=typename S::engine_t;
span<E,false>s=s_;
assert(s[0]==deg);
trunc<E>r(s.begin(),s.end());
r[0]=0;
for(int i=1;i<r.len();i++)r[i]=-r[i];
return ps_exp(integ_shift(std::move(r)));
}
template<trunc_like S>
trunc<typename S::engine_t>euler_transform(const S&a){
using E=typename S::engine_t;
trunc<E>r=deriv_shift(a);
std::vector<bool>is_prime(size_t(r.len()),true);
for(int p=2;p<r.len();p++){
if(!is_prime[p])continue;
for(int i=1;i*p<r.len();i++){
r[i*p]+=r[i];
is_prime[i*p]=false;
}
}
return ps_exp(integ_shift(r));
}
template<trunc_like S>
trunc<typename S::engine_t>inverse_euler_transform(const S&a){
using E=typename S::engine_t;
trunc<E>r=deriv_shift(ps_log(a));
std::vector<bool>is_prime(size_t(r.len()),true);
for(int p=2;p<r.len();p++){
if(!is_prime[p])continue;
for(int i=(r.len()-1)/p;i>=1;i--){
r[i*p]-=r[i];
is_prime[i*p]=false;
}
}
return integ_shift(r);
}
template<fft::engine E>struct packed_bivariate{
using T=typename E::value_type;
int L,l;
std::vector<T>c;
packed_bivariate(int L_,std::span<const T>g):L(L_),l(0),c(size_t(4)<<L){
c[0]=T(1);
for(int i=0;i<sz(g);i++)c[(2<<L)+i]=-g[i];
}
fft::transformed<E>advance(){
int B=4<<L;
auto tq=E::transform(std::span<const T>(c),B);
auto tn=E::negate_arg(tq,B);
E::finish(
E::downsample(E::mul(tq,tn,B),B/2,false),
std::span<T>(c).first(B/2)
);
l++;
for(int i=0;i<(2<<(L-l));i++){
c[(2<<L)+i]=c[i];
c[i]=T(0);
}
c[2<<L]-=T(1);
c[0]=T(1);
std::fill(c.begin()+(2<<L)+(1<<(L-l)),c.end(),T(0));
for(int i=0;i<(2<<L);i+=2<<(L-l)){
for(int j=0;j<(1<<(L-l));j++){
c[i+(1<<(L-l))+j]=T(0);
}
}
return tn;
}
};
template<trunc_like SF,trunc_like SG>requires fft::same_engine<SF,SG>
trunc<typename SF::engine_t>ps_compose(const SF&f_,const SG&g_){
using E=typename SF::engine_t;
using T=typename E::value_type;
span<E,false>f=f_;
span<E,false>g=g_;
if(g.len()==0)return{};
int m=f.len();
int n=g.len();
int L=__builtin_ctz(unsigned(nextPow2(n)));
int B=4<<L;
packed_bivariate<E>Q(L,g.coeffs());
std::vector<fft::transformed<E>>tneg;
tneg.reserve(L);
for(int l=1;l<=L;l++)tneg.push_back(Q.advance());
trunc<E>P;
{
P=trunc<E>(f.begin(),f.end());
std::reverse(P.begin(),P.end());
trunc<E>QL((1<<L)+1);
for(int i=0;i<=(1<<L);i++){
QL[i]=Q.c[2*i];
}
QL.resize(size_t(m),T(0));
P*=ps_inv(QL);
std::reverse(P.begin(),P.end());
P.resize(size_t(1)<<L,T(0));
std::reverse(P.begin(),P.end());
P.resize(size_t(B),T(0));
for(int i=(1<<L)-1;i>0;i--){
P[2*i]=P[i];
P[i]=T(0);
}
}
for(int l=L-1;l>=0;l--){
for(int i=(2<<L)-1;i>0;i--){
T v=P[i];
P[2*i]=((2*i)&(1<<(L-l)))?T(0):v;
P[i]=T(0);
}
auto tp=E::transform(std::span<const T>(P),B);
E::finish(E::mul(tneg[l],tp,B),std::span<T>(P));
for(int i=0;i<(2<<L);i++){
P[i]=P[(2<<L)+i];
P[(2<<L)+i]=T(0);
}
}
return trunc<E>(P.begin(),P.begin()+n);
}
template<exact_like P,exact_like Q>requires fft::same_engine<P,Q>
P::engine_t::value_type kth_term_of_rational_function(
const P&p,
const Q&q,
uint64_t k
){
using E=P::engine_t;
using T=E::value_type;
assert(q.len()>0&&q[0]!=T(0));
if(p.len()==0)return T(0);
int d=std::max(p.len()+1,q.len());
assert(d>=2);
int n=nextPow2((d-1)+d-1);
fft::transformed<E>tq,tp;
if(auto cq=detail::cache_of(q)){E::extend_to(cq->get(),n,q);tq=cq->get();}
if(auto cp=detail::cache_of(p)){E::extend_to(cp->get(),n,p);tp=cp->get();}
std::vector<T>p_buf(d-1,T(0));
std::ranges::copy(std::span<const T>(p),p_buf.begin());
std::vector<T>q_buf(d,T(0));
std::ranges::copy(std::span<const T>(q),q_buf.begin());
while(k>0){
E::extend_to(tq,n,q_buf);
auto tnq=E::negate_arg(tq,n);
E::extend_to(tp,n,p_buf);
auto ntp=E::downsample(E::mul(tp,tnq,n),n/2,bool(k&1));
assert(ntp.size()==n/2);
if constexpr(std::same_as<typename E::product,typename E::transformed>){
tp=ntp;
}else{
tp={};
}
E::finish(std::move(ntp),std::span(p_buf));
k>>=1;
if(!k){
q_buf[0]*=q_buf[0];
break;
}
auto ntq=E::downsample(E::mul(tq,tnq,n),n/2,false);
assert(ntq.size()==n/2);
if constexpr(std::same_as<typename E::product,typename E::transformed>){
tq=ntq;
}else{
tq={};
}
if(n/2==d-1){
T v0=q_buf[0]*q_buf[0];
E::finish(std::move(ntq),std::span(q_buf).first(d-1));
q_buf[d-1]=std::exchange(q_buf[0],v0)-v0;
}else{
E::finish(std::move(ntq),std::span(q_buf));
}
}
return p_buf[0]*inv(q_buf[0]);
}
template<trunc_like S,exact_like Q>requires fft::same_engine<S,Q>
S::engine_t::value_type kth_term_of_linear_recurrence(
const S&s,
const Q&q,
uint64_t k
){
using E=S::engine_t;
using T=E::value_type;
assert(q.len()>0&&q[0]!=T(0));
assert(s.len()>=q.len()-1);
fft::transformed<E>tq;
auto q_cached=detail::as_cached_span(q,tq);
span<E,false>sv=s;
auto p=exact<E>(sv.first(q.len()-1)*q_cached);
return kth_term_of_rational_function(p,q_cached,k);
}
}
#pragma GCC diagnostic pop
// clang-format on
// @formatter:on

#pragma once

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

#include "fft/series_core.hpp"

// ==== analytic ops ====
// Free functions over series-like operands; each borrows the operand's span
// and writes a fresh result.
// TODO: reuse/populate the operands' whole/prefix transform caches

namespace ecnerwala::series {

template <like S>
vec<typename S::engine_t, S::exact_v> stretch(const S& a_, int n) {
	using E = typename S::engine_t;
	span<E, S::exact_v> a = a_;
	vec<E, S::exact_v> r(size_t(a.len()));
	for (int i = 0; i*n < a.len(); i++) {
		r[i*n] = a[i];
	}
	return r;
}
template <like S>
vec<typename S::engine_t, S::exact_v> deriv_shift(const S& a_) {
	using E = typename S::engine_t;
	span<E, S::exact_v> a = a_;
	vec<E, S::exact_v> r(a.begin(), a.end());
	for (int i = 0; i < r.len(); i++) {
		r[i] *= i;
	}
	return r;
}
template <like S>
vec<typename S::engine_t, S::exact_v> integ_shift(const S& a_) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, S::exact_v> a = a_;
	assert(a[0] == 0);
	vec<E, S::exact_v> r(a.begin(), a.end());
	T f = 1;
	for (int i = 1; i < r.len(); i++) {
		r[i] *= f;
		f *= i;
	}
	f = inv(f);
	for (int i = r.len() - 1; i > 0; i--) {
		r[i] *= f;
		f *= i;
	}
	return r;
}
template <like S>
vec<typename S::engine_t, S::exact_v> integ_shift_offset(const S& a_, int offset) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, S::exact_v> a = a_;
	vec<E, S::exact_v> r(a.begin(), a.end());
	T f = 1;
	for (int i = 0; i < r.len(); i++) {
		r[i] *= f;
		f *= i + offset;
	}
	assert(f != 0);
	f = inv(f);
	for (int i = r.len() - 1; i >= 0; i--) {
		r[i] *= f;
		f *= i + offset;
	}
	return r;
}
template <trunc_like S>
trunc<typename S::engine_t> deriv_shift_log(const S& a) {
	return deriv_shift(a) * ps_inv(a);
}
template <trunc_like S>
trunc<typename S::engine_t> ps_log(const S& a) {
	assert(a[0] == 1);
	return integ_shift(deriv_shift_log(a));
}
template <trunc_like S>
trunc<typename S::engine_t> ps_exp(const S& a_) {
	// See https://mathexp.eu/bostan/publications/BoSc09a.pdf for details
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, false> a = a_;
	assert(a.len() >= 1);
	assert(a[0] == 0);
	trunc<E> r(1, T(1)); r.reserve(size_t(a.len()));
	trunc<E> invR(1, T(1)); invR.reserve(size_t(a.len()));
	while (r.len() < a.len()) {
		int o_sz = r.len();
		int n_sz = std::min(o_sz * 2, a.len());
		trunc<E> t = deriv_shift(trunc<E>(a.begin(), a.begin() + o_sz));
		fft::multiply_circular<E>(std::span<const T>(t), std::span<const T>(r).first(o_sz), std::span<T>(t), o_sz);
		t = deriv_shift(r) - t;
		t *= invR;
		t.resize(size_t(n_sz - o_sz));
		trunc<E> v(a.begin() + o_sz, a.begin() + n_sz);
		v -= integ_shift_offset(t, o_sz);
		v *= r;
		r.resize(size_t(n_sz));
		std::copy(v.begin(), v.end(), r.begin() + o_sz);
		if (r.len() < a.len()) {
			// double invR via a Newton step
			assert(r.len() == 2 * invR.len());
			int n = invR.len();
			int nn = r.len();
			trunc<E> tmp(size_t(4) * n);
			fft::square<E>(std::span<const T>(invR).first(n), std::span<T>(tmp));
			fft::multiply<E>(std::span<const T>(tmp).first(nn), std::span<const T>(r).first(nn), std::span<T>(tmp));
			invR.resize(size_t(nn));
			for (int i = n; i < nn; i++) invR[i] = -tmp[i];
		}
	}
	return r;
}
template <trunc_like S>
trunc<typename S::engine_t> ps_pow_monic(const S& a_, typename S::engine_t::value_type k) {
	using E = typename S::engine_t;
	span<E, false> a = a_;
	if (a.len() == 0) return {};
	assert(a[0] == 1);
	trunc<E> l = ps_log(a_);
	l *= k;
	return ps_exp(l);
}
template <trunc_like S>
trunc<typename S::engine_t> ps_pow(const S& a_, int64_t k) {
	using E = typename S::engine_t;
	using T = typename E::value_type;
	span<E, false> a = a_;
	assert(k >= 0);
	if (k == 0) {
		trunc<E> r(size_t(a.len()), T(0));
		if (r.len() > 0) r[0] = T(1);
		return r;
	}

	int st = 0;
	while (st < a.len() && a[st] == 0) st++;

	if (st > 0 && k > (a.len() - 1) / st) {
		return trunc<E>(size_t(a.len()), T(0));
	}

	trunc<E> r(a.begin() + st, a.end() - (st * (k-1)));
	T leading_coeff = r[0];
	r *= inv(leading_coeff);
	r = ps_pow_monic(r, T(k));
	r *= power(leading_coeff, k);
	r.insert(r.begin(), size_t(st * k), T(0));
	assert(r.len() == a.len());
	return r;
}

template <trunc_like S>
trunc<typename S::engine_t> to_newton_sums(const S& a, int deg) {
	auto r = deriv_shift_log(a);
	r[0] = deg;
	for (int i = 1; i < r.len(); i++) r[i] = -r[i];
	return r;
}
template <trunc_like S>
trunc<typename S::engine_t> from_newton_sums(const S& s_, int deg) {
	using E = typename S::engine_t;
	span<E, false> s = s_;
	assert(s[0] == deg);
	trunc<E> r(s.begin(), s.end());
	r[0] = 0;
	for (int i = 1; i < r.len(); i++) r[i] = -r[i];
	return ps_exp(integ_shift(std::move(r)));
}

// Calculates prod 1/(1-x^i)^{a[i]}
template <trunc_like S>
trunc<typename S::engine_t> euler_transform(const S& a) {
	using E = typename S::engine_t;
	trunc<E> r = deriv_shift(a);
	std::vector<bool> is_prime(size_t(r.len()), true);
	for (int p = 2; p < r.len(); p++) {
		if (!is_prime[p]) continue;
		for (int i = 1; i*p < r.len(); i++) {
			r[i*p] += r[i];
			is_prime[i*p] = false;
		}
	}
	return ps_exp(integ_shift(r));
}
template <trunc_like S>
trunc<typename S::engine_t> inverse_euler_transform(const S& a) {
	using E = typename S::engine_t;
	trunc<E> r = deriv_shift(ps_log(a));
	std::vector<bool> is_prime(size_t(r.len()), true);
	for (int p = 2; p < r.len(); p++) {
		if (!is_prime[p]) continue;
		for (int i = (r.len()-1)/p; i >= 1; i--) {
			r[i*p] -= r[i];
			is_prime[i*p] = false;
		}
	}
	return integ_shift(r);
}

// Helper packed bivariate buffer for Kinoshita-Li composition (arXiv:2404.05177).
//
// The motivation is performing Bostan-Mori (Graeffe root-squaring) to compute
// something like [x^n] P / Q_0(x, y) with deg_y(Q_0) = 1 and deg_x(Q_0) = n.
//
// In each step, we want to compute Q_{i+1}(x^2, y) = Q_i(x, y) * Q_i(-x, y).
// This doubles the degree of y and also lets us truncate x at half the previous
// degree, leaving the total size invariant.
//
// We will store Q as a packed buffer with x as the inner dimension to facilitate easy Q(-x) substitution.
// The inner span will be 2*deg(x), and the outer span will be 2*deg(y).
// As we advance, we will also return the cached transform of Q_i(-x, y) for the caller to use in the numerator.
template <fft::engine E> struct packed_bivariate {
	using T = typename E::value_type;
	int L, l;
	std::vector<T> c;

	// Q_0 = 1 - y g(x), deg g < n <= 2^L
	packed_bivariate(int L_, std::span<const T> g) : L(L_), l(0), c(size_t(4) << L) {
		c[0] = T(1);
		for (int i = 0; i < sz(g); i++) c[(2 << L) + i] = -g[i];
	}

	fft::transformed<E> advance() {
		int B = 4 << L;
		auto tq = E::transform(std::span<const T>(c), B);
		auto tn = E::negate_arg(tq, B);
		E::finish(
			E::downsample(E::mul(tq, tn, B), B/2, false),
			std::span<T>(c).first(B/2)
		);
		l++;
		// undo the circular wraparound using monicity in y
		for (int i = 0; i < (2 << (L - l)); i++) {
			c[(2 << L) + i] = c[i];
			c[i] = T(0);
		}
		c[2 << L] -= T(1);
		c[0] = T(1);
		// zero x coefficients beyond the level's truncation mod x^(2^(L-l))
		std::fill(c.begin() + (2 << L) + (1 << (L - l)), c.end(), T(0));
		for (int i = 0; i < (2 << L); i += 2 << (L - l)) {
			for (int j = 0; j < (1 << (L - l)); j++) {
				c[i + (1 << (L - l)) + j] = T(0);
			}
		}
		return tn;
	}
};

// Calculates f(g(x)) mod x^n where deg(g) == n
template <trunc_like SF, trunc_like SG> requires fft::same_engine<SF, SG>
trunc<typename SF::engine_t> ps_compose(const SF& f_, const SG& g_) {
	using E = typename SF::engine_t;
	using T = typename E::value_type;
	span<E, false> f = f_;
	span<E, false> g = g_;
	if (g.len() == 0) return {};

	int m = f.len();
	int n = g.len();

	// https://arxiv.org/pdf/2404.05177
	// Consider P(y) = f(1/y) has terms from y^{-(m-1)}...y^0 (Laurent series)
	// We want [y^0] P(y) / (1 - y g(x))
	// Let Q_0 = 1 - yg(x)
	// Q_{i+1}(x^2, y) = Q_i(x, y) * Q_i(-x, y) mod x^{ceil(n / 2^i)}
	// deg_y(Q_i) = 2^i, deg_x(Q_i) = ceil(n / 2^i) - 1
	//
	// [y^0] P(y) / Q_l(x^2^l, y) * Q_{l-1}(-x^2^{l-1}, y) * Q_{l-2}(-x^2^{l-2}, y) * ... * Q_0(-x, y)
	// The total y deg of Q_{k-1} ... Q_0 is 2^k-1
	int L = __builtin_ctz(unsigned(nextPow2(n)));
	int B = 4 << L;
	packed_bivariate<E> Q(L, g.coeffs());
	// tneg[l] is the transform of Q_l(-x, y), reused by the pushdown pass below
	std::vector<fft::transformed<E>> tneg;
	tneg.reserve(L);
	for (int l = 1; l <= L; l++) tneg.push_back(Q.advance());
	trunc<E> P;
	{
		P = trunc<E>(f.begin(), f.end());
		std::reverse(P.begin(), P.end());
		trunc<E> QL((1 << L) + 1);
		for (int i = 0; i <= (1 << L); i++) {
			QL[i] = Q.c[2 * i];
		}
		QL.resize(size_t(m), T(0));
		P *= ps_inv(QL);
		std::reverse(P.begin(), P.end());
		P.resize(size_t(1) << L, T(0));
		std::reverse(P.begin(), P.end());
		P.resize(size_t(B), T(0));
		for (int i = (1 << L) - 1; i > 0; i--) {
			P[2*i] = P[i];
			P[i] = T(0);
		}
	}
	for (int l = L-1; l >= 0; l--) {
		// Spread it out, clear the high terms
		for (int i = (2 << L) - 1; i > 0; i--) {
			T v = P[i];
			P[2*i] = ((2*i) & (1 << (L-l))) ? T(0) : v;
			P[i] = T(0);
		}
		auto tp = E::transform(std::span<const T>(P), B);
		E::finish(E::mul(tneg[l], tp, B), std::span<T>(P));
		for (int i = 0; i < (2 << L); i++) {
			P[i] = P[(2 << L) + i];
			P[(2 << L) + i] = T(0);
		}
	}
	return trunc<E>(P.begin(), P.begin() + n);
}

// [x^k] p(x)/q(x) (Bostan-Mori) for an exact rational function.
template <exact_like P, exact_like Q> requires fft::same_engine<P, Q>
P::engine_t::value_type kth_term_of_rational_function(
	const P& p,
	const Q& q,
	uint64_t k
) {
	using E = P::engine_t;
	using T = E::value_type;

	assert(q.len() > 0 && q[0] != T(0));
	// Check this here so we avoid accessing p[0]
	if (p.len() == 0) return T(0);

	// Size up in a pretty conservative way
	int d = std::max(p.len() + 1, q.len());
	assert(d >= 2);

	int n = nextPow2((d-1) + d - 1); // >= d

	// Seed the loop transforms from any whole caches; the buffers below hold the
	// current p, q (zero-padded, which extend_to tolerates).
	fft::transformed<E> tq, tp;
	if (auto cq = detail::cache_of(q)) { E::extend_to(cq->get(), n, q); tq = cq->get(); }
	if (auto cp = detail::cache_of(p)) { E::extend_to(cp->get(), n, p); tp = cp->get(); }

	std::vector<T> p_buf(d-1, T(0));
	std::ranges::copy(std::span<const T>(p), p_buf.begin());
	std::vector<T> q_buf(d, T(0));
	std::ranges::copy(std::span<const T>(q), q_buf.begin());

	while (k > 0) {
		E::extend_to(tq, n, q_buf);
		auto tnq = E::negate_arg(tq, n);
		E::extend_to(tp, n, p_buf);

		// P <- downsample(P(x) * Q(-x))
		auto ntp = E::downsample(E::mul(tp, tnq, n), n/2, bool(k & 1));
		assert(ntp.size() == n/2);
		if constexpr (std::same_as<typename E::product, typename E::transformed>) {
			tp = ntp;
		} else {
			tp = {};
		}
		E::finish(std::move(ntp), std::span(p_buf));
		k >>= 1;

		// Save the last iteration if we're done
		if (!k) {
			// HACK: fix the constant coefficient of q only
			q_buf[0] *= q_buf[0];
			break;
		}

		// Q <- downsample(Q(x) * Q(-x))
		auto ntq = E::downsample(E::mul(tq, tnq, n), n/2, false);
		assert(ntq.size() == n/2);
		if constexpr (std::same_as<typename E::product, typename E::transformed>) {
			tq = ntq;
		} else {
			tq = {};
		}
		if (n/2 == d-1) {
			// Fix the wraparound
			T v0 = q_buf[0] * q_buf[0];
			E::finish(std::move(ntq), std::span(q_buf).first(d-1));
			q_buf[d-1] = std::exchange(q_buf[0], v0) - v0;
		} else {
			E::finish(std::move(ntq), std::span(q_buf));
		}
	}
	return p_buf[0] * inv(q_buf[0]);
}

// Find the kth term of linearly recurrent sequence S with char poly Q and len(S) >= len(Q)-1
template <trunc_like S, exact_like Q> requires fft::same_engine<S, Q>
S::engine_t::value_type kth_term_of_linear_recurrence(
	const S& s,
	const Q& q,
	uint64_t k
) {
	using E = S::engine_t;
	using T = E::value_type;

	assert(q.len() > 0 && q[0] != T(0));
	assert(s.len() >= q.len()-1);

	// Don't even bother with P so we don't have to do truncation checks
	// TODO: Could use generic multiply for this whole part?
	fft::transformed<E> tq;
	auto q_cached = detail::as_cached_span(q, tq);

	// Compute the prefix and then hard-cast it to exact
	span<E, false> sv = s;
	auto p = exact<E>(sv.first(q.len()-1) * q_cached);
	return kth_term_of_rational_function(p, q_cached, k);
}

/* namespace ecnerwala::series */ }

Back to top page