最古の遺跡解説

この問題です

ジャッジ:C - 最古の遺跡
問題:https://www.ioi-jp.org/joi/2006/2007-ho-prob_and_sol/2007-ho.pdf

オーダー的にO(N3)が辛くてO(N2)が通るので、O(N2)を考えます。 二点を固定して、残りの二頂点C,Dに点があるかを判定すれば良いです。
f:id:eiya5498513:20171008203034j:plain

まずは、C1,D1に頂点があるかを確認します。
この後はC1,D1をただC,Dと呼びます。

ここで予めある座標に頂点があるかどうかをbool配列でもっておくと、あとはC1,D1の座標を求めるだけであるなしが判定できます。

なんとなく正方形ABCDを回転していない正方形で囲ってみます。(なぜこの発想になるのかの説明は余裕があったら追記します...)
f:id:eiya5498513:20171008203947j:plain

x = Bx - Ax、y = By - Ayとなるようにx,yを取ります。そうするとx,yは上図のようになります。
ここで、三角形ABO1≡三角形CAO2なので、CO2==x、O2A==yです。
f:id:eiya5498513:20171008204448j:plain

よって、Cx = Ax - y、Cy = Ay + xで求められることが分かります。
又、DはDx = Cx + x、Dy = Cy + yで求まります。

この図はA->Bが右上がりの時の図なので、A->Bが右下がりの時どうなるか考えてみます。

どんどんBを下げていくと、こうなります。
f:id:eiya5498513:20171008205141j:plain
この時は三角形が作れないですが、無理やりy==0と考えてやると、上手く行くことが分かります。

Bをもっと下げてみます。
f:id:eiya5498513:20171008205336j:plain
図の形が変わりました。
xとyをそのまま負数に拡張してやると、図のようになり、式を変えずにそのまま計算することが出来ます。

xが負の時も同じように考えると、式を全く変形せずにCx = Ax - y、Cy = Ay + x、Dx = Cx + x、Dy = Cy + yの式で良いことが分かります。

次にC2,D2について考えます。
Aを中心にxとyを反転させる(Aを中心に紙を180度回す)と図のようになります。 回転後のBをA'、回転後のAをB'と考えると、C1',D1'がそれぞれC2,D2に対応します。
よって、「AとBを入れ替えてもう一度C1,D1を求める処理をする」とC2,D2の座標を求めることが出来ます。

ここで、C2,D2を求めるときに行っているのはAとBを交換する作業であることから、AとBの選び方をいつものように組み合わせを考えるときの以下のfor分

for(int a=0;a < N;++a)for(int b = a+1;b < N;++b)

を使うと思いますが、今回は

for(int a=0;a < N;++a)for(int b = 0;b < N;++b){
if(a==b){continue;}
//処理
}

とa,bの順番も考慮してやることで、C1,D1の座標を求めるコードのみにしてやることが出来ます。

ところで、a==bの時は何が起こるのでしょう。
この時は同じ頂点二つを使って正方形を作ろうとし、その4頂点は全て同じ座標に存在するので、正方形が作れてしまいます。
しかし実際は4頂点を選んで正方形を作らないといけないので、作ってはいけない正方形を作ってしまいます。これは弾く必要が有ります。
しかしよく考えると、この時正方形の面積は0、つまり正方形の最小の面積1よりも小さいので、実は「正方形が存在する場合は」弾いてやらなくても面積のmaxを取る段階で弾くことが出来ます。
逆に面積のmaxが0のときは正方形が一つも存在しないということなので、例外処理である0を出力する必要が有ります。
が、値が0の時に0を出力するという処理になっているので、結局そのまま出力すれば良いことが分かります。

この方針でコードを実装してやるとこうなります。短く書けるのでバグを埋め込みにくくなります。
Submission #1671065 - 第6回日本情報オリンピック 本選(オンライン)

tupleを展開して関数に渡す関数

tupleを展開して関数に渡す関数
boostを使えば出来るらしいけど、知らない。
出来栄えはかなり雑です。

template <class F, class... Ts, class... Us>
typename std::enable_if<
	sizeof...(Us) == sizeof...(Ts),
	std::result_of_t<F(Ts...)>>::type
	apply_impl(F&& fun, std::tuple<Ts...>&, Us&... us)
{
	return fun(us...);
}
template <class F, class... Ts, class... Us>
typename std::enable_if<
	sizeof...(Us) < sizeof...(Ts),
	std::result_of_t<F(Ts...)>>::type
	apply_impl(F&& fun, std::tuple<Ts...>& args, Us&... us)
{
	return apply_impl(fun, args, us..., std::get<sizeof...(Us)>(args));
}
template <class F, class... Ts>
std::result_of_t<F(Ts...)> apply(F&& fun, std::tuple<Ts...>& args)
{
	return apply_impl(fun, args);
}

使用例

template <class... Ts>
auto pop_front(std::tuple<Ts...>& args)
{
	return apply([](auto&&, auto&&... v) {return std::make_tuple(v...); }, args);
}

プリム法の証明

かなり乱暴なことをしている気がします。厳密な証明は他をあたってください。

プリム法とは

最小全域木を構成するアルゴリズムの一つです。
以下の手順で行います。

  1. ある頂点を選び、木Tに追加します
  2. 木Tから最も近い頂点を木Tに追加します
  3. Tが全域木になるまで2を繰り返します
  4. Tは最小全域木です

プリム法の証明

1が終わった段階で、木Tは頂点を一つ含む木です。このとき明らかに、木Tは最小全域木の部分木です。

木Tが最小全域木の部分木であるとします。
木Tから、Tに含まれない頂点に対する辺の内、最もコストが小さい辺をe、その辺の結ぶ頂点の内Tに含まれない方をvとします。
f:id:eiya5498513:20170830211615p:plain
操作2では木Tに辺eを追加します。これは、最小全域木に辺eが含まれるからです。
これを背理法で証明します。
最小全域木に辺eが含まれないと仮定します。
このとき、辺eを含まないある全域木Uが最小全域木だとします。
f:id:eiya5498513:20170830210508p:plain
Uにはvが含まれているので、Uには「Tに含まれる頂点->(Tに含まれない頂点Xnが0個以上)->v」という辺eを使わない経路が存在します。(図の緑+黄色)
(しない場合は、vへの辺が無いことになるので、Uが全域木であることに矛盾します。)
f:id:eiya5498513:20170830211352p:plain
この「Tに含まれる頂点->X1(又はv)」の辺をfとします。
最小全域木からfを取り除いて、eを追加した木は、明らかに全域木です。
又、この全域木は、Uよりもコストの合計が低いです。
これは、Uが最小全域木であるので矛盾します。
よって、「最小全域木に辺eが含まれない」という仮定は誤りです。
よって、最小全域木に辺eが含まれます。
よって、木Tに辺eを追加した木も最小全域木の部分木です。

おまけ

O(V^2)プリム法のコード(ほぼ蟻本のコピペ)

int cost[MAX_V][MAX_V]; // cost[u][v]は辺e=(u,v)のコスト(存在しない場合はINF)
int mincost[MAX_V]; // 集合Xからの辺の最小コスト
bool used[MAX_V]; // 頂点iがXに含まれているか
int V; // 頂点数
int prim() {
	for (int i = 0; i < V; ++i) {
		mincost[i] = INF;
		used[i] = false;
	}
	mincost[0] = 0;
	int res = 0;
	while (true) {
		int v = -1;
		// Xに属さない頂点のうちXからの辺のコストが最小になる頂点を探す
		for (int u = 0; u < V; u++) {
			if (!used[u] && (v == -1 || mincost[u] < mincost[v])) v = u;
		}
		if (v == -1) break;
		used[v] = true; // 頂点vをXに追加
		res += mincost[v]; // 辺のコストを加える
		for (int u = 0; u < V; u++) {
			mincost[u] = std::min(mincost[u], cost[v][u]);
		}
	}
	return res;
}

O(ElogE)プリム法(自作)
O(ElogV)が存在するはずなので、誰か教えてください...
恐らく、二重辺を認めないor消去して、O(ElogE)→O(Elog(V^2))→O(ElogV)とやると思われる。(Luzhiledプロありがとうございます)

#include <queue>
#include <vector>
#include <functional>
#include <utility>
#include <algorithm>
#include <iterator>

using COST_T = uint32_t;
constexpr uint32_t MAX_V = 変える;
constexpr COST_T INF = 変える;//std::numeric_limits<double>::infinity()

#if defined(_MSC_VER) && defined(_DEBUG)
static_assert(false, "リリースでコンパイルしないと遅いよ!!");
#endif

struct edge {
	uint32_t to;
	COST_T cost;
	edge() {}
	edge(uint32_t to_, COST_T cost_)
		:to(to_), cost(cost_) {}
};
std::vector<edge> graph[MAX_V];// 頂点vからの辺はgraph[v]
//O(V+ElogE)
COST_T prim(uint32_t s = 0) {
	static bool used[MAX_V]; // 頂点iがTに含まれているか
	for (int i = 0; i < MAX_V; ++i) {
		used[i] = false;
	}
	using Pair = std::pair<COST_T, uint32_t>;
	std::priority_queue<Pair, std::vector<Pair>, std::greater<Pair>> queue;
	queue.emplace(0, s);
	COST_T res = 0;
	while (!queue.empty()) {
		// Tからの辺のコストが最小になる頂点
		COST_T cost = queue.top().first;
		uint32_t v = queue.top().second; queue.pop();
		if (used[v]) { continue; }

		used[v] = true; // 頂点vをTに追加
		res += cost; // 辺のコストを加える
		for (auto& next : graph[v]) {
			if (used[next.to]) { continue; }
			queue.emplace(next.cost, next.to);
		}
	}
	return res;
}

O(ElogV)ダイクストラのテンプレ

O(ElogV)ダイクストラのテンプレ

#include <queue>
#include <vector>
#include <functional>
#include <utility>
#include <algorithm>
#include <iterator>

using COST_T = uint32_t;
constexpr uint32_t N_MAX = 変える;
constexpr COST_T INF = 変える;//std::numeric_limits<double>::infinity()

#if defined(_MSC_VER) && defined(_DEBUG)
static_assert(false, "リリースでコンパイルしないと遅いよ!!");
#endif

struct edge {
	uint32_t to;
	COST_T cost;
	edge() {}
	edge(uint32_t to_, COST_T cost_)
		:to(to_), cost(cost_) {}
};
std::vector<edge> graph[N_MAX];

//O(ElogV)ダイクストラ
COST_T D[N_MAX];
void Dijkstra(uint32_t s)
{
	using P = std::pair<COST_T, uint32_t>;//cost pos
	std::priority_queue<P, std::vector<P>, std::greater<P> > que;
	std::fill(std::begin(D), std::end(D), INF);

	D[s] = 0;
	que.emplace(0, s);
	while (!que.empty())
	{
		auto p = que.top(); que.pop();
		const auto& nowpos = p.second;
		const auto& nowcost = p.first;
		if (D[nowpos] < nowcost) { continue; }
		for (const auto& e : graph[nowpos])
		{
			auto cost = nowcost + e.cost;
			if (cost < D[e.to]) {
				D[e.to] = cost;
				que.emplace(cost, e.to);
			}
		}
	}
}

使用例

Dijkstra(0);

N次元配列を同じ値で埋めるテンプレ2

N次元配列を同じ値で埋めるテンプレ
キャストが要らないバージョン
シンプルなのはこっち:N次元配列を同じ値で埋めるテンプレ - 永夜の記録

#include <type_traits>
template<typename T, typename U>
typename std::enable_if<std::rank<T>::value == 0>::type fill_all(T& arr, const U& v) {
	arr = v;
}
template<typename ARR, typename U>
typename std::enable_if<std::rank<ARR>::value!=0>::type fill_all(ARR& arr, const U& v) {
	for (auto& i : arr) {
		fill_all(i, v);
	}
}

使用例

long long dp[1000][1000];
fill_all(dp,-1);//キャストが要らない

C++14用mod_int

C++14用mod_intです。勝手にmodを取ります。
型名はmintです。

using mint = mint_base<1000000007>;

の1000000007がとるmodになります。
+-*/やインクリメント,デクリメントが出来ます。
~で逆元を取れますが、そのままでも割り算が使えます。(どちらもlogMODで、偶に定数倍でTLEするので、その時は事前計算してください)
ストリーム演算子があるので、>>や<<で入出力が出来ます。
繰り返し二乗法(ipow),階乗(fact_set),組み合わせ(combination)も一緒に入ってます。
ユーザー定義リテラル_miがあるので、mint(1)の代わりに1_miと出来ます。

#include <algorithm>
#include <utility>
#include <type_traits>
#include <cstdint>
template<typename Arithmetic, typename Integral>
std::enable_if_t< std::is_unsigned<Integral>::value, Arithmetic>
ipow(Arithmetic bace, Integral n)
{
	//繰り返し二条法
	auto res = (Arithmetic)(1);
	while (n > 0) {
		if (n & 1) res *= bace;
		bace *= bace;
		n >>= 1;
	}
	return res;
}
constexpr bool is_prime(uint32_t N)
{
	if (N <= 1) {
		return false;
	}
	for (size_t i = 2; i*i <= N; ++i)
	{
		if (N%i == 0) {
			return false;
		}
	}
	return true;
}
template <uint64_t MOD> class mint_base;
//mint_base_base型用の累乗関数
template <uint64_t MOD> constexpr mint_base<MOD> m_pow(mint_base<MOD> x, uint64_t n)noexcept;
//mod計算を自動で行う整数テンプレートクラス
template <uint64_t MOD_ = 1000000007>
class mint_base
{
public:
	static constexpr auto MOD = MOD_;
	static_assert(!(MOD <= 2), "MOD cannot be below 2.");
	static_assert(MOD <= (0xFFFFFFFFFFFFFFFF / 2), "MOD is too big");//加算してオーバーフローしない
	static_assert(MOD <= 0xFFFFFFFF, "MOD is too big");//乗算してオーバーフローしない
	constexpr mint_base<MOD> operator+(const mint_base<MOD> &other)const noexcept
	{
		auto v = *this;
		return v += other;
	}
	constexpr mint_base<MOD> operator-(const mint_base<MOD> &other)const noexcept
	{
		auto v = *this;
		return v -= other;
	}
	constexpr mint_base<MOD> operator*(const mint_base<MOD> &other)const noexcept
	{
		auto v = *this;
		return v *= other;
	}
	constexpr auto operator/(const mint_base<MOD> &other)const noexcept
	{
		auto v = *this;
		return v /= other;
	}
	constexpr mint_base<MOD>& operator+=(const mint_base<MOD> &other) noexcept
	{
		a += other.a;
		if (MOD <= a) { a -= MOD; };
		return *this;
	}
	constexpr mint_base<MOD>& operator-=(const mint_base<MOD> &other) noexcept
	{
		if (a >= other.a) {
			a -= other.a;
		}
		else {
			a = (a + MOD) - other.a;
		}
		return *this;
	}
	constexpr mint_base<MOD>& operator*=(const mint_base<MOD> &other) noexcept
	{
#if 1
		a *= other.a;
		a %= MOD;
#else
		//MOD <= (MAXUINT64 / 2)条件下
		uint64_t b = other.a, v = 0;
		while (b > 0) {
			if (b & 1) {
				v += a;
				if (v >= MOD)v -= MOD;
			}
			a += a;
			if (MOD <= a)a -= MOD;
			b >>= 1;
		}
		a = v;
#endif
		return *this;
	}
	constexpr mint_base<MOD>& operator/=(const mint_base<MOD> &other) noexcept
	{
		return *this *= ~other;
	}
	constexpr mint_base<MOD> operator+()const noexcept { return *this; }
	constexpr mint_base<MOD> operator-()const noexcept
	{
		return{ MOD - a, mod_value_tag{} };
	}
	constexpr mint_base<MOD>& operator++() noexcept
	{
		if (MOD <= ++a) { a = 0; };
		return *this;
	}
	constexpr mint_base<MOD>& operator--() noexcept
	{
		if (a <= 0) { a = MOD; };
		--a;
		return *this;
	}
	constexpr mint_base<MOD> operator++(int) noexcept
	{
		auto tmp = *this;
		++*this;
		return tmp;
	}
	constexpr mint_base<MOD> operator--(int) noexcept
	{
		auto tmp = *this;
		--*this;
		return tmp;
	}
	constexpr mint_base<MOD> operator~()const noexcept
	{
		return ipow(*this, e_phi - 1);
	}
	constexpr mint_base<MOD>& operator=(const mint_base<MOD> &other) noexcept
	{
		a = other.a;
		return *this;
	}
	constexpr explicit operator uint64_t()const noexcept
	{
		return a;
	}
	constexpr explicit operator unsigned()const noexcept
	{
		return (unsigned)a;
	}
	static constexpr uint64_t getmod() noexcept
	{
		return MOD;
	}
	constexpr mint_base(uint64_t a_) noexcept :a(a_ % MOD) {}
	constexpr mint_base()noexcept : a(0) {}
	struct mod_value_tag {};
	constexpr mint_base(uint64_t a_, mod_value_tag) :a(a_) {}
private:
	static constexpr uint64_t get_e_phi()noexcept {
		//オイラー値の導出
		uint64_t temp = MOD;
		uint64_t m_ = MOD;
		for (uint64_t i = 2; i * i <= m_; ++i)
		{
			if (m_ % i == 0)
			{
				temp = temp / i * (i - 1);
				for (; m_ % i == 0; m_ /= i);
			}
		}
		if (m_ != 1)temp = temp / m_ * (m_ - 1);
		return temp;
	}
	static constexpr uint64_t e_phi = get_e_phi();//オイラー値
	uint64_t a;
};
//mint_base型用の累乗関数
template<uint64_t MOD>constexpr mint_base<MOD> m_pow(mint_base<MOD> x, uint64_t n)noexcept
{
	mint_base<MOD> res = 1;
	while (n > 0)
	{
		if (n & 1)res *= x;
		x *= x;
		n >>= 1;
	}
	return res;
}
//mint_baseの階乗計算
//O(x)時間が必要のため、fact_set関数を推奨する。
template<uint64_t MOD>constexpr mint_base<MOD> fact(mint_base<MOD> x)noexcept
{
	mint_base<MOD> res(1);
	for (uint64_t i = 1; i <= (uint64_t)x; ++i)
	{
		res *= i;
	}
	return res;
}
//mint_baseの階乗計算
//0からxまでの階乗を返す
//O(x)時間が必要
template<uint64_t MOD>std::vector<mint_base<MOD>> fact_set(mint_base<MOD> x = mint_base<MOD>(-1))
{
	mint_base<MOD> res(1);
	std::vector<mint_base<MOD>> set((uint64_t)(x)+1);
	set[0] = 1;
	for (uint64_t i = 1; i <= (uint64_t)x; ++i)
	{
		res *= i;
		set[i] = res;
	}
	return res;
}
//mint_base型のstreamへの出力
template<uint64_t MOD> std::ostream& operator<<(std::ostream& os, mint_base<MOD> i)
{
	os << (uint64_t)i;
	return os;
}
//mint_base型のstreamからの入力
template<uint64_t MOD> std::istream& operator >> (std::istream& is, mint_base<MOD>& i)
{
	uint64_t tmp;
	is >> tmp;
	i = tmp;
	return is;
}
typedef mint_base<1000000007> mint;
namespace mint_literal {
	constexpr mint operator""_mi(unsigned long long x)noexcept {
		return mint(x);
	}
}
using namespace mint_literal;
//mint_baseの階乗計算
//0からxまでの階乗を返す
//O(N)
template<int32_t X, uint64_t MOD = mint::MOD>
/*constexpr*/ std::array<mint_base<MOD>, X + 1> fact_set_c()
{
	mint_base<MOD> res(1);
	std::array<mint_base<MOD>, X + 1> set;
	set[0] = 1;
	for (int32_t i = 1; i <= X; ++i)
	{
		res *= i;
		set[i] = res;
	}
	return set;
}
template<typename RET = mint, typename Integral>
RET combination(Integral all, Integral get)
{
	assert(all >= get);
	get = std::min(all - get, get);
#if 1
	//時間計算量O(1)+初期化O(NlogMOD)
	static const auto fact_v = fact_set_c<要素数+1>();
	static const auto fact_div_v = [&]() {
		auto tmp = fact_v;
		for (auto& i : tmp) { i = ~i; }
		return tmp;
	}();
	//return fact_v[all] / (fact_v[get] * fact_v[all - get]);
	return fact_v[all] * fact_div_v[get] * fact_div_v[all - get];
#elif 0
	//時間計算量O(1)
	//空間計算量、初期化時間計算量O(N^2)
	constexpr int32_t ALL_MAX = 要素数;// 10'000;
	static std::vector<RET> DP_comb[ALL_MAX + 1];
	if (!DP_comb[all].empty())
	{
		return DP_comb[all][get];
	}

	if (DP_comb[0].empty())
	{
		DP_comb[0].resize(1);
		DP_comb[0][0] = (RET)1;
		DP_comb[1].resize(1);
		DP_comb[1][0] = (RET)1;
	}
	for (int32_t i = 2; i <= all; i++)
	{
		if (DP_comb[i].empty())
		{
			int32_t size = i / 2 + 1;
			DP_comb[i].resize(size);
			DP_comb[i][0] = (RET)1;
			for (int32_t j = 1; j < size - 1; j++)
			{
				DP_comb[i][j] = DP_comb[i - 1][j - 1] + DP_comb[i - 1][j];
			}
			DP_comb[i][size - 1] = DP_comb[i - 1][size - 2] + DP_comb[i - 1][(i & 1) ? (size - 1) : (size - 2)];
		}
	}
	return DP_comb[all][get];
#else
	//時間計算量O(get * logMOD)
	RET ret = (RET)1;
	for (Integral i = 1; i <= get; ++i)
	{
		ret *= all + 1 - i;
		ret /= i;
	}
	return ret;
#endif
}

N次元配列を同じ値で埋めるテンプレ

N次元配列を同じ値で埋めるテンプレを置いておきます

template<typename T>
void fill_all(T& arr, const T& v) {
	arr = v;
}
template<typename ARR, typename U>
void fill_all(ARR& arr, const U& v) {
	for (auto& i : arr) { fill_all(i, v); }
}

使用例

int dp[1000][1000];
fill_all(dp,-1);//dpの全ての要素に-1を代入
long long dp[1000][1000];
fill_all(dp,(long long)-1);//型が同じでないといけないのでキャストしてください

キャストが要らないバージョン:N次元配列を同じ値で埋めるテンプレ2 - 永夜の記録