#include "dirichlet_series.hpp"

View this file on GitHub · Last update: 2026-07-26 16:57:25-07:00

Required by

• 🧪 dirichlet_series.test.cpp

Verified with

• verify/dirichlet_convolution_and_prefix_sums.test.cpp
• verify/dirichlet_inverse_and_prefix_sums.test.cpp
• verify/sum_of_totient_function.test.cpp

Code

	Coverage	Exec / Excl / Total
Lines	77.6%	318 / 0 / 410
Functions	98.5%	67 / 0 / 68
Branches	77.8%	535 / 0 / 688
Full report

#pragma once

#include <cstdint>
#include <algorithm>
#include <cassert>
#include <type_traits>
#include <array>

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);
	}
};

}

#include <cstdint>
#include <algorithm>
#include <cassert>
#include <type_traits>
#include <array>
#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);
	}
};

}

// 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/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);
}
};
}
#pragma GCC diagnostic pop
// clang-format on
// @formatter:on

#pragma once

#include <cstdint>
#include <algorithm>
#include <cassert>
#include <type_traits>
#include <array>

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);
	}
};

}

Back to top page