gcdConvlution (添字gcd畳み込み・約数系の高速ゼータ/メビウス変換)
(src/math/gcdConvolution.hpp)

• View this file on GitHub
• Last update: 2023-03-05 22:25:53+09:00
• Include: #include "src/math/gcdConvolution.hpp"

概要

列f、gに対して $\displaystyle h_i\ =\ \sum_{\gcd(j, k)=i} f_j\times g_k$ を満たす列 $h$ を求めます。

内部で約数集合の高速ゼータ変換、約数集合における高速メビウス変換をしているのでそちらを使いたい際もこれをincludeしてください。

約数集合の高速ゼータ変換とは列fに対して $\displaystyle F_i\ =\ \sum_{i\mid j} f_j$ を満たすFを求めることです。

約数集合の高速メビウス変換とは列Fに対して $f_i\ =\ \sum_{i\mid j} \mu (\frac{j}{i}) F_j$ を求めることです。 $\mu$ はメビウス関数を表します。この操作は約数集合の高速ゼータ変換の逆変換です。

機能

コンストラクタ

zawa::gcdConvolution<T>(std::size_t n)

内部で使用するエラトステネスの篩を初期化します。変換したい列のサイズの最大数を超える値を引数nに入れてください

計算量

$O(n\log (\log n))$

メンバ

fastZetaTransform

std::vector<T> fastZetaTransform(const std::vector<T>& f) const

約数集合の高速ゼータ変換を行います。変換後の列がstd::vector<T>の型で返されます。変換を行いたいstd::vectorの参照を引数にとります。

計算量

$O(n\log (\log n))$

fastMobiusTransform

std::vector<T> fastMobiusTransform(const std::vector<T>& F) const

約数集合の高速メビウス変換を行います。変換後の列がstd::vector<T>の型で返されます。変換を行いたいstd::vectorの参照を引数にとります。

計算量

$O(n\log (\log n))$

convolution

std::vector<T> convolution(const std::vector<T>& f, const std::vector<T>& f) const

添字gcd畳み込みを行います。結果をstd::vector<T>の型で返します。畳み込みを行いたい二つの列の参照を引数にとります。

計算量

$O(n\log (\log n))$

これは結局なんなの

TODO: NOTE

参考

# エラトステネスの篩の活用法を総特集！ 〜 高速素因数分解・メビウスの反転公式 〜

Depends on

• primeSieve (エラトステネスの篩） (src/math/primeSieve.hpp)

Verified with

• test/gcdConvolution.test.cpp

Code

#pragma once

#include "./primeSieve.hpp"

#include <vector>
#include <algorithm>

namespace zawa {

template <typename T>
class gcdConvolution {
private: 
	primeSieve sieve;

public:
	gcdConvolution(std::size_t n) : sieve(n) {}
	
	std::vector<T> fastZetaTransform(const std::vector<T>& f) const {
		std::vector<T> res(f.begin(), f.end());
		for (int i = 1 ; i <= (int)(f.size()) ; i++) {
			if (sieve[i]) {
				for (int j = (int)res.size() / i ; j >= 1 ; j--) {
					res[j - 1] += res[j * i - 1];
				}
			}
		}
		return res;
	}

	std::vector<T> fastMobiusTransform(const std::vector<T>& F) const {
		std::vector<T> res(F.begin(), F.end());
		for (int i = 1 ; i <= (int)(F.size()) ; i++) {
			if (sieve[i]) {
				for (int j = 1 ; j * i <= (int)(F.size()) ; j++) {
					res[j - 1] -= res[j * i - 1];
				}
			}
		}
		return res;
	}

	std::vector<T> convolution(const std::vector<T>& f, const std::vector<T>& g) const {
		std::vector<T> F = fastZetaTransform(f);
		std::vector<T> G = fastZetaTransform(g);
		std::vector<T> H(std::min(F.size(), G.size()));
		for (int i = 0 ; i < (int)(H.size()) ; i++) {
			H[i] = F[i] * G[i];
		}
		return fastMobiusTransform(H);
	}
};

}// namespace zawa

#line 2 "src/math/gcdConvolution.hpp"

#line 2 "src/math/primeSieve.hpp"

#include <vector>

namespace zawa {

class primeSieve {
private:
	std::vector<bool> sieve;

public:
	primeSieve() {}
	primeSieve(std::size_t n) : sieve(n + 1, true) {
		if (n >= 0) {
			sieve[0] = false;
		}
		if (n >= 1) {
			sieve[1] = false;
		}
		for (std::size_t i = 2 ; i <= n ; i++) {
			if (sieve[i]) {
				for (std::size_t j = i * 2 ; j <= n ; j += i) {
					sieve[j] = false;
				}
			}
		}
	}

	inline bool operator[](std::size_t i) const {
		return sieve[i];
	}

	inline std::size_t size() const {
		return sieve.size();
	}
};

}// namespace zawa
#line 4 "src/math/gcdConvolution.hpp"

#line 6 "src/math/gcdConvolution.hpp"
#include <algorithm>

namespace zawa {

template <typename T>
class gcdConvolution {
private: 
	primeSieve sieve;

public:
	gcdConvolution(std::size_t n) : sieve(n) {}
	
	std::vector<T> fastZetaTransform(const std::vector<T>& f) const {
		std::vector<T> res(f.begin(), f.end());
		for (int i = 1 ; i <= (int)(f.size()) ; i++) {
			if (sieve[i]) {
				for (int j = (int)res.size() / i ; j >= 1 ; j--) {
					res[j - 1] += res[j * i - 1];
				}
			}
		}
		return res;
	}

	std::vector<T> fastMobiusTransform(const std::vector<T>& F) const {
		std::vector<T> res(F.begin(), F.end());
		for (int i = 1 ; i <= (int)(F.size()) ; i++) {
			if (sieve[i]) {
				for (int j = 1 ; j * i <= (int)(F.size()) ; j++) {
					res[j - 1] -= res[j * i - 1];
				}
			}
		}
		return res;
	}

	std::vector<T> convolution(const std::vector<T>& f, const std::vector<T>& g) const {
		std::vector<T> F = fastZetaTransform(f);
		std::vector<T> G = fastZetaTransform(g);
		std::vector<T> H(std::min(F.size(), G.size()));
		for (int i = 0 ; i < (int)(H.size()) ; i++) {
			H[i] = F[i] * G[i];
		}
		return fastMobiusTransform(H);
	}
};

}// namespace zawa

Back to top page