tkmst201's Library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub tkmst201/Library

:heavy_check_mark: Test/zeta_moebius_transform.set.1.test.cpp

Depends on

Code

#define PROBLEM "https://onlinejudge.u-aizu.ac.jp/challenges/sources/JAG/Summer/2446?year=2012"
#define ERROR "1e-7"

#include "Mathematics/zeta_moebius_transform.hpp"
#include "Mathematics/euclid.hpp"

#include <vector>
#include <iostream>
#include <cmath>

int main() {
	using ll = long long;
	
	int N;
	ll M;
	std::cin >> N >> M;
	std::vector<ll> A(N);
	std::vector<int> P(N);
	for (int i = 0; i < N; ++i) std::cin >> A[i];
	for (int i = 0; i < N; ++i) std::cin >> P[i];
	
	auto solve1 = [&]() {
		std::vector<ll> v((1 << N));
		
		for (int s = 1; s < 1 << N; ++s) {
			ll l = 1;
			bool ng = false;
			for (int i = 0; i < N; ++i) {
				if (~s >> i & 1) continue;
				ll g = tk::gcd<ll>(l, A[i]);
				l /= g;
				if (l > M / A[i]) { ng = true; break; } // l A[i] <= M

				else l *= A[i];
			}
			if (ng) continue;
			v[s] = M / l;
		}
		
		tk::moebius_transform_set_lower(v);
		
		double ans = 0;
		
		for (int s = 0; s < 1 << N; ++s) {
			double p = 1;
			for (int i = 0; i < N; ++i) {
				if (s >> i & 1) p *= P[i] * 0.01;
				else p *= (100 - P[i]) * 0.01;
			}
			ans += std::abs(v[s]) * p;
		}
		return ans;
	};
	
	auto solve2 = [&]() {
		std::vector<ll> v((1 << N));
		for (int s = 0; s < 1 << N; ++s) {
			ll l = 1;
			bool ng = false;
			for (int i = 0; i < N; ++i) {
				if (~s >> i & 1) continue;
				ll g = tk::gcd<ll>(l, A[i]);
				l /= g;
				if (l > M / A[i]) { ng = true; break; } // l A[i] <= M

				else l *= A[i];
			}
			if (ng) continue;
			v[s] = M / l;
		}
		
		tk::moebius_transform_set_upper(v);
		tk::zeta_tranform_set_lower(v);
		
		double ans = 0;
		
		for (int s = 0; s < 1 << N; ++s) {
			double p = 1;
			for (int i = 0; i < N; ++i) {
				if (s >> i & 1) p *= P[i] * 0.01;
				else p *= (100 - P[i]) * 0.01;
			}
			ans += (M - v[(1 << N) - 1 - s]) * p;
		}
		return ans;
	};
	
	double ans1 = solve1();
	double ans2 = solve2();
	printf("%.16f\n", (ans1 + ans2) / 2.0);
}
#line 1 "Test/zeta_moebius_transform.set.1.test.cpp"
#define PROBLEM "https://onlinejudge.u-aizu.ac.jp/challenges/sources/JAG/Summer/2446?year=2012"
#define ERROR "1e-7"

#line 1 "Mathematics/zeta_moebius_transform.hpp"



/*
last-updated: 2020/11/19

TODO: AOJ, yukicoder から verify 問題を探す

# 仕様
## 集合
void zeta_tranform_set_lower(std::vector<T> & v)
	時間計算量: Θ(|v|log|v|)
	配列 v を配列 v' に書き換える(i の下位集合について総和を取る)
	v'[i] := \Sum_{j \subset i} v[j]

void moebius_transform_set_lower(std::vector<T> & v)
	時間計算量: Θ(|v|log|v|)
	配列 v を次の式を満たす配列 v' に書き換える(i の下位集合について総和を取る前に戻す)
	v[i] = \Sum_{j \subset i} v'[j]

void zeta_tranform_set_upper(std::vector<T> & v)
	時間計算量: Θ(|v|log|v|)
	配列 v を配列 v' に書き換える(i の上位集合について総和を取る)
	v'[i] := \Sum_{i \subset j} v[j]

void moebius_transform_set_upper(std::vector<T> & v)
	時間計算量: Θ(|v|log|v|)
	配列 v を次の式を満たす配列 v' に書き換える(i の上位集合について総和を取る前に戻す)
	v[i] = \Sum_{i \subset j} v'[j]

## 約数倍数
void zeta_tranform_divisor_lower(std::vector<T> & v)
	時間計算量: O(|v|loglog|v|)
	配列 v を配列 v' に書き換える(i の約数について総和を取る)
	v'[i] := \Sum_{j|i} v[j]

void moebius_tranform_divisor_lower(std::vector<T> & v)
	時間計算量: O(|v|loglog|v|)
	配列 v を次の式を満たす配列 v' に書き換える(i の約数について総和を取る前に戻す)
	v[i] = \Sum_{j|i} v'[j]

void zeta_tranform_divisor_upper(std::vector<T> & v)
	時間計算量: O(|v|loglog|v|)
	配列 v を配列 v' に書き換える(i の倍数について総和)
	v'[i] := \Sum_{i|j} v[j]

void moebius_tranform_divisor_upper(std::vector<T> & v)
	時間計算量: O(|v|loglog|v|)
	配列 v を次の式を満たす配列 v' に書き換える(i の約数について総和を取る前に戻す)
	v[i] = \Sum_{i|j} v'[j]

# 参考
2020/11/12: https://qiita.com/convexineq/items/afc84dfb9ee4ec4a67d5
2020/11/14: https://aprilganmo.hatenablog.com/entry/2020/02/27/170239
2020/11/14: https://noshi91.hatenablog.com/entry/2018/12/27/121649
2020/11/18: https://aprilganmo.hatenablog.com/entry/2020/07/24/190816

# verify
void zeta_tranform_set_lower(std::vector<T> & v)
void moebius_transform_set_lower(std::vector<T> & v)
void zeta_tranform_set_upper(std::vector<T> & v)
- OK

void moebius_transform_set_upper(std::vector<T> & v)
- TODO: メビウス(上位集合) の verify

void zeta_tranform_divisor_lower(std::vector<T> & v)
- https://atcoder.jp/contests/abc172/submissions/18200046

void moebius_tranform_divisor_lower(std::vector<T> & v)
- TODO: メビウス(倍数) の verify

void zeta_tranform_divisor_upper(std::vector<T> & v)
void moebius_tranform_divisor_upper(std::vector<T> & v)
- https://atcoder.jp/contests/abc162/submissions/18200657
*/

#include <vector>
#include <cassert>

namespace tk {
template<typename T>
void zeta_tranform_set_lower(std::vector<T> & v) {
	assert(!v.empty());
	using size_type = std::size_t;
	for (size_type i = 1; i < v.size(); i <<= 1) {
		for (size_type j = 0; j < v.size(); ++j) {
			if (j & i) v[j] += v[j ^ i];
		}
	}
}

template<typename T>
void moebius_transform_set_lower(std::vector<T> & v) {
	assert(!v.empty());
	using size_type = std::size_t;
	for (size_type i = 1; i < v.size(); i <<= 1) {
		for (size_type j = 0; j < v.size(); ++j) {
			if (j & i) v[j] -= v[j ^ i];
		}
	}
}

template<typename T>
void zeta_tranform_set_upper(std::vector<T> & v) {
	assert(!v.empty());
	using size_type = std::size_t;
	for (size_type i = 1; i < v.size(); i <<= 1) {
		for (size_type j = 0; j < v.size(); ++j) {
			if (~j & i) v[j] += v[j ^ i];
		}
	}
}

template<typename T>
void moebius_transform_set_upper(std::vector<T> & v) {
	assert(!v.empty());
	using size_type = std::size_t;
	for (size_type i = 1; i < v.size(); i <<= 1) {
		for (size_type j = 0; j < v.size(); ++j) {
			if (~j & i) v[j] -= v[j ^ i];
		}
	}
}
} // namespace tk


namespace tk {
template<typename T>
void zeta_tranform_divisor_lower(std::vector<T> & v) {
	assert(!v.empty());
	using size_type = std::size_t;
	std::vector<bool> sieve(v.size(), true);
	for (size_type p = 2; p < v.size(); ++p) {
		if (!sieve[p]) continue;
		for (size_type i = 1, t = p; t < v.size(); ++i, t += p) {
			v[t] += v[i];
			sieve[t] = false;
		}
	}
}

template<typename T>
void moebius_tranform_divisor_lower(std::vector<T> & v) {
	assert(!v.empty());
	using size_type = std::size_t;
	std::vector<bool> sieve(v.size(), true);
	for (size_type p = 2; p < v.size(); ++p) {
		if (!sieve[p]) continue;
		for (size_type i = (v.size() - 1) / p, t = i * p; i > 0; --i, t -= p) {
			v[t] -= v[i];
			sieve[t] = false;
		}
	}
}

template<typename T>
void zeta_tranform_divisor_upper(std::vector<T> & v) {
	assert(!v.empty());
	using size_type = std::size_t;
	std::vector<bool> sieve(v.size(), true);
	for (size_type p = 2; p < v.size(); ++p) {
		if (!sieve[p]) continue;
		for (size_type i = (v.size() - 1) / p, t = i * p; i > 0; --i, t -= p) {
			v[i] += v[t];
			sieve[t] = false;
		}
	}
}

template<typename T>
void moebius_tranform_divisor_upper(std::vector<T> & v) {
	assert(!v.empty());
	using size_type = std::size_t;
	std::vector<bool> sieve(v.size(), true);
	for (size_type p = 2; p < v.size(); ++p) {
		if (!sieve[p]) continue;
		for (size_type i = 1, t = p; t < v.size(); ++i, t += p) {
			v[i] -= v[t];
			sieve[t] = false;
		}
	}
}
} // namespace tk



#line 1 "Mathematics/euclid.hpp"



#line 5 "Mathematics/euclid.hpp"
#include <utility>
#include <tuple>
#include <type_traits>
#include <cmath>

/**
 * @brief https://tkmst201.github.io/Library/Mathematics/euclid.hpp
 */
namespace tk {
template<typename T>
constexpr T gcd(T a, T b) noexcept {
	static_assert(std::is_integral<T>::value);
	assert(a >= 0);
	assert(b >= 0);
	while (b != 0) {
		const T t = a % b;
		a = b; b = t;
	}
	return a;
}

template<typename T>
constexpr T lcm(T a, T b) noexcept {
	static_assert(std::is_integral<T>::value);
	assert(a >= 0);
	assert(b >= 0);
	if (a == 0 || b == 0) return 0;
	return a / gcd(a, b) * b;
}

template<typename T>
constexpr std::tuple<T, T, T> ext_gcd(T a, T b) noexcept {
	static_assert(std::is_integral<T>::value);
	static_assert(std::is_signed<T>::value);
	assert(a != 0);
	assert(b != 0);
	T a1 = (a > 0) * 2 - 1, a2 = 0, b1 = 0, b2 = (b > 0) * 2 - 1;
	a = std::abs(a);
	b = std::abs(b);
	while (b > 0) {
		const T q = a / b;
		T tmp = a - q * b; a = b; b = tmp;
		tmp = a1 - q * b1; a1 = b1; b1 = tmp;
		tmp = a2 - q * b2; a2 = b2; b2 = tmp;
	}
	return {a, a1, a2};
}
} // namespace tk



#line 6 "Test/zeta_moebius_transform.set.1.test.cpp"

#line 8 "Test/zeta_moebius_transform.set.1.test.cpp"
#include <iostream>
#line 10 "Test/zeta_moebius_transform.set.1.test.cpp"

int main() {
	using ll = long long;
	
	int N;
	ll M;
	std::cin >> N >> M;
	std::vector<ll> A(N);
	std::vector<int> P(N);
	for (int i = 0; i < N; ++i) std::cin >> A[i];
	for (int i = 0; i < N; ++i) std::cin >> P[i];
	
	auto solve1 = [&]() {
		std::vector<ll> v((1 << N));
		
		for (int s = 1; s < 1 << N; ++s) {
			ll l = 1;
			bool ng = false;
			for (int i = 0; i < N; ++i) {
				if (~s >> i & 1) continue;
				ll g = tk::gcd<ll>(l, A[i]);
				l /= g;
				if (l > M / A[i]) { ng = true; break; } // l A[i] <= M

				else l *= A[i];
			}
			if (ng) continue;
			v[s] = M / l;
		}
		
		tk::moebius_transform_set_lower(v);
		
		double ans = 0;
		
		for (int s = 0; s < 1 << N; ++s) {
			double p = 1;
			for (int i = 0; i < N; ++i) {
				if (s >> i & 1) p *= P[i] * 0.01;
				else p *= (100 - P[i]) * 0.01;
			}
			ans += std::abs(v[s]) * p;
		}
		return ans;
	};
	
	auto solve2 = [&]() {
		std::vector<ll> v((1 << N));
		for (int s = 0; s < 1 << N; ++s) {
			ll l = 1;
			bool ng = false;
			for (int i = 0; i < N; ++i) {
				if (~s >> i & 1) continue;
				ll g = tk::gcd<ll>(l, A[i]);
				l /= g;
				if (l > M / A[i]) { ng = true; break; } // l A[i] <= M

				else l *= A[i];
			}
			if (ng) continue;
			v[s] = M / l;
		}
		
		tk::moebius_transform_set_upper(v);
		tk::zeta_tranform_set_lower(v);
		
		double ans = 0;
		
		for (int s = 0; s < 1 << N; ++s) {
			double p = 1;
			for (int i = 0; i < N; ++i) {
				if (s >> i & 1) p *= P[i] * 0.01;
				else p *= (100 - P[i]) * 0.01;
			}
			ans += (M - v[(1 << N) - 1 - s]) * p;
		}
		return ans;
	};
	
	double ans1 = solve1();
	double ans2 = solve2();
	printf("%.16f\n", (ans1 + ans2) / 2.0);
}
Back to top page