Senの競技プログラミング備忘録

こけた問題を自分用の解説で載せる。けんちょんさんのブログを目指したい。質的にも量的にも。こけた問題だけに限定するけど

ARC155 - D - Pairs

atcoder.jp

概要

長さNの数列Aが与えられ、それらの中の2つの要素の組はN*(N-1)/2個あるが、それら積のうち、小さい順にK番にあたる値を求めよ。

制約: 1 <= N <= 2 * 10^5, -10^9 <= A_i <= 10^9

ACするまで

1. 億マス計算じゃん。(直前のバチャで億マス計算を難なくACして勝ったと確信)
2. え・・・要素マイナスもあるし、同じ要素同士の組選べないからいろいろ除外しないといけないのか・・・。
3. とりあえず書いてみるも、なかなか合わずにEに行く。ちなみにEもDPだと早いところで気づいたくせに最後まで正しい実装できなかった。
4. コンテスト終了
5. 2019/2/16終了
6. 朝起きて、一から書き直してみるもサンプル3で死ぬ。
7. 1時間にわたる解析の結果、どうやらlong doubleでも演算の精度が足りないことがバグの原因らしい。
8. 除算を使わないように書き換えて、AC!

ACするまで

この問題は、Atcoderの過去問にあたる億マス計算という問題とほぼ同じであるが、この問題では配列の要素に負や0も含まれていて、かつ同じ要素2つによるペアは禁止というさらに厳しい制約になっている。

根幹をなす考え方へ至るのは億マス計算の解説スライドを見るとわかりやすいので、ここではその結論を述べる。


f(x):=x以下の要素がK個以上存在するか?
を定めて、そのxで二分探索をする。

さて、判定関数のf(x)はどうすればいいのか?今回はxもA[i]も負、ゼロになりうるので、細心の注意を払って、考えていこう。

x>=0 and A[i]>0

これは、一番スタンダードなものであり、億マス計算のような書き方そのままでもよい。つまり、昇順に並んだAから、x / A[i]を切り捨てたもの(以後floor()関数と呼ぶ)以下の最右端の場所を二分探索する。それまでのAの要素は、A[i]とかけ合わせると必ずx以下になるので、x以下のペアの個数の増分は二分探索からも止まった位置からわかる。
ただし、ここで同じものを選んでペアを作ってはいけない、というルールがある。上記の二分探索でfloor(div / A[i]) >= A[i]ならば、(A[i], A[i])というペアもカウントしたことになるので、1つ分引く必要がある。

x >= 0 and A[i]==0

この時、A[i]と自身を除くすべてのペアはすべて積が0になるため、条件を満たす。よって、増分はN-1である。

x>=0 and A[i]<0

例として、x=7, A[i]=-3の場合を考える。7/-3=-2.333より、Aのうち、-2.333より大きい要素ならば、すべてx以下になる。
つまり、切り上げ関数roundを使って、round(x / A[i])と書くと、降順に並んだAのうち、round(x / A[i])以上のものの最右端を二分探索して、その場合計算すればよい。
この時、重複カウントについては、round(x / A[i]) <= A[i]であれば、重複カウントしてるということになるので、1つ引く必要がある。

x<0 and A[i]>0

例として、x = -5, A[i]=2の場合を考える。-5/2=-2.5より、Aのうち、-2.5以下のものであれば、x=-5以下のペアを作れる。
一般化すると、x>=0 and A[i]>0の場合と同じく、昇順に並んだAから、x / A[i]を切り捨てたもの(以後floor()関数と呼ぶ)以下の最右端の場所を二分探索する。
この条件下での重複を考えてみる。xは負で、A[i]は正なので、割った値は負であるが、対象になる区間の右端が負であるのにA[i]は正なので、重複することはあり得ない。

x<0 and A[i]==0

この時、どのようにペアを作っても0になるが、いずれもx以下ではないので、増分には寄与しない。

x<0 and A[i]<0

例として、x=-7, A[i]=-3である。-7/-3=2.333であり、これはAの要素のうち、2.333以上のものと-3がペアを作れば、必ず-7以下になるということ。
一般化すると、x>=0 and A[i]<0の時と同じように、降順に並んだAのうち、round(x / A[i])以上のものの最右端を二分探索して、その場合計算すればよい。
この条件下でも、条件を満たすのはx / a[i]以上であるが、x / a[i]は正であり、a[i]は負なので、重複することはあり得ない。

ここまでの考察をコードに落とし込むと、次のようなものになる。

typedef long long ll;
//f(x) := x以下のペアはK個以上存在するのか?
//aは昇順、revaは降順にソート済
bool f(ll key) {
	ll cnt = 0;
	for (int i = 0; i < N; i++) {
		if (key >= 0 && a[i] > 0) {
			cnt += lower_bound(a, a + N, floor((long double)key / (long double)a[i]) + 1) - a;
			if (floor((long double)key / (long double)a[i]) >= (long double)a[i])cnt--;
			/*
			if(key == 448283280358331041 && i == 24){
			cout << round((long double)key / (long double)a[i]) << endl;
			}
			*/
		}
		else if (key >= 0 && a[i] == 0) {
			cnt += N - 1;
		}
		else if (key >= 0 && a[i] < 0) {
			cnt += lower_bound(reva, reva + N, ceil((long double)key / (long double)a[i]) - 1, greater<ll>()) - reva;
			if (ceil((long double)key / (long double)a[i]) <= (long double)a[i])cnt--;
			
		}
		else if (key < 0 && a[i] > 0) {
			cnt += lower_bound(a, a + N, floor((long double)key / (long double)a[i]) + 1) - a;
			//重複することはない
		}
		else if (key < 0 && a[i] == 0) {
			cnt += 0;
		}
		else if (key < 0 && a[i] < 0) {
			cnt += lower_bound(reva, reva + N, ceil((long double)key / (long double)a[i]) - 1, greater<ll>()) - reva;
			//重複することはない
		}
	}
	return cnt >= 2 * K;
}

しかし、この判定関数はサンプルケース3で落ちる。

僕の手元では、448283280358331064と表示されるはずが、最後の2ケタだけが違う448283280358331041となるのである。
がんばって解析した結果、どうやら、(0-idxで)A[24]==537203564の時に、上の関数は1ペア分多くカウントするらしい。これは次のが原因らしい。

x == 448283280358331064, A[i] == 537203564の時に、
x / A[i] ==
834475625.99999995718......
であるが、どうやら、floor((long double)x / (long double)A[i])の結果は、
834475626
である。そして、A[27] == 834475626であるため、これを誤ってカウントするみたいである。(余談だが、A[24] * A[27]がサンプル3の答え)
つまり、残念ながら、long doubleでさえ、10^18 / 10^9オーダーの計算では、上のような実装法では精度が足りないということである。(long doubleで通してる人もいたので、やり方を改めればできると思われるが、不安要素を取り除くために下のような書き方がいいと思う。値を埋め込んで計算したら大丈夫だったが、変数に代入するとダメだったりと、実験結果は不安定だとわかる。)

これを回避するには、除法を使わない実装にするしかない。

これは、既存のlower_boundを使わずに、自分で二分探索を書くしかない。僕はめぐる式で次のように実装した。
また、上の6つの場合分けも、乗算にすることでx/A[i]で発生するゼロ除算がなくなり、場合分けも2つで済む。

int my_b_search(ll key, ll val, bool mode) {
	//mode -> trueならば昇順での二分探索 falseなら降順での二分探索
	int ok = mode ? -1 : N;
	int ng = mode ? N : -1;
	while (abs(ok - ng) > 1) {
		int mid = (ok + ng) / 2;
		if (val * a[mid] <= key)ok = mid;
		else ng = mid;
	}
	return ok;
}
 
bool f(ll key) {
	ll cnt = 0;
	for (int i = 0; i < N; i++) {
		if (a[i] >= 0) {
			cnt += my_b_search(key, a[i], true) + 1;
		}
		else {
			cnt += N - my_b_search(key, a[i], false);
		}
		if (key >= a[i] * a[i])cnt--;
	}
	return cnt >= 2 * K;
}

これでACすることができる。
ACソースコード
Submission #10180135 - AtCoder Beginner Contest 155

一言

Dに置くものにしては難しすぎるし、解説が正直何もわからない人にとっては意味不明だし、実装するだけと言って誤差死してるので実装例はほしかった。久しぶりにAtcoderらしくない解説だと感じた。(上から目線)
誤差死したのでレートには痛かったけど、ここで転んでよかった。

Special Thanks

Submission #10176434 - AtCoder Beginner Contest 155
このtakoyakimaster65さんのソースコードを大いに参考にした。ありがとう!
バグについて的を射たアドバイスをくれたhotman78(@hotmanww)、ありがとう!