読者です 読者をやめる 読者になる 読者になる

ABC033 D-三角形の分類

 この問題
 http://abc033.contest.atcoder.jp/tasks/abc033_d
 xy平面上の座標がN個与えられて、そのうち三点を選んで作れる三角形\begin{eqnarray}{}_N C _3\end{eqnarray}個のうちの鋭角三角形、直角三角形、鈍角三角形の個数を答える問題。愚直解は自明(というか何しても解ける)ですが、まともな速度のものを書こうとすると少し難しいです。というか全然分かりませんでした……orz。ということで色々メモしておきます。
 まずはatan2関数の使い方メモ。
 atan2関数→atan2(b, a)でxy平面上の点(a, b)と原点(0, 0)がなす角を求めることができます。これを偏角と呼び、極座標(r, \theta)\thetaと同じです。
 (a, b)極座標(r, \theta)が同じ点を表すとします。このとき、a = r \cos \theta , b = r \sin \thetaが成り立ちますから、\frac{b}{a} = \tan \thetaです。このとき、tangent逆関数arctanとすると\theta = \arctan \frac{b}{a} と表記できます。この\thetaを勝手に計算してくれるのがatan関数で、\frac{b}{a}を引数に持つのがatan、bとaの二つの引数を持つのがatan2関数らしいです。また、戻り値は両方ともdoubleですが、値の範囲が異なり、- \frac{\pi}{2}  \leq atan \leq  \frac{\pi}{2}である一方、-\pi \leq atan2 \leq \piとなっており、atan2の方がユークリッド座標での使用には適しているようです。
 さて、解法ですが、とりあえずある点iを選び、その点と別の点jを結んだ傾きをatan2関数で偏角に変換して配列に入れ、ソートしておきます。そして、偏角が小さい方から点jを見ていきます。このとき、直線ijと直線ikがなす角が90度以上になるようなkの数d、90度になるようなkの数cをlower_boundや尺取法で求めます。そしてそのままd、cを鈍角三角形の数、直角三角形の数に足せばいいわけです。
 というのも、直角三角形に直角は一つしかありませんし、鈍角も同様です。ですからこのように全ての角をなめるように見ていけば鈍角/直角三角形の数がわかります。そして総数から引けば鋭角三角形の数もわかると。
 なかなか頭いい感じの解法ですね。尺とりが苦手なのでまずはlower_boundで実装。lbは二分探索なのでオーダーはO(n^2 \log n)くらいでしょうか。
 実装で注意する点がいくつかあります。
 ある直線ijについて、直線ijと直線ikがなす角が180度以下になるような最小のkを求める場合を考えましょう。この時、ij偏角が大きく、たとえば3.1などだった場合はkは存在しないことになってしまいます。しかし、偏角(弧度法、小数で表記)が-0.1であるikがあれば、ijとの間の角は177.6度程度で、kは存在します。このようなケースに備えるため、偏角をいれてあるソート済み配列に偏角+360度した角を入れ、偏角を二周分持つことが必要です。実際は0以上の値で一番小さい偏角までをいれれば十分なのですが、面倒なので全部入れてしまっています。全く知りませんでしたがわりと普通のテクらしいです。また、lower_boundをintにキャストするには

int x = lower_bound(hoge.begin(), hoge.end(), x) - hoge.begin();

で大丈夫っぽい。あとは90度付近の誤差丸め込みですけどこれはそんなに難しくないかな。解説等を見ながら実装しました。

#include <bits/stdc++.h>
using namespace std;
typedef long long int ll;
const double EPS = 1e-13, pi = 3.1415926535897932384626;
#define rep(i,n) for(int i=0;i<(int)(n);i++)
int x[2001], y[2001];
int main(){
    int n;
    scanf("%d", &n);
    rep(i, n) scanf("%d %d", x + i, y + i);
    ll ccnt = 0, dcnt = 0;//ccntが直角、dcntが鈍角
    rep(i, n){
        vector <double> da;//drift angle(偏角)
        rep(j, n){
            if(j == i) continue;
            da.push_back(atan2(y[j] - y[i], x[j] - x[i]));
        }
        sort(da.begin(), da.end());
        rep(j, n - 1) da.push_back(da[j] + pi * 2);//円環を二周ぶん持つ
        rep(j, n - 1){
            int x = lower_bound(da.begin(), da.end(), da[j] + pi / 2 - EPS) - da.begin();
            int y = upper_bound(da.begin(), da.end(), da[j] + pi / 2 + EPS) - da.begin();
            int z = lower_bound(da.begin(), da.end(), da[j] + pi) - da.begin();
            ccnt += y - x;
            dcnt += z - y;
        }
    }
    ll ecnt = (ll)n * (n - 1) * (n - 2) / 6 - (ccnt + dcnt);
    printf("%lld %lld %lld\n", ecnt, ccnt, dcnt);
    return 0;
}

内積で尺とりも試しましたがバグって無理でした……orz