ecnerwala's competitive programming library
// competitive-verifier: PROBLEM https://judge.yosupo.jp/problem/dirichlet_inverse_and_prefix_sums
#include <bits/stdc++.h>
#include <cassert>
#include "dirichlet_series.hpp"
#include "modnum.hpp"
int main() {
std::ios_base::sync_with_stdio(false), std::cin.tie(nullptr);
int T; std::cin >> T;
while (T--) {
int64_t N; std::cin >> N;
static dirichlet_series::div_vector_layout layout;
layout = N;
using num = modnum<998244353>;
using ds_prefix = dirichlet_series::prefix<layout, num>;
ds_prefix F;
for (int i = 1; i < layout.len; i++) std::cin >> F.st[i];
ds_prefix delta([&](int64_t) { return 1; });
ds_prefix H = delta / F;
for (int i = 1; i < layout.len; i++) std::cout << H.st[i] << " \n"[i+1==layout.len];
}
return 0;
}
#include <bits/stdc++.h>
#line 1 "verify/dirichlet_inverse_and_prefix_sums.test.cpp"
// competitive-verifier: PROBLEM https://judge.yosupo.jp/problem/dirichlet_inverse_and_prefix_sums
#line 5 "verify/dirichlet_inverse_and_prefix_sums.test.cpp"
#line 2 "src/dirichlet_series.hpp"
#line 8 "src/dirichlet_series.hpp"
namespace dirichlet_series {
inline int inv(int v) {
assert(v == 1);
return 1;
}
inline int64_t inv(int64_t v) {
assert(v == 1);
return 1;
}
constexpr int64_t floor_sqrt(int64_t N) {
assert(N >= 0);
if (N == 0) return 0;
int64_t a = N;
while (true) {
int64_t b = N/a;
assert(a >= b);
if (a-b <= 1) return b;
a = (a+b+1)>>1;
}
}
class div_vector_layout {
public:
int64_t N;
int64_t rt = floor_sqrt(N);
int len = int(2 * rt + (rt * (rt+1) <= N));
constexpr div_vector_layout(int64_t N_ = 1) : N(N_) {}
constexpr int get_value_bucket(int64_t a) const {
return a <= rt ? int(a) : len - int(N/a);
}
constexpr int64_t get_bucket_bound(int i) const {
return i <= rt ? i : N/(len-i);
}
};
template <const div_vector_layout& layout, typename T> class div_vector {
public:
// Let's just make everything public, getters and setters are too much work
T* st = new T[layout.len+1]{}; // Allocate one extra on each side
T* en = st + layout.len;
div_vector() = default;
/* Rule of 5 declarations */
div_vector(div_vector const& o) {
std::copy(o.st, o.en, st);
}
div_vector& operator = (div_vector const& o) {
std::copy(o.st, o.en, st);
return *this;
}
friend void swap(div_vector& a, div_vector& b) {
std::swap(a.st, b.st);
std::swap(a.en, b.en);
}
div_vector(div_vector && o) : st(nullptr), en(nullptr) {
swap(*this, o);
}
div_vector& operator = (div_vector && o) {
swap(*this, o);
return *this;
}
~div_vector() { delete[] st; }
T& operator [] (int64_t v) { return st[layout.get_value_bucket(v)]; }
T& operator [] (int64_t v) const { return st[layout.get_value_bucket(v)]; }
};
template <div_vector_layout const& layout, typename T, typename Derived> class vectorspace_mixin {
private:
Derived& underlying() {
return static_cast<Derived&>(*this);
}
Derived const& underlying() const {
return static_cast<Derived const&>(*this);
}
public:
friend Derived operator + (Derived&& a) {
for (int64_t i = 1; i < layout.len; i++) {
a.st[i] = +a.st[i];
}
return a;
}
friend Derived operator + (Derived const& a) { return +Derived(a); }
friend Derived operator - (Derived && a) {
for (int64_t i = 1; i < layout.len; i++) {
a.st[i] = -a.st[i];
}
return a;
}
friend Derived operator - (Derived const& a) { return -Derived(a); }
Derived& operator += (Derived const& o) {
for (int64_t i = 1; i < layout.len; i++) {
underlying().st[i] += o.st[i];
}
return underlying();
}
friend Derived operator + (Derived && a, Derived const& b) { return a += b; }
friend Derived operator + (Derived const& a, Derived && b) {
for (int64_t i = 1; i < layout.len; i++) {
b.st[i] = a.st[i] + b.st[i];
}
return b;
}
friend Derived operator + (Derived && a, Derived && b) { return std::move(a) + b; }
friend Derived operator + (Derived const& a, Derived const& b) { return Derived(a) + b; }
template <typename F> Derived& operator += (F f) {
for (int64_t i = 1; i < layout.len; i++) {
underlying().st[i] += f(layout.get_bucket_bound(i));
}
return underlying();
}
Derived& operator -= (Derived const& o) {
for (int64_t i = 1; i < layout.len; i++) {
underlying().st[i] -= o.st[i];
}
return underlying();
}
friend Derived operator - (Derived && a, Derived const& b) { return a -= b; }
friend Derived operator - (Derived const& a, Derived && b) {
for (int64_t i = 1; i < layout.len; i++) {
b.st[i] = a.st[i] - b.st[i];
}
return b;
}
friend Derived operator - (Derived && a, Derived && b) { return std::move(a) - b; }
friend Derived operator - (Derived const& a, Derived const& b) { return Derived(a) - b; }
template <typename F> Derived& operator -= (F f) {
for (int64_t i = 1; i < layout.len; i++) {
underlying().st[i] -= f(layout.get_bucket_bound(i));
}
return underlying();
}
Derived& operator *= (T const& t) {
for (int64_t i = 1; i < layout.len; i++) {
underlying().st[i] *= t;
}
return underlying();
}
friend Derived operator * (Derived && a, T const& t) { return a *= t; }
friend Derived operator * (Derived const& a, T const& t) { return Derived(a) * t; }
// Just in case, don't assume multiplication is commutative.
friend Derived operator * (T const& t, Derived && a) {
for (int64_t i = 1; i < layout.len; i++) {
a.st[i] = t * a.st[i];
}
return a;
}
friend Derived operator * (T const& t, Derived const& a) { return t * Derived(a); }
Derived& operator /= (T const& t) {
for (int64_t i = 1; i < layout.len; i++) {
underlying().st[i] /= t;
}
return underlying();
}
friend Derived operator / (Derived && a, T const& t) { return a /= t; }
friend Derived operator / (Derived const& a, T const& t) { return Derived(a) / t; }
};
template <div_vector_layout const& layout, typename T> class values;
template <div_vector_layout const& layout, typename T> class prefix;
template <div_vector_layout const& layout, typename T> class bit;
template <div_vector_layout const& layout, typename T> class values : public div_vector<layout, T>, public vectorspace_mixin<layout, T, values<layout, T>> {
public:
values() = default;
template <typename F, std::enable_if_t<std::is_invocable_r_v<T, F, int64_t, int64_t>, bool> = true>
values(F f) {
for (int i = 1; i < layout.len; i++) {
this->st[i] = f(layout.get_bucket_bound(i-1), layout.get_bucket_bound(i));
}
}
template <typename U> explicit values(values<layout, U> const& o) {
for (int i = 1; i < layout.len; i++) {
this->st[i] = T(o.st[i]);
}
}
explicit values(prefix<layout, T> && o) : div_vector<layout, T>(static_cast<div_vector<layout, T>&&>(std::move(o))) {
for (int i = layout.len - 1; i > 1; i--) {
this->st[i] -= this->st[i-1];
}
}
explicit values(prefix<layout, T> const& o) {
for (int i = layout.len - 1; i > 1; i--) {
this->st[i] = o.st[i] - o.st[i-1];
}
this->st[1] = o.st[1];
}
};
template <div_vector_layout const& layout, typename T> class prefix : public div_vector<layout, T>, public vectorspace_mixin<layout, T, prefix<layout, T>> {
public:
prefix() = default;
template <typename F, std::enable_if_t<std::is_invocable_r_v<T, F, int64_t>, bool> = true>
prefix(F f) {
for (int i = 1; i < layout.len; i++) {
this->st[i] = f(layout.get_bucket_bound(i));
}
}
template <typename U> explicit prefix(prefix<layout, U> const& o) {
for (int i = 1; i < layout.len; i++) {
this->st[i] = T(o.st[i]);
}
}
explicit prefix(values<layout, T> && o) : div_vector<layout, T>(static_cast<div_vector<layout, T>&&>(std::move(o))) {
for (int i = 2; i < layout.len; i++) {
this->st[i] += this->st[i-1];
}
}
explicit prefix(values<layout, T> const& o) {
T pref = this->st[1] = o.st[1];
for (int i = 2; i < layout.len; i++) {
this->st[i] = (pref += o.st[i]);
}
}
private:
// This essentially runs *this += a * b, except it doesn't convolve any
// terms involving 1*i and leaves those for the user-provided function f.
// (f is called for each i in [2, layout.len-1].) This allows us to
// easily implement multiplication or division or sqrt. (Note that a or b
// are allowed to be equal to this.)
template <typename F>
void convolve_helper(prefix const& a, prefix const& b, F f) {
// We roughly want to apply this[N/z] += a_val[x] * b_val[y] for all xyz <= N
//
// We'll split into the following cases (WLOG x <= y):
// 0a. x = 1 or y = 1
// 0b. x = y > 1
// 1. x < y <= z <= N/x/y
// 2. max(x, z) < y <= N/x/z
T cur_sum = a.st[1] * b.st[1];
for (int i = 2; i < layout.len; i++) {
cur_sum += this->st[i];
// Case 2: max(x, z) < y <= N/x/z
// x^2 <= N / z
// x <= N / z / z
if (i >= layout.len - layout.rt) {
int z = int(layout.len - i);
assert(z <= layout.rt);
int64_t rt_over_z = layout.rt/z;
int64_t N_over_z = layout.N/z;
int64_t x_max = N_over_z/(z+1);
T tot_val = T();
for (int64_t x = 2; x * (x+1) <= N_over_z && x <= x_max; x++) {
// ylo = std::max(x, z)
int ylo_idx = std::max(int(x), z);
// yhi = N / x / z
int yhi_idx = int(x <= rt_over_z ? layout.len - x * z : N_over_z / x);
assert(ylo_idx < yhi_idx);
T ax = a.st[x] - a.st[x-1];
T bx = b.st[x] - b.st[x-1];
T ay = a.st[yhi_idx] - a.st[ylo_idx];
T by = b.st[yhi_idx] - b.st[ylo_idx];
T v = ax * by + ay * bx;
tot_val += v;
}
cur_sum += tot_val;
if (i+1 < layout.len) {
this->st[i+1] -= tot_val;
}
}
this->st[i] = f(i, cur_sum);
T ai = a.st[i] - a.st[i-1];
T bi = b.st[i] - b.st[i-1];
// Case 0a: x = 1
cur_sum += ai * b.st[1] + a.st[1] * bi;
if (i <= layout.rt) {
// Case 1: x < y <= z <= N/x/y (y = i)
// xy <= z <= N/y
int64_t rt_over_i = layout.rt / i;
int64_t N_over_i = layout.N / i;
int x_max = int(std::min<int64_t>(N_over_i / i, i-1));
T tot_sub = T();
for (int x = 2; x <= x_max; x++) {
T v;
v = ai * (b.st[x] - b.st[x-1]) + (a.st[x] - a.st[x-1]) * bi;
int zlo_idx = int(x <= rt_over_i ? x * i : layout.len - (N_over_i / x));
this->st[zlo_idx] += v;
tot_sub += v;
}
this->en[-(i-1)] -= tot_sub;
// Case 0b: x = y > 1
{
int zlo_idx = int(i <= rt_over_i ? i * i : layout.len - (N_over_i / i));
this->st[zlo_idx] += ai * bi;
}
}
}
}
public:
friend prefix operator * (prefix const& a, prefix const& b) {
prefix r;
r.st[1] = a.st[1] * b.st[1];
r.convolve_helper(a, b, [&](int i, T cur_sum) -> T {
return cur_sum + (a.st[i] - a.st[i-1]) * b.st[1] + a.st[1] * (b.st[i] - b.st[i-1]);
});
return r;
}
prefix& operator *= (const prefix& o) { return *this = *this * o; }
friend T get_conv_N(prefix const& a, prefix const& b) {
T ans = a.st[1] * b.en[-1];
for (int i = 2; i <= layout.len; i++) {
ans += (a.st[i] - a.st[i-1]) * b.en[-i];
}
return ans;
}
friend prefix operator / (prefix const& a, prefix const& b) {
prefix r;
T inv_b1 = inv(b.st[1]);
r.st[1] = a.st[1] * inv_b1;
r.convolve_helper(r, b, [&](int i, T cur_sum) -> T {
return (a.st[i] - (cur_sum + r.st[1] * (b.st[i] - b.st[i-1]))) * inv_b1 + r.st[i-1];
});
return r;
}
prefix& operator /= (const prefix& o) { return *this = *this / o; }
friend prefix sqrt(const prefix& a) {
prefix r;
// assert(a.st[1] == 1);
r.st[1] = 1;
T inv_2 = inv(T(2));
r.convolve_helper(r, r, [&](int i, T cur_sum) -> T {
return (a.st[i] - cur_sum) * inv_2 + r.st[i-1];
});
return r;
}
// This computes a pseudo-Euler transform of the sequence.
//
// Formally, given a Dirichlet series
// A = sum a_i / i^s,
// we output the Dirichlet series corresponding to
// B = prod 1 / (1 - a_i i^{-s})
//
// Note: strictly speaking, the standard Euler transform over a generating function should be
// A = sum a_i / i^s -> B = prod 1 / (1 - i^{-s})^a_i
// but our defintion is better suited for totally multiplicative functions,
// and always works over general rings. Also, the two definitions match
// when the a_i are always 0/1.
//
// This runs in $O(n^{2/3})$ time, but requires small inverses (up to 1/120).
friend prefix euler_transform_fraction(prefix a_pref) {
values<layout, T> a(std::move(a_pref));
// assert(a.st[1] == 0);
// Phase 0: stash away values up to the 6th root of N
int x;
for (x = 2; layout.rt / x / x / x > 0; x++) { }
// Phase 1: adjust the values and insert the necessary extra powers
std::array<T, 6> invs{T{}, T(1), inv(T(2)), inv(T(3)), inv(T(4)), inv(T(5))};
for (int i = int(layout.rt); i >= x; i--) {
T v = a.st[i];
int e = 1;
T pv = v;
int64_t pi = i;
while (pi <= layout.N/i) {
e++;
pi *= i;
pv *= v;
a.st[layout.get_value_bucket(pi)] += pv * invs[e];
}
}
// Phase 2: now we take exp of the adjusted version
// In particular, we take e^a = 1 + a + a^2 / 2 + a^3 / 6 + a^4 / 24 + a^5 / 120
prefix v;
for (int i = x; i < layout.len; i++) {
v.st[i] = v.st[i-1] + a.st[i];
}
prefix r = v * v;
for (int i = x; i < layout.len; i++) {
r.st[i] = r.st[i] * invs[5] + v.st[i];
}
r *= v;
for (int i = x; i < layout.len; i++) {
r.st[i] = r.st[i] * invs[4] + v.st[i];
}
r *= v;
for (int i = x; i < layout.len; i++) {
r.st[i] = r.st[i] * invs[3] + v.st[i];
}
r *= v;
for (int i = x; i < layout.len; i++) {
r.st[i] = r.st[i] * invs[2] + v.st[i];
}
for (int i = 1; i < layout.len; i++) {
r.st[i] += T(1);
}
// Phase 3: apply the extra below x
for (x--; x >= 2; x--) {
T ax = a.st[x];
if (ax == 0) continue;
for (int i = x; i < layout.len; i++) {
r.st[i] += r.st[layout.get_value_bucket(layout.get_bucket_bound(i) / x)] * ax;
}
}
return r;
}
// This computes the inverse of the pseudo-Euler transformation. See the
// comment on euler_transform() for more details.
friend prefix inverse_euler_transform_fraction(prefix a) {
values<layout, T> r;
// assert(a.st[1] == 1);
// Phase 1: manually eliminate values up to the 6th root of a
int x;
for (x = 2; layout.rt / x / x / x > 0; x++) {
T v = a.st[x] - T(1);
if (v == 0) continue; // Small optimization, good for prime counting in particular
r.st[x] = v;
for (int i = layout.len - 1; i > x; i--) {
a.st[i] -= a.st[layout.get_value_bucket(layout.get_bucket_bound(i) / x)] * v;
}
a.st[x] = T(1);
}
for (int i = 1; i < x; i++) {
a.st[i] = T();
}
for (int i = x; i < layout.len; i++) {
a.st[i] -= T(1);
}
std::array<T, 6> invs{T{}, T(1), inv(T(2)), inv(T(3)), inv(T(4)), inv(T(5))};
// Phase 2: now we take log of the remaining thing, using just the first few terms.
// In particular, we take log_a = a^5 / 5 - a^4 / 4 + a^3 / 3 - a^2 / 2 + a
prefix log_a;
for (int i = x; i < layout.len; i++) {
log_a.st[i] = a.st[i] * invs[5];
}
log_a *= a;
for (int i = x; i < layout.len; i++) {
log_a.st[i] -= a.st[i] * invs[4];
}
log_a *= a;
for (int i = x; i < layout.len; i++) {
log_a.st[i] += a.st[i] * invs[3];
}
log_a *= a;
for (int i = x; i < layout.len; i++) {
log_a.st[i] -= a.st[i] * invs[2];
}
log_a *= a;
for (int i = x; i < layout.len; i++) {
log_a.st[i] += a.st[i] * invs[1];
}
// Phase 3: correct log_a; we need to get rid of the extra powers.
for (int i = x; i < layout.len; i++) {
r.st[i] = log_a.st[i] - log_a.st[i-1];
}
for (; x <= layout.rt; x++) {
T v = r.st[x];
int e = 1;
T pv = v;
int64_t px = x;
while (px <= layout.N/x) {
e++;
px *= x;
pv *= v;
r.st[layout.get_value_bucket(px)] -= pv * invs[e];
}
}
return prefix(std::move(r));
}
friend prefix euler_transform_binary_indexed_tree(prefix a_pref) {
int x = 2;
while (x <= layout.N / x / x) x++;
prefix r_pref;
for (int i = x; i < layout.len; i++) {
r_pref.st[i] = a_pref.st[i] - a_pref.st[x-1];
}
for (int i = x; i <= layout.rt; i++) {
T vi = a_pref.st[i] - a_pref.st[i-1];
if (vi == 0) continue;
int64_t N_over_i = layout.N / i;
int64_t rt_over_i = layout.rt / i;
int64_t max_z = N_over_i / i;
for (int z = 1; z <= max_z; z++) {
int jlo_idx = i-1;
int jhi_idx = int(z <= rt_over_i ? layout.len - i * z : N_over_i / z);
assert(jlo_idx < jhi_idx);
T v = vi * (a_pref.st[jhi_idx] - a_pref.st[jlo_idx]);
r_pref.st[layout.len-z] += v;
}
}
bit<layout, T> r(std::move(r_pref));
r.increment_bucket_suffix(1, T(1));
for (int i = x-1; i >= 2; i--) {
T cur = a_pref.st[i] - a_pref.st[i-1];
if (cur == 0) continue;
r.sparse_mul_unlimited(i, cur);
}
return prefix(std::move(r));
}
friend prefix inverse_euler_transform_binary_indexed_tree(prefix a_pref) {
// assert(a_pref.st[1] == 1);
bit<layout, T> a_bit(std::move(a_pref));
values<layout, T> r;
// First, use the BIT to clear up to N^1/3
int x;
for (x = 2; x <= layout.N / x / x; x++) {
T cur = a_bit.get_bucket_prefix(x) - T(1);
if (cur == 0) continue;
r.st[x] = cur;
a_bit.sparse_div_unlimited(x, cur);
}
a_pref = prefix<layout, T>(std::move(a_bit));
// Now, a_pref contains terms of the form r[i] or r[i] * r[j], so let's
// subtract out the semiprimes.
// Note that N^1/3 < x <= i <= j, so N/i/j <= N^1/3 < x
for (int i = x; i < layout.len; i++) {
T vi = a_pref.st[i] - a_pref.st[i-1];
r.st[i] = vi;
}
// We want i <= j <= N/i/z
for (int i = x; i <= layout.rt; i++) {
T vi = r.st[i];
if (vi == 0) continue;
int64_t N_over_i = layout.N / i;
int64_t rt_over_i = layout.rt / i;
int64_t max_z = N_over_i / i;
for (int z = 1; z <= max_z; z++) {
int jlo_idx = i-1;
int jhi_idx = int(z <= rt_over_i ? layout.len - i * z : N_over_i / z);
assert(jlo_idx < jhi_idx);
T v = vi * (a_pref.st[jhi_idx] - a_pref.st[jlo_idx]);
r.en[-z] -= v;
if (z > 0) r.en[-(z-1)] += v;
}
}
return prefix(std::move(r));
}
};
// TODO: This will be useful for sparse convolution, which is nice for e.g. exp/log/prime counting
template <div_vector_layout const& layout, typename T> class bit : public div_vector<layout, T>, public vectorspace_mixin<layout, T, bit<layout, T>> {
public:
bit() = default;
template <typename U> explicit bit(bit<layout, U> const& o) {
for (int i = 1; i < layout.len; i++) {
this->st[i] = T(o.st[i]);
}
}
explicit bit(prefix<layout, T> && o) : div_vector<layout, T>(static_cast<div_vector<layout, T>&&>(std::move(o))) {
for (int i = layout.len - 1; i >= 1; i--) {
this->st[i] -= this->st[i & (i-1)];
}
}
explicit bit(prefix<layout, T> const& o) {
for (int i = layout.len - 1; i >= 1; i--) {
this->st[i] = o.st[i] - o.st[i & (i-1)];
}
}
explicit operator prefix<layout, T> () && {
prefix<layout, T> r;
swap(static_cast<div_vector<layout, T>&>(r), static_cast<div_vector<layout, T>&>(*this));
for (int i = 1; i < layout.len; i++) {
r.st[i] += r.st[i & (i-1)];
}
return r;
}
explicit operator prefix<layout, T> () const& {
prefix<layout, T> r;
for (int i = 1; i < layout.len; i++) {
r.st[i] = this->st[i] + r.st[i & (i-1)];
}
return r;
}
T get_bucket_prefix(int a) const {
T r = T();
for (; a > 0; a -= a & -a) {
r += this->st[a];
}
return r;
}
T get_prefix(int64_t v) const {
return get_bucket_prefix(layout.get_value_bucket(v));
}
void increment_bucket_suffix(int a, T d) {
for (; a < layout.len; a += a & -a) {
this->st[a] += d;
}
}
void increment_suffix(int64_t v, T d) const {
return increment_bucket_suffix(layout.get_value_bucket(v), d);
}
// These 4 functions facilitate some simple sparse convolution.
// They each take O(sqrt(N/x) log(N)) time.
// multiply by (1 + w x^s)
void sparse_mul_at_most_one(int64_t x, T w) {
assert(x > 1);
int64_t j = 1;
T cur = get_bucket_prefix(int(j <= layout.rt / x ? layout.len - j * x : layout.N / x / j));
for (; (j+1) <= layout.N / x / (j+1); j++) {
T nxt = get_bucket_prefix(int(j + 1 <= layout.rt / x ? layout.len - (j+1) * x : layout.N / x / (j+1)));
if (cur != nxt) {
increment_bucket_suffix(int(layout.len - j), w * (cur - nxt));
cur = nxt;
}
}
for (int64_t i = layout.N / x / j; i > 0; i--) {
T nxt = get_bucket_prefix(int(i-1));
if (cur != nxt) {
increment_bucket_suffix(int(i <= layout.rt / x ? i * x : layout.len - layout.N / x / i), w * (cur - nxt));
cur = nxt;
}
}
}
// multiply by 1/(1 - w x^s) = 1 + wx^s + w^2 (x^2)^s + ...
void sparse_mul_unlimited(int64_t x, T w) {
assert(x > 1);
T prv = T();
int64_t i;
for (i = 1; i <= layout.N / x / i; i++) {
T cur = get_bucket_prefix(int(i));
if (cur != prv) {
increment_bucket_suffix(int(i <= layout.rt / x ? i * x : layout.len - layout.N / x / i), w * (cur - prv));
prv = cur;
}
}
for (int64_t j = layout.N / x / i; j >= 1; j--) {
T cur = get_bucket_prefix(int(j <= layout.rt / x ? layout.len - j * x : layout.N / x / j));
if (cur != prv) {
increment_bucket_suffix(int(layout.len - j), w * (cur - prv));
prv = cur;
}
}
}
// divide by (1 + w x^s)
void sparse_div_at_most_one(int64_t x, T w) {
return sparse_mul_unlimited(x, -w);
}
// divide by 1/(1 - w x^s) = 1 + wx^s + w^2 (x^2)^s + ...
void sparse_div_unlimited(int64_t x, T w) {
return sparse_mul_at_most_one(x, -w);
}
};
}
#line 2 "src/modnum.hpp"
#line 9 "src/modnum.hpp"
template <typename T> T mod_inv_in_range(T a, T m) {
// assert(0 <= a && a < m);
T x = a, y = m;
// abs coeff of a in x and y (they're always opposite sign)
T vx = 1, vy = 0;
bool swap = false;
while (x) {
T k = y / x;
y %= x;
vy += k * vx;
std::swap(x, y);
std::swap(vx, vy);
swap ^= 1;
}
assert(y == 1);
return swap ? vy : m - vy;
}
template <typename T> struct extended_gcd_result {
T gcd;
T coeff_a, coeff_b;
};
template <typename T> extended_gcd_result<T> extended_gcd(T a, T b) {
T x = a, y = b;
// coeff of a and b in x and y
T ax = 1, ay = 0;
T bx = 0, by = 1;
while (x) {
T k = y / x;
y %= x;
ay -= k * ax;
by -= k * bx;
std::swap(x, y);
std::swap(ax, ay);
std::swap(bx, by);
}
return {y, ay, by};
}
template <typename T> T mod_inv(T a, T m) {
a %= m;
a = a < 0 ? a + m : a;
return mod_inv_in_range(a, m);
}
// Derives the boilerplate operator surface of a number type from its compound
// ops, ==, neg(), and inv().
// Bodies are only instantiated on use, so a type may omit some of the
// underlying pieces if the corresponding derived ops are never called.
template <typename Self>
struct num_ops {
Self operator+ () const { return static_cast<const Self&>(*this); }
Self operator- () const { return static_cast<const Self&>(*this).neg(); }
friend Self operator ++ (Self& a, int) { Self r = a; ++a; return r; }
friend Self operator -- (Self& a, int) { Self r = a; --a; return r; }
friend Self operator + (const Self& a, const Self& b) { return Self(a) += b; }
friend Self operator - (const Self& a, const Self& b) { return Self(a) -= b; }
friend Self operator * (const Self& a, const Self& b) { return Self(a) *= b; }
friend Self operator / (const Self& a, const Self& b) { return Self(a) /= b; }
friend bool operator != (const Self& a, const Self& b) { return !(a == b); }
friend Self neg(const Self& a) { return a.neg(); }
friend Self inv(const Self& a) { return a.inv(); }
};
// Storage and arithmetic for numbers mod Self::MOD, as a reduced
// representative v in [0, MOD) of unsigned type V.
// The type provides static MOD (of type V), reduce (value -> representative),
// and *=;
// everything else is derived here, valid for any MOD up to V's full range
// (sums and differences are tracked mod 2^bits, so no headroom is needed).
// Hooks may be overridden in the type's own body (e.g. a faster += / -=).
template <typename Self, typename V>
struct mod_ops : num_ops<Self> {
static_assert(std::unsigned_integral<V>);
V v;
struct is_reduced_tag {};
mod_ops() : v(0) {}
mod_ops(V v_, is_reduced_tag) : v(v_) { assert(v < Self::MOD); }
template <std::integral I> mod_ops(I x) : v(Self::reduce(x)) {}
static Self from_reduced(V v) { return Self(v, is_reduced_tag{}); }
// A negative value reduces via its nonnegative complement: x = -1 - ~x.
static V reduce(std::signed_integral auto x) {
using U = std::make_unsigned_t<decltype(x)>;
return x < 0 ? V(Self::MOD - 1 - Self::reduce(U(~x))) : Self::reduce(U(x));
}
explicit operator V() const { return v; }
std::make_signed_t<V> balanced() const {
return std::make_signed_t<V>(Self::MOD-v > v ? v : v - Self::MOD);
}
friend bool operator == (const Self& a, const Self& b) { return a.v == b.v; }
friend std::ostream& operator << (std::ostream& out, const Self& n) { return out << n.v; }
friend std::istream& operator >> (std::istream& in, Self& n) { int64_t v_; in >> v_; n = Self(v_); return in; }
Self& operator ++ () {
++v;
if (v == Self::MOD) v = 0;
return self();
}
Self& operator -- () {
if (v == 0) v = Self::MOD;
--v;
return self();
}
Self& operator += (const Self& o) { v = Self::sub_mod_raw(v, Self::MOD - o.v); return self(); }
Self& operator -= (const Self& o) { v = Self::sub_mod_raw(v, o.v); return self(); }
Self& operator /= (const Self& o) { return self() *= o.inv(); }
// Returns a - b mod MOD, for b in [0, MOD]; wraparound detects the underflow.
static V sub_mod_raw(V a, V b) { return a < b ? a - b + Self::MOD : a - b; }
Self neg() const { return from_reduced(v ? Self::MOD - v : 0); }
Self inv() const { return from_reduced(mod_inv_in_range(v, Self::MOD)); }
private:
Self& self() { return static_cast<Self&>(*this); }
};
template <auto MOD_> struct modnum : mod_ops<modnum<MOD_>, std::make_unsigned_t<decltype(MOD_)>> {
using Self = modnum;
static_assert(MOD_ > 0, "MOD must be positive");
using V = std::make_unsigned_t<decltype(MOD_)>;
static constexpr V MOD = V(MOD_);
using base = mod_ops<modnum, V>;
using base::base;
using base::v;
using base::reduce;
static V reduce(std::unsigned_integral auto x) { return V(x % MOD); }
explicit operator std::make_signed_t<V>() const
requires (MOD <= V(std::numeric_limits<std::make_signed_t<V>>::max()))
{
return std::make_signed_t<V>(v);
}
Self& operator *= (const Self& o) {
if constexpr (sizeof(V) <= 4) v = V(uint64_t(v) * o.v % MOD);
else v = V(__uint128_t(v) * o.v % MOD);
return *this;
}
};
struct mod_goldilocks : mod_ops<mod_goldilocks, uint64_t> {
using Self = mod_goldilocks;
static constexpr uint64_t MOD = 0xffffffff00000001ull;
static constexpr uint64_t EPS = -MOD;
// We have 2^32 is a primitive 6th root of unity.
// Note that omega_8 + omega_8^7 == 2^24 - 2^72 == sqrt(2)
// We'll pick the root so that 2^24 - 2^72 is our primitive 384th root of unity.
static constexpr uint64_t PRIMITIVE_ROOT = 2717;
using base = mod_ops<mod_goldilocks, uint64_t>;
using base::base;
using base::reduce;
mod_goldilocks() = default;
mod_goldilocks(__int128_t a) : base(a < 0 ? uint64_t(MOD - 1 - __uint128_t(~a) % MOD) : uint64_t(__uint128_t(a) % MOD), is_reduced_tag{}) {}
mod_goldilocks(__uint128_t a) : base(uint64_t(a % MOD), is_reduced_tag{}) {}
// Avoids the division: any uint64_t is within MOD of reduced.
static uint64_t reduce(std::unsigned_integral auto x) {
static_assert(sizeof(x) <= 8);
uint64_t a = x;
return a >= MOD ? a - MOD : a;
}
// returns a-b, assuming -MOD <= a-b, e.g. b <= MOD
static uint64_t sub_mod_raw(uint64_t a, uint64_t b) {
#if defined(__x86_64__)
// TODO: We could try to write this using intrinsics, but GCC sometimes produces the wrong code.
uint64_t res_wrapped = a;
uint64_t adjustment = b;
asm (
// AT&T syntax: SRC DST
"sub %[y], %[x]\n\t"
// Trick from plonky2 implementation:
// After sub, flag CF is set iff we underflowed. We want to correct by EPS == 2^32 - 1 iff C is set.
// sbb (subtract with borrow) computes DST <- DST - SRC - CF
// Thus, we can use the 32-bit form of sbb on a dummy register to load CF ? EPS : 0.
// Here, we'll just reuse the original register holding b.
"sbb %k[y], %k[y]\n\t"
: [x] "+r"(res_wrapped),
[y] "+r"(adjustment)
:
: "cc"
);
#else
uint64_t res_wrapped = a - b;
uint64_t adjustment = (res_wrapped > a) ? EPS : 0;
#endif
return res_wrapped - adjustment;
}
// Reduce lo + 2^64 * mi + 2^96 * hi, where hi <= MOD
static uint64_t reduce_u160_raw(uint64_t lo, uint32_t mi, uint64_t hi) {
// result = lo - hi + EPS * mi
// 0 <= lo <= 2^64 - 1 = MOD + EPS - 1
// 0 <= EPS * mi <= (2^32 - 1) * EPS = MOD - 1 - EPS
// 0 <= hi <= MOD
// -MOD <= lo - hi + EPS * mi <= 2*MOD-2
// so we do have some leeway
return sub_mod_raw(sub_mod_raw(lo, hi), MOD-(uint64_t(mi)<<32)+mi);
}
static uint64_t reduce_u128_raw(__uint128_t v) {
uint64_t hi = uint64_t(v >> 64);
uint64_t lo = uint64_t(v);
uint32_t hi_hi = uint32_t(hi >> 32);
uint32_t hi_lo = uint32_t(hi);
return reduce_u160_raw(lo, hi_lo, hi_hi);
}
Self& operator *= (Self o) {
v = reduce_u128_raw(__uint128_t(v) * __uint128_t(o.v));
return *this;
}
};
template <typename T> T power(T a, long long b) {
assert(b >= 0);
T r = 1; while (b) { if (b & 1) r *= a; b >>= 1; a *= a; } return r;
}
template <typename U, typename V> struct pairnum : num_ops<pairnum<U, V>> {
using Self = pairnum;
U u;
V v;
pairnum() : u(0), v(0) {}
pairnum(long long val) : u(val), v(val) {}
pairnum(const U& u_, const V& v_) : u(u_), v(v_) {}
friend std::ostream& operator << (std::ostream& out, const Self& n) { return out << '(' << n.u << ',' << ' ' << n.v << ')'; }
friend std::istream& operator >> (std::istream& in, Self& n) { long long val; in >> val; n = Self(val); return in; }
friend bool operator == (const Self& a, const Self& b) { return a.u == b.u && a.v == b.v; }
Self inv() const {
return Self(u.inv(), v.inv());
}
Self neg() const {
return Self(u.neg(), v.neg());
}
Self& operator ++ () {
++u, ++v;
return *this;
}
Self& operator -- () {
--u, --v;
return *this;
}
Self& operator += (const Self& o) {
u += o.u;
v += o.v;
return *this;
}
Self& operator -= (const Self& o) {
u -= o.u;
v -= o.v;
return *this;
}
Self& operator *= (const Self& o) {
u *= o.u;
v *= o.v;
return *this;
}
Self& operator /= (const Self& o) {
u /= o.u;
v /= o.v;
return *this;
}
};
template <typename tag> struct dynamic_modnum : mod_ops<dynamic_modnum<tag>, uint32_t> {
using Self = dynamic_modnum;
private:
inline static uint32_t MOD_ = 0;
inline static uint64_t BARRETT_M = 0;
public:
// Make only the const-reference public, to force the use of set_mod
static constexpr uint32_t const& MOD = MOD_;
using base = mod_ops<dynamic_modnum, uint32_t>;
using base::base;
using base::v;
using base::reduce;
// Barret reduction taken from KACTL:
/**
* Author: Simon Lindholm
* Date: 2020-05-30
* License: CC0
* Source: https://en.wikipedia.org/wiki/Barrett_reduction
* Description: Compute $a \% b$ about 5 times faster than usual, where $b$ is constant but not known at compile time.
* Returns a value congruent to $a \pmod b$ in the range $[0, 2b)$.
* Status: proven correct, stress-tested
* Measured as having 4 times lower latency, and 8 times higher throughput, see stress-test.
* Details:
* More precisely, it can be proven that the result equals 0 only if $a = 0$,
* and otherwise lies in $[1, (1 + a/2^64) * b)$.
*/
static void set_mod(int mod) {
assert(mod > 0);
MOD_ = uint32_t(mod);
BARRETT_M = (uint64_t(-1) / MOD);
}
static uint32_t barrett_reduce_partial(uint64_t a) {
return uint32_t(a - uint64_t((__uint128_t(BARRETT_M) * a) >> 64) * MOD);
}
static uint32_t barrett_reduce(uint64_t a) {
int32_t res = int32_t(barrett_reduce_partial(a) - MOD);
return uint32_t((res < 0) ? res + int32_t(MOD) : res);
}
struct mod_reader {
friend std::istream& operator >> (std::istream& i, mod_reader) {
int mod; i >> mod;
Self::set_mod(mod);
return i;
}
};
static mod_reader MOD_READER() {
return mod_reader();
}
static uint32_t reduce(std::unsigned_integral auto x) {
static_assert(sizeof(x) <= 8);
return barrett_reduce(x);
}
explicit operator int() const { return int(v); }
Self& operator *= (const Self& o) {
v = barrett_reduce(uint64_t(v) * o.v);
return *this;
}
};
template <typename T> struct mod_constraint {
T v, mod;
friend mod_constraint operator & (mod_constraint a, mod_constraint b) {
if (a.mod < b.mod) std::swap(a, b);
if (b.mod == 1) return a;
extended_gcd_result<T> egcd = extended_gcd<T>(a.mod, b.mod);
assert(a.v % egcd.gcd == b.v % egcd.gcd);
T extra = b.v - a.v % b.mod;
extra /= egcd.gcd;
extra *= egcd.coeff_a;
extra %= b.mod / egcd.gcd;
extra += (extra < 0) ? b.mod / egcd.gcd : 0;
return mod_constraint{
a.v + extra * a.mod,
a.mod * (b.mod / egcd.gcd)
};
}
};
#line 8 "verify/dirichlet_inverse_and_prefix_sums.test.cpp"
int main() {
std::ios_base::sync_with_stdio(false), std::cin.tie(nullptr);
int T; std::cin >> T;
while (T--) {
int64_t N; std::cin >> N;
static dirichlet_series::div_vector_layout layout;
layout = N;
using num = modnum<998244353>;
using ds_prefix = dirichlet_series::prefix<layout, num>;
ds_prefix F;
for (int i = 1; i < layout.len; i++) std::cin >> F.st[i];
ds_prefix delta([&](int64_t) { return 1; });
ds_prefix H = delta / F;
for (int i = 1; i < layout.len; i++) std::cout << H.st[i] << " \n"[i+1==layout.len];
}
return 0;
}
// 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>
// src/dirichlet_series.hpp
namespace dirichlet_series{
inline int inv(int v){
assert(v==1);
return 1;
}
inline int64_t inv(int64_t v){
assert(v==1);
return 1;
}
constexpr int64_t floor_sqrt(int64_t N){
assert(N>=0);
if(N==0)return 0;
int64_t a=N;
while(true){
int64_t b=N/a;
assert(a>=b);
if(a-b<=1)return b;
a=(a+b+1)>>1;
}
}
class div_vector_layout{
public:
int64_t N;
int64_t rt=floor_sqrt(N);
int len=int(2*rt+(rt*(rt+1)<=N));
constexpr div_vector_layout(int64_t N_=1):N(N_){}
constexpr int get_value_bucket(int64_t a)const{
return a<=rt?int(a):len-int(N/a);
}
constexpr int64_t get_bucket_bound(int i)const{
return i<=rt?i:N/(len-i);
}
};
template<const div_vector_layout&layout,typename T>class div_vector{
public:
T*st=new T[layout.len+1]{};
T*en=st+layout.len;
div_vector()=default;
div_vector(div_vector const&o){
std::copy(o.st,o.en,st);
}
div_vector&operator=(div_vector const&o){
std::copy(o.st,o.en,st);
return*this;
}
friend void swap(div_vector&a,div_vector&b){
std::swap(a.st,b.st);
std::swap(a.en,b.en);
}
div_vector(div_vector&&o):st(nullptr),en(nullptr){
swap(*this,o);
}
div_vector&operator=(div_vector&&o){
swap(*this,o);
return*this;
}
~div_vector(){delete[]st;}
T&operator[](int64_t v){return st[layout.get_value_bucket(v)];}
T&operator[](int64_t v)const{return st[layout.get_value_bucket(v)];}
};
template<div_vector_layout const&layout,typename T,typename Derived>class vectorspace_mixin{
private:
Derived&underlying(){
return static_cast<Derived&>(*this);
}
Derived const&underlying()const{
return static_cast<Derived const&>(*this);
}
public:
friend Derived operator+(Derived&&a){
for(int64_t i=1;i<layout.len;i++){
a.st[i]=+a.st[i];
}
return a;
}
friend Derived operator+(Derived const&a){return+Derived(a);}
friend Derived operator-(Derived&&a){
for(int64_t i=1;i<layout.len;i++){
a.st[i]=-a.st[i];
}
return a;
}
friend Derived operator-(Derived const&a){return-Derived(a);}
Derived&operator+=(Derived const&o){
for(int64_t i=1;i<layout.len;i++){
underlying().st[i]+=o.st[i];
}
return underlying();
}
friend Derived operator+(Derived&&a,Derived const&b){return a+=b;}
friend Derived operator+(Derived const&a,Derived&&b){
for(int64_t i=1;i<layout.len;i++){
b.st[i]=a.st[i]+b.st[i];
}
return b;
}
friend Derived operator+(Derived&&a,Derived&&b){return std::move(a)+b;}
friend Derived operator+(Derived const&a,Derived const&b){return Derived(a)+b;}
template<typename F>Derived&operator+=(F f){
for(int64_t i=1;i<layout.len;i++){
underlying().st[i]+=f(layout.get_bucket_bound(i));
}
return underlying();
}
Derived&operator-=(Derived const&o){
for(int64_t i=1;i<layout.len;i++){
underlying().st[i]-=o.st[i];
}
return underlying();
}
friend Derived operator-(Derived&&a,Derived const&b){return a-=b;}
friend Derived operator-(Derived const&a,Derived&&b){
for(int64_t i=1;i<layout.len;i++){
b.st[i]=a.st[i]-b.st[i];
}
return b;
}
friend Derived operator-(Derived&&a,Derived&&b){return std::move(a)-b;}
friend Derived operator-(Derived const&a,Derived const&b){return Derived(a)-b;}
template<typename F>Derived&operator-=(F f){
for(int64_t i=1;i<layout.len;i++){
underlying().st[i]-=f(layout.get_bucket_bound(i));
}
return underlying();
}
Derived&operator*=(T const&t){
for(int64_t i=1;i<layout.len;i++){
underlying().st[i]*=t;
}
return underlying();
}
friend Derived operator*(Derived&&a,T const&t){return a*=t;}
friend Derived operator*(Derived const&a,T const&t){return Derived(a)*t;}
friend Derived operator*(T const&t,Derived&&a){
for(int64_t i=1;i<layout.len;i++){
a.st[i]=t*a.st[i];
}
return a;
}
friend Derived operator*(T const&t,Derived const&a){return t*Derived(a);}
Derived&operator/=(T const&t){
for(int64_t i=1;i<layout.len;i++){
underlying().st[i]/=t;
}
return underlying();
}
friend Derived operator/(Derived&&a,T const&t){return a/=t;}
friend Derived operator/(Derived const&a,T const&t){return Derived(a)/t;}
};
template<div_vector_layout const&layout,typename T>class values;
template<div_vector_layout const&layout,typename T>class prefix;
template<div_vector_layout const&layout,typename T>class bit;
template<div_vector_layout const&layout,typename T>class values:public div_vector<layout,T>,public vectorspace_mixin<layout,T,values<layout,T>>{
public:
values()=default;
template<typename F,std::enable_if_t<std::is_invocable_r_v<T,F,int64_t,int64_t>,bool> =true>
values(F f){
for(int i=1;i<layout.len;i++){
this->st[i]=f(layout.get_bucket_bound(i-1),layout.get_bucket_bound(i));
}
}
template<typename U>explicit values(values<layout,U>const&o){
for(int i=1;i<layout.len;i++){
this->st[i]=T(o.st[i]);
}
}
explicit values(prefix<layout,T>&&o):div_vector<layout,T>(static_cast<div_vector<layout,T>&&>(std::move(o))){
for(int i=layout.len-1;i>1;i--){
this->st[i]-=this->st[i-1];
}
}
explicit values(prefix<layout,T>const&o){
for(int i=layout.len-1;i>1;i--){
this->st[i]=o.st[i]-o.st[i-1];
}
this->st[1]=o.st[1];
}
};
template<div_vector_layout const&layout,typename T>class prefix:public div_vector<layout,T>,public vectorspace_mixin<layout,T,prefix<layout,T>>{
public:
prefix()=default;
template<typename F,std::enable_if_t<std::is_invocable_r_v<T,F,int64_t>,bool> =true>
prefix(F f){
for(int i=1;i<layout.len;i++){
this->st[i]=f(layout.get_bucket_bound(i));
}
}
template<typename U>explicit prefix(prefix<layout,U>const&o){
for(int i=1;i<layout.len;i++){
this->st[i]=T(o.st[i]);
}
}
explicit prefix(values<layout,T>&&o):div_vector<layout,T>(static_cast<div_vector<layout,T>&&>(std::move(o))){
for(int i=2;i<layout.len;i++){
this->st[i]+=this->st[i-1];
}
}
explicit prefix(values<layout,T>const&o){
T pref=this->st[1]=o.st[1];
for(int i=2;i<layout.len;i++){
this->st[i]=(pref+=o.st[i]);
}
}
private:
template<typename F>
void convolve_helper(prefix const&a,prefix const&b,F f){
T cur_sum=a.st[1]*b.st[1];
for(int i=2;i<layout.len;i++){
cur_sum+=this->st[i];
if(i>=layout.len-layout.rt){
int z=int(layout.len-i);
assert(z<=layout.rt);
int64_t rt_over_z=layout.rt/z;
int64_t N_over_z=layout.N/z;
int64_t x_max=N_over_z/(z+1);
T tot_val=T();
for(int64_t x=2;x*(x+1)<=N_over_z&&x<=x_max;x++){
int ylo_idx=std::max(int(x),z);
int yhi_idx=int(x<=rt_over_z?layout.len-x*z:N_over_z/x);
assert(ylo_idx<yhi_idx);
T ax=a.st[x]-a.st[x-1];
T bx=b.st[x]-b.st[x-1];
T ay=a.st[yhi_idx]-a.st[ylo_idx];
T by=b.st[yhi_idx]-b.st[ylo_idx];
T v=ax*by+ay*bx;
tot_val+=v;
}
cur_sum+=tot_val;
if(i+1<layout.len){
this->st[i+1]-=tot_val;
}
}
this->st[i]=f(i,cur_sum);
T ai=a.st[i]-a.st[i-1];
T bi=b.st[i]-b.st[i-1];
cur_sum+=ai*b.st[1]+a.st[1]*bi;
if(i<=layout.rt){
int64_t rt_over_i=layout.rt/i;
int64_t N_over_i=layout.N/i;
int x_max=int(std::min<int64_t>(N_over_i/i,i-1));
T tot_sub=T();
for(int x=2;x<=x_max;x++){
T v;
v=ai*(b.st[x]-b.st[x-1])+(a.st[x]-a.st[x-1])*bi;
int zlo_idx=int(x<=rt_over_i?x*i:layout.len-(N_over_i/x));
this->st[zlo_idx]+=v;
tot_sub+=v;
}
this->en[-(i-1)]-=tot_sub;
{
int zlo_idx=int(i<=rt_over_i?i*i:layout.len-(N_over_i/i));
this->st[zlo_idx]+=ai*bi;
}
}
}
}
public:
friend prefix operator*(prefix const&a,prefix const&b){
prefix r;
r.st[1]=a.st[1]*b.st[1];
r.convolve_helper(a,b,[&](int i,T cur_sum)->T{
return cur_sum+(a.st[i]-a.st[i-1])*b.st[1]+a.st[1]*(b.st[i]-b.st[i-1]);
});
return r;
}
prefix&operator*=(const prefix&o){return*this=*this*o;}
friend T get_conv_N(prefix const&a,prefix const&b){
T ans=a.st[1]*b.en[-1];
for(int i=2;i<=layout.len;i++){
ans+=(a.st[i]-a.st[i-1])*b.en[-i];
}
return ans;
}
friend prefix operator/(prefix const&a,prefix const&b){
prefix r;
T inv_b1=inv(b.st[1]);
r.st[1]=a.st[1]*inv_b1;
r.convolve_helper(r,b,[&](int i,T cur_sum)->T{
return(a.st[i]-(cur_sum+r.st[1]*(b.st[i]-b.st[i-1])))*inv_b1+r.st[i-1];
});
return r;
}
prefix&operator/=(const prefix&o){return*this=*this/o;}
friend prefix sqrt(const prefix&a){
prefix r;
r.st[1]=1;
T inv_2=inv(T(2));
r.convolve_helper(r,r,[&](int i,T cur_sum)->T{
return(a.st[i]-cur_sum)*inv_2+r.st[i-1];
});
return r;
}
friend prefix euler_transform_fraction(prefix a_pref){
values<layout,T>a(std::move(a_pref));
int x;
for(x=2;layout.rt/x/x/x>0;x++){}
std::array<T,6>invs{T{},T(1),inv(T(2)),inv(T(3)),inv(T(4)),inv(T(5))};
for(int i=int(layout.rt);i>=x;i--){
T v=a.st[i];
int e=1;
T pv=v;
int64_t pi=i;
while(pi<=layout.N/i){
e++;
pi*=i;
pv*=v;
a.st[layout.get_value_bucket(pi)]+=pv*invs[e];
}
}
prefix v;
for(int i=x;i<layout.len;i++){
v.st[i]=v.st[i-1]+a.st[i];
}
prefix r=v*v;
for(int i=x;i<layout.len;i++){
r.st[i]=r.st[i]*invs[5]+v.st[i];
}
r*=v;
for(int i=x;i<layout.len;i++){
r.st[i]=r.st[i]*invs[4]+v.st[i];
}
r*=v;
for(int i=x;i<layout.len;i++){
r.st[i]=r.st[i]*invs[3]+v.st[i];
}
r*=v;
for(int i=x;i<layout.len;i++){
r.st[i]=r.st[i]*invs[2]+v.st[i];
}
for(int i=1;i<layout.len;i++){
r.st[i]+=T(1);
}
for(x--;x>=2;x--){
T ax=a.st[x];
if(ax==0)continue;
for(int i=x;i<layout.len;i++){
r.st[i]+=r.st[layout.get_value_bucket(layout.get_bucket_bound(i)/x)]*ax;
}
}
return r;
}
friend prefix inverse_euler_transform_fraction(prefix a){
values<layout,T>r;
int x;
for(x=2;layout.rt/x/x/x>0;x++){
T v=a.st[x]-T(1);
if(v==0)continue;
r.st[x]=v;
for(int i=layout.len-1;i>x;i--){
a.st[i]-=a.st[layout.get_value_bucket(layout.get_bucket_bound(i)/x)]*v;
}
a.st[x]=T(1);
}
for(int i=1;i<x;i++){
a.st[i]=T();
}
for(int i=x;i<layout.len;i++){
a.st[i]-=T(1);
}
std::array<T,6>invs{T{},T(1),inv(T(2)),inv(T(3)),inv(T(4)),inv(T(5))};
prefix log_a;
for(int i=x;i<layout.len;i++){
log_a.st[i]=a.st[i]*invs[5];
}
log_a*=a;
for(int i=x;i<layout.len;i++){
log_a.st[i]-=a.st[i]*invs[4];
}
log_a*=a;
for(int i=x;i<layout.len;i++){
log_a.st[i]+=a.st[i]*invs[3];
}
log_a*=a;
for(int i=x;i<layout.len;i++){
log_a.st[i]-=a.st[i]*invs[2];
}
log_a*=a;
for(int i=x;i<layout.len;i++){
log_a.st[i]+=a.st[i]*invs[1];
}
for(int i=x;i<layout.len;i++){
r.st[i]=log_a.st[i]-log_a.st[i-1];
}
for(;x<=layout.rt;x++){
T v=r.st[x];
int e=1;
T pv=v;
int64_t px=x;
while(px<=layout.N/x){
e++;
px*=x;
pv*=v;
r.st[layout.get_value_bucket(px)]-=pv*invs[e];
}
}
return prefix(std::move(r));
}
friend prefix euler_transform_binary_indexed_tree(prefix a_pref){
int x=2;
while(x<=layout.N/x/x)x++;
prefix r_pref;
for(int i=x;i<layout.len;i++){
r_pref.st[i]=a_pref.st[i]-a_pref.st[x-1];
}
for(int i=x;i<=layout.rt;i++){
T vi=a_pref.st[i]-a_pref.st[i-1];
if(vi==0)continue;
int64_t N_over_i=layout.N/i;
int64_t rt_over_i=layout.rt/i;
int64_t max_z=N_over_i/i;
for(int z=1;z<=max_z;z++){
int jlo_idx=i-1;
int jhi_idx=int(z<=rt_over_i?layout.len-i*z:N_over_i/z);
assert(jlo_idx<jhi_idx);
T v=vi*(a_pref.st[jhi_idx]-a_pref.st[jlo_idx]);
r_pref.st[layout.len-z]+=v;
}
}
bit<layout,T>r(std::move(r_pref));
r.increment_bucket_suffix(1,T(1));
for(int i=x-1;i>=2;i--){
T cur=a_pref.st[i]-a_pref.st[i-1];
if(cur==0)continue;
r.sparse_mul_unlimited(i,cur);
}
return prefix(std::move(r));
}
friend prefix inverse_euler_transform_binary_indexed_tree(prefix a_pref){
bit<layout,T>a_bit(std::move(a_pref));
values<layout,T>r;
int x;
for(x=2;x<=layout.N/x/x;x++){
T cur=a_bit.get_bucket_prefix(x)-T(1);
if(cur==0)continue;
r.st[x]=cur;
a_bit.sparse_div_unlimited(x,cur);
}
a_pref=prefix<layout,T>(std::move(a_bit));
for(int i=x;i<layout.len;i++){
T vi=a_pref.st[i]-a_pref.st[i-1];
r.st[i]=vi;
}
for(int i=x;i<=layout.rt;i++){
T vi=r.st[i];
if(vi==0)continue;
int64_t N_over_i=layout.N/i;
int64_t rt_over_i=layout.rt/i;
int64_t max_z=N_over_i/i;
for(int z=1;z<=max_z;z++){
int jlo_idx=i-1;
int jhi_idx=int(z<=rt_over_i?layout.len-i*z:N_over_i/z);
assert(jlo_idx<jhi_idx);
T v=vi*(a_pref.st[jhi_idx]-a_pref.st[jlo_idx]);
r.en[-z]-=v;
if(z>0)r.en[-(z-1)]+=v;
}
}
return prefix(std::move(r));
}
};
template<div_vector_layout const&layout,typename T>class bit:public div_vector<layout,T>,public vectorspace_mixin<layout,T,bit<layout,T>>{
public:
bit()=default;
template<typename U>explicit bit(bit<layout,U>const&o){
for(int i=1;i<layout.len;i++){
this->st[i]=T(o.st[i]);
}
}
explicit bit(prefix<layout,T>&&o):div_vector<layout,T>(static_cast<div_vector<layout,T>&&>(std::move(o))){
for(int i=layout.len-1;i>=1;i--){
this->st[i]-=this->st[i&(i-1)];
}
}
explicit bit(prefix<layout,T>const&o){
for(int i=layout.len-1;i>=1;i--){
this->st[i]=o.st[i]-o.st[i&(i-1)];
}
}
explicit operator prefix<layout,T>()&&{
prefix<layout,T>r;
swap(static_cast<div_vector<layout,T>&>(r),static_cast<div_vector<layout,T>&>(*this));
for(int i=1;i<layout.len;i++){
r.st[i]+=r.st[i&(i-1)];
}
return r;
}
explicit operator prefix<layout,T>()const&{
prefix<layout,T>r;
for(int i=1;i<layout.len;i++){
r.st[i]=this->st[i]+r.st[i&(i-1)];
}
return r;
}
T get_bucket_prefix(int a)const{
T r=T();
for(;a>0;a-=a&-a){
r+=this->st[a];
}
return r;
}
T get_prefix(int64_t v)const{
return get_bucket_prefix(layout.get_value_bucket(v));
}
void increment_bucket_suffix(int a,T d){
for(;a<layout.len;a+=a&-a){
this->st[a]+=d;
}
}
void increment_suffix(int64_t v,T d)const{
return increment_bucket_suffix(layout.get_value_bucket(v),d);
}
void sparse_mul_at_most_one(int64_t x,T w){
assert(x>1);
int64_t j=1;
T cur=get_bucket_prefix(int(j<=layout.rt/x?layout.len-j*x:layout.N/x/j));
for(;(j+1)<=layout.N/x/(j+1);j++){
T nxt=get_bucket_prefix(int(j+1<=layout.rt/x?layout.len-(j+1)*x:layout.N/x/(j+1)));
if(cur!=nxt){
increment_bucket_suffix(int(layout.len-j),w*(cur-nxt));
cur=nxt;
}
}
for(int64_t i=layout.N/x/j;i>0;i--){
T nxt=get_bucket_prefix(int(i-1));
if(cur!=nxt){
increment_bucket_suffix(int(i<=layout.rt/x?i*x:layout.len-layout.N/x/i),w*(cur-nxt));
cur=nxt;
}
}
}
void sparse_mul_unlimited(int64_t x,T w){
assert(x>1);
T prv=T();
int64_t i;
for(i=1;i<=layout.N/x/i;i++){
T cur=get_bucket_prefix(int(i));
if(cur!=prv){
increment_bucket_suffix(int(i<=layout.rt/x?i*x:layout.len-layout.N/x/i),w*(cur-prv));
prv=cur;
}
}
for(int64_t j=layout.N/x/i;j>=1;j--){
T cur=get_bucket_prefix(int(j<=layout.rt/x?layout.len-j*x:layout.N/x/j));
if(cur!=prv){
increment_bucket_suffix(int(layout.len-j),w*(cur-prv));
prv=cur;
}
}
}
void sparse_div_at_most_one(int64_t x,T w){
return sparse_mul_unlimited(x,-w);
}
void sparse_div_unlimited(int64_t x,T w){
return sparse_mul_at_most_one(x,-w);
}
};
}
// src/modnum.hpp
template<typename T>T mod_inv_in_range(T a,T m){
T x=a,y=m;
T vx=1,vy=0;
bool swap=false;
while(x){
T k=y/x;
y%=x;
vy+=k*vx;
std::swap(x,y);
std::swap(vx,vy);
swap^=1;
}
assert(y==1);
return swap?vy:m-vy;
}
template<typename T>struct extended_gcd_result{
T gcd;
T coeff_a,coeff_b;
};
template<typename T>extended_gcd_result<T>extended_gcd(T a,T b){
T x=a,y=b;
T ax=1,ay=0;
T bx=0,by=1;
while(x){
T k=y/x;
y%=x;
ay-=k*ax;
by-=k*bx;
std::swap(x,y);
std::swap(ax,ay);
std::swap(bx,by);
}
return{y,ay,by};
}
template<typename T>T mod_inv(T a,T m){
a%=m;
a=a<0?a+m:a;
return mod_inv_in_range(a,m);
}
template<typename Self>
struct num_ops{
Self operator+()const{return static_cast<const Self&>(*this);}
Self operator-()const{return static_cast<const Self&>(*this).neg();}
friend Self operator++(Self&a,int){Self r=a;++a;return r;}
friend Self operator--(Self&a,int){Self r=a;--a;return r;}
friend Self operator+(const Self&a,const Self&b){return Self(a)+=b;}
friend Self operator-(const Self&a,const Self&b){return Self(a)-=b;}
friend Self operator*(const Self&a,const Self&b){return Self(a)*=b;}
friend Self operator/(const Self&a,const Self&b){return Self(a)/=b;}
friend bool operator!=(const Self&a,const Self&b){return!(a==b);}
friend Self neg(const Self&a){return a.neg();}
friend Self inv(const Self&a){return a.inv();}
};
template<typename Self,typename V>
struct mod_ops:num_ops<Self>{
static_assert(std::unsigned_integral<V>);
V v;
struct is_reduced_tag{};
mod_ops():v(0){}
mod_ops(V v_,is_reduced_tag):v(v_){assert(v<Self::MOD);}
template<std::integral I>mod_ops(I x):v(Self::reduce(x)){}
static Self from_reduced(V v){return Self(v,is_reduced_tag{});}
static V reduce(std::signed_integral auto x){
using U=std::make_unsigned_t<decltype(x)>;
return x<0?V(Self::MOD-1-Self::reduce(U(~x))):Self::reduce(U(x));
}
explicit operator V()const{return v;}
std::make_signed_t<V>balanced()const{
return std::make_signed_t<V>(Self::MOD-v>v?v:v-Self::MOD);
}
friend bool operator==(const Self&a,const Self&b){return a.v==b.v;}
friend std::ostream&operator<<(std::ostream&out,const Self&n){return out<<n.v;}
friend std::istream&operator>>(std::istream&in,Self&n){int64_t v_;in>>v_;n=Self(v_);return in;}
Self&operator++(){
++v;
if(v==Self::MOD)v=0;
return self();
}
Self&operator--(){
if(v==0)v=Self::MOD;
--v;
return self();
}
Self&operator+=(const Self&o){v=Self::sub_mod_raw(v,Self::MOD-o.v);return self();}
Self&operator-=(const Self&o){v=Self::sub_mod_raw(v,o.v);return self();}
Self&operator/=(const Self&o){return self()*=o.inv();}
static V sub_mod_raw(V a,V b){return a<b?a-b+Self::MOD:a-b;}
Self neg()const{return from_reduced(v?Self::MOD-v:0);}
Self inv()const{return from_reduced(mod_inv_in_range(v,Self::MOD));}
private:
Self&self(){return static_cast<Self&>(*this);}
};
template<auto MOD_>struct modnum:mod_ops<modnum<MOD_>,std::make_unsigned_t<decltype(MOD_)>>{
using Self=modnum;
static_assert(MOD_>0,"MOD must be positive");
using V=std::make_unsigned_t<decltype(MOD_)>;
static constexpr V MOD=V(MOD_);
using base=mod_ops<modnum,V>;
using base::base;
using base::v;
using base::reduce;
static V reduce(std::unsigned_integral auto x){return V(x%MOD);}
explicit operator std::make_signed_t<V>()const
requires(MOD<=V(std::numeric_limits<std::make_signed_t<V>>::max()))
{
return std::make_signed_t<V>(v);
}
Self&operator*=(const Self&o){
if constexpr(sizeof(V)<=4)v=V(uint64_t(v)*o.v%MOD);
else v=V(__uint128_t(v)*o.v%MOD);
return*this;
}
};
struct mod_goldilocks:mod_ops<mod_goldilocks,uint64_t>{
using Self=mod_goldilocks;
static constexpr uint64_t MOD=0xffffffff00000001ull;
static constexpr uint64_t EPS=-MOD;
static constexpr uint64_t PRIMITIVE_ROOT=2717;
using base=mod_ops<mod_goldilocks,uint64_t>;
using base::base;
using base::reduce;
mod_goldilocks()=default;
mod_goldilocks(__int128_t a):base(a<0?uint64_t(MOD-1-__uint128_t(~a)%MOD):uint64_t(__uint128_t(a)%MOD),is_reduced_tag{}){}
mod_goldilocks(__uint128_t a):base(uint64_t(a%MOD),is_reduced_tag{}){}
static uint64_t reduce(std::unsigned_integral auto x){
static_assert(sizeof(x)<=8);
uint64_t a=x;
return a>=MOD?a-MOD:a;
}
static uint64_t sub_mod_raw(uint64_t a,uint64_t b){
#if defined(__x86_64__)
uint64_t res_wrapped=a;
uint64_t adjustment=b;
asm(
"sub %[y], %[x]\n\t"
"sbb %k[y], %k[y]\n\t"
:[x]"+r"(res_wrapped),
[y]"+r"(adjustment)
:
:"cc"
);
#else
uint64_t res_wrapped=a-b;
uint64_t adjustment=(res_wrapped>a)?EPS:0;
#endif
return res_wrapped-adjustment;
}
static uint64_t reduce_u160_raw(uint64_t lo,uint32_t mi,uint64_t hi){
return sub_mod_raw(sub_mod_raw(lo,hi),MOD-(uint64_t(mi)<<32)+mi);
}
static uint64_t reduce_u128_raw(__uint128_t v){
uint64_t hi=uint64_t(v>>64);
uint64_t lo=uint64_t(v);
uint32_t hi_hi=uint32_t(hi>>32);
uint32_t hi_lo=uint32_t(hi);
return reduce_u160_raw(lo,hi_lo,hi_hi);
}
Self&operator*=(Self o){
v=reduce_u128_raw(__uint128_t(v)*__uint128_t(o.v));
return*this;
}
};
template<typename T>T power(T a,long long b){
assert(b>=0);
T r=1;while(b){if(b&1)r*=a;b>>=1;a*=a;}return r;
}
template<typename U,typename V>struct pairnum:num_ops<pairnum<U,V>>{
using Self=pairnum;
U u;
V v;
pairnum():u(0),v(0){}
pairnum(long long val):u(val),v(val){}
pairnum(const U&u_,const V&v_):u(u_),v(v_){}
friend std::ostream&operator<<(std::ostream&out,const Self&n){return out<<'('<<n.u<<','<<' '<<n.v<<')';}
friend std::istream&operator>>(std::istream&in,Self&n){long long val;in>>val;n=Self(val);return in;}
friend bool operator==(const Self&a,const Self&b){return a.u==b.u&&a.v==b.v;}
Self inv()const{
return Self(u.inv(),v.inv());
}
Self neg()const{
return Self(u.neg(),v.neg());
}
Self&operator++(){
++u,++v;
return*this;
}
Self&operator--(){
--u,--v;
return*this;
}
Self&operator+=(const Self&o){
u+=o.u;
v+=o.v;
return*this;
}
Self&operator-=(const Self&o){
u-=o.u;
v-=o.v;
return*this;
}
Self&operator*=(const Self&o){
u*=o.u;
v*=o.v;
return*this;
}
Self&operator/=(const Self&o){
u/=o.u;
v/=o.v;
return*this;
}
};
template<typename tag>struct dynamic_modnum:mod_ops<dynamic_modnum<tag>,uint32_t>{
using Self=dynamic_modnum;
private:
inline static uint32_t MOD_=0;
inline static uint64_t BARRETT_M=0;
public:
static constexpr uint32_t const&MOD=MOD_;
using base=mod_ops<dynamic_modnum,uint32_t>;
using base::base;
using base::v;
using base::reduce;
static void set_mod(int mod){
assert(mod>0);
MOD_=uint32_t(mod);
BARRETT_M=(uint64_t(-1)/MOD);
}
static uint32_t barrett_reduce_partial(uint64_t a){
return uint32_t(a-uint64_t((__uint128_t(BARRETT_M)*a)>>64)*MOD);
}
static uint32_t barrett_reduce(uint64_t a){
int32_t res=int32_t(barrett_reduce_partial(a)-MOD);
return uint32_t((res<0)?res+int32_t(MOD):res);
}
struct mod_reader{
friend std::istream&operator>>(std::istream&i,mod_reader){
int mod;i>>mod;
Self::set_mod(mod);
return i;
}
};
static mod_reader MOD_READER(){
return mod_reader();
}
static uint32_t reduce(std::unsigned_integral auto x){
static_assert(sizeof(x)<=8);
return barrett_reduce(x);
}
explicit operator int()const{return int(v);}
Self&operator*=(const Self&o){
v=barrett_reduce(uint64_t(v)*o.v);
return*this;
}
};
template<typename T>struct mod_constraint{
T v,mod;
friend mod_constraint operator&(mod_constraint a,mod_constraint b){
if(a.mod<b.mod)std::swap(a,b);
if(b.mod==1)return a;
extended_gcd_result<T>egcd=extended_gcd<T>(a.mod,b.mod);
assert(a.v%egcd.gcd==b.v%egcd.gcd);
T extra=b.v-a.v%b.mod;
extra/=egcd.gcd;
extra*=egcd.coeff_a;
extra%=b.mod/egcd.gcd;
extra+=(extra<0)?b.mod/egcd.gcd:0;
return mod_constraint{
a.v+extra*a.mod,
a.mod*(b.mod/egcd.gcd)
};
}
};
// verify/dirichlet_inverse_and_prefix_sums.test.cpp
int main(){
std::ios_base::sync_with_stdio(false),std::cin.tie(nullptr);
int T;std::cin>>T;
while(T--){
int64_t N;std::cin>>N;
static dirichlet_series::div_vector_layout layout;
layout=N;
using num=modnum<998244353>;
using ds_prefix=dirichlet_series::prefix<layout,num>;
ds_prefix F;
for(int i=1;i<layout.len;i++)std::cin>>F.st[i];
ds_prefix delta([&](int64_t){return 1;});
ds_prefix H=delta/F;
for(int i=1;i<layout.len;i++)std::cout<<H.st[i]<<" \n"[i+1==layout.len];
}
return 0;
}
#pragma GCC diagnostic pop
// clang-format on
// @formatter:on
| Env | Name | Status | Elapsed | Memory |
|---|---|---|---|---|
| g++-sanitizer | example_00 |
|
13 ms | 8 MB |
| g++-sanitizer | example_01 |
|
12 ms | 8 MB |
| g++-sanitizer | many_00 |
|
451 ms | 37 MB |
| g++-sanitizer | many_01 |
|
445 ms | 37 MB |
| g++-sanitizer | many_02 |
|
464 ms | 37 MB |
| g++-sanitizer | max_00 |
|
10458 ms | 38 MB |
| g++-sanitizer | max_01 |
|
10566 ms | 38 MB |
| g++-sanitizer | max_02 |
|
10342 ms | 38 MB |
| g++-sanitizer | max_03 |
|
10306 ms | 38 MB |
| g++-sanitizer | max_04 |
|
10483 ms | 38 MB |
| g++-sanitizer | random_00 |
|
3835 ms | 24 MB |
| g++-sanitizer | random_01 |
|
6161 ms | 30 MB |
| g++-sanitizer | random_02 |
|
1493 ms | 18 MB |
| g++-sanitizer | random_03 |
|
9960 ms | 37 MB |
| g++-sanitizer | random_04 |
|
8269 ms | 34 MB |
| g++-sanitizer | small_00 |
|
18 ms | 10 MB |
| g++-sanitizer | small_01 |
|
16 ms | 9 MB |
| g++-sanitizer | small_02 |
|
15 ms | 9 MB |
| g++-sanitizer | small_03 |
|
17 ms | 10 MB |
| g++-sanitizer | small_04 |
|
15 ms | 10 MB |
| g++-sanitizer | very_small_00 |
|
14 ms | 8 MB |
| g++-sanitizer | very_small_01 |
|
14 ms | 8 MB |
| g++-sanitizer | very_small_02 |
|
13 ms | 8 MB |
| g++-sanitizer | very_small_03 |
|
11 ms | 8 MB |
| g++-sanitizer | very_small_04 |
|
10 ms | 8 MB |
| g++ | example_00 |
|
2 ms | 4 MB |
| g++ | example_01 |
|
2 ms | 4 MB |
| g++ | many_00 |
|
171 ms | 4 MB |
| g++ | many_01 |
|
171 ms | 4 MB |
| g++ | many_02 |
|
167 ms | 4 MB |
| g++ | max_00 |
|
3038 ms | 27 MB |
| g++ | max_01 |
|
3051 ms | 27 MB |
| g++ | max_02 |
|
3187 ms | 27 MB |
| g++ | max_03 |
|
3176 ms | 27 MB |
| g++ | max_04 |
|
3213 ms | 27 MB |
| g++ | random_00 |
|
1102 ms | 15 MB |
| g++ | random_01 |
|
1839 ms | 19 MB |
| g++ | random_02 |
|
424 ms | 9 MB |
| g++ | random_03 |
|
3092 ms | 26 MB |
| g++ | random_04 |
|
2446 ms | 23 MB |
| g++ | small_00 |
|
4 ms | 4 MB |
| g++ | small_01 |
|
3 ms | 4 MB |
| g++ | small_02 |
|
2 ms | 4 MB |
| g++ | small_03 |
|
3 ms | 4 MB |
| g++ | small_04 |
|
3 ms | 4 MB |
| g++ | very_small_00 |
|
2 ms | 4 MB |
| g++ | very_small_01 |
|
2 ms | 4 MB |
| g++ | very_small_02 |
|
2 ms | 4 MB |
| g++ | very_small_03 |
|
2 ms | 4 MB |
| g++ | very_small_04 |
|
2 ms | 3 MB |