This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub maspypy/library
#include "nt/gcd_convolution.hpp"
#include "nt/zeta.hpp" template <typename T> vc<T> gcd_convolution(vc<T> A, vc<T>& B) { assert(len(A) == len(B)); multiplier_zeta(A); multiplier_zeta(B); FOR(i, len(A)) A[i] *= B[i]; multiplier_mobius(A); return A; }
#line 2 "nt/primetable.hpp" template <typename T = int> vc<T> primetable(int LIM) { ++LIM; const int S = 32768; static int done = 2; static vc<T> primes = {2}, sieve(S + 1); if (done < LIM) { done = LIM; primes = {2}, sieve.assign(S + 1, 0); const int R = LIM / 2; primes.reserve(int(LIM / log(LIM) * 1.1)); vc<pair<int, int>> cp; for (int i = 3; i <= S; i += 2) { if (!sieve[i]) { cp.eb(i, i * i / 2); for (int j = i * i; j <= S; j += 2 * i) sieve[j] = 1; } } for (int L = 1; L <= R; L += S) { array<bool, S> block{}; for (auto& [p, idx]: cp) for (int i = idx; i < S + L; idx = (i += p)) block[i - L] = 1; FOR(i, min(S, R - L)) if (!block[i]) primes.eb((L + i) * 2 + 1); } } int k = LB(primes, LIM + 1); return {primes.begin(), primes.begin() + k}; } #line 3 "nt/zeta.hpp" template <typename T> void divisor_zeta(vc<T>& A) { assert(A[0] == 0); int N = len(A) - 1; auto P = primetable(N); for (auto&& p: P) { FOR3(x, 1, N / p + 1) A[p * x] += A[x]; } } template <typename T> void divisor_mobius(vc<T>& A) { assert(A[0] == 0); int N = len(A) - 1; auto P = primetable(N); for (auto&& p: P) { FOR3_R(x, 1, N / p + 1) A[p * x] -= A[x]; } } template <typename T> void multiplier_zeta(vc<T>& A) { assert(A[0] == 0); int N = len(A) - 1; auto P = primetable(N); for (auto&& p: P) { FOR3_R(x, 1, N / p + 1) A[x] += A[p * x]; } } template <typename T> void multiplier_mobius(vc<T>& A) { assert(A[0] == 0); int N = len(A) - 1; auto P = primetable(N); for (auto&& p: P) { FOR3(x, 1, N / p + 1) A[x] -= A[p * x]; } } #line 2 "nt/gcd_convolution.hpp" template <typename T> vc<T> gcd_convolution(vc<T> A, vc<T>& B) { assert(len(A) == len(B)); multiplier_zeta(A); multiplier_zeta(B); FOR(i, len(A)) A[i] *= B[i]; multiplier_mobius(A); return A; }