複素数による平面幾何プログラミング

C++

ICPCが近いので、ライブラリを作らないといかん…。
というわけで、平面幾何。
平面幾何を扱うのにcomplexを使うのはもはや常識(?)
なのかもしれないが、そこをちょっとおさらいしてみることにする。
最初このテクニックを見たとき、
自分でstruct point{ int x,y; }; とか、
演算子の定義とかを書く必要がないのか、と思ったのだが
私の考えが浅かった。複素数で点を扱うと、
複素数平面上で平面幾何を扱うことになり、
便利な関数がたくさん利用できるということになる。
さらに、回転が簡単にかけたりするので、問題を回転させて
より簡単な問題に帰着したりできるかもしれない。

  • 加算、スカラ倍
typedef complex<double> pt;

pt x,y;
pt z=x+y;
double c;
pt w=c*z;

このようなことが何もせずに使える。

  • 二点間の距離
abs(x-y);
  • ベクトルのなす角

ベクトルvのx軸となす角は、

arg(v);

ベクトルuとvのなす角は

arg(u/v);

と書ける。


ベクトルuをベクトルvに移す写像fを

pt f=polar(abs(u)/abs(v),arg(u/v));

簡潔に書ける。

ベクトルu,vの内積を

real(u*conj(v));

で求めることができる。
(成分を計算することにより簡単に確かめられる)

ベクトルu,vに対して行列[u,v]の行列式を

imag(conj(u)*v);

で求めることができる。
(同上)

  • 三角形の符号付面積
double area(pt a,pt b,pt c)
{
  return imag(conj(b-a)*(c-a))*0.5;
}

点aを原点に平行移動させて、行列式より計算できる。

  • 二直線の交点

直線 a+su,b+tv (a,b,u,v はベクトル、s,tはスカラ変数)
の交点。(s,tを返す)

pair<double,double> cross_lines(pt a,pt u,pt b,pt v)
{
  double det=imag(conj(u)*v);
  if (abs(det)<epsilon) throw "交差しない";
  return make_pair(
    -imag(conj(v)*(b-a))/det,
    -imag(conj(u)*(b-a))/det);
}

証明は…成分を計算すると確かに正しいのだが、
式の意味するところは良く分からんです…。
交差しないときに例外を投げてしまっているが、
これはあんまり良くないので、実際に使うときは
うまいインターフェースを考えよう。

  • 円と直線の交点

直線 x=a+su と円 |x-c|=r の交点を求める。
(sを二つ返す)

pair<double,double> cross_circle_and_line(pt a,pt u,pt c,double r)
{
  pt ac=a-c;
  double b=real(u*conj(ac));
  double d=b*b-norm(u)*(norm(ac)-r*r);
  if (d<0) throw "交差しない";
  return make_pair(
    (-b+sqrt(d))/norm(u),
    (-b-sqrt(d))/norm(u));
}

単純に直線と円の式を解くと上の式は得られる。

  • 二円の交点

円 |x-c|=r と |x-d|=s の交点を求める。

pair<pt,pt> cross_circles(pt c,double r,pt d,double s)
{
  double di=abs(c-d);
  if (di>r+s || di<abs(r-s))
    throw "交差しない";
  double l=((r*r-s*s)/di+di)/2.0;
  double n=sqrt(r*r-l*l);
  pt e=c+(d-c)*l/di;
  pt p=(d-c)*n/di*pt(0,-1);
  return make_pair(e+p,e-p);
}

二つの円の交点と、二つの円の中心のなす凧形を考えれば
式が導出できるが、ややこしいので省略。
複素数ならではなことはできていないが、
短くする効果は十分だろう。