１．計算幾何の基本　　　シリーズ：「基礎的計算幾何ライブラリの作り方」

シリーズの構成

１．計算幾何の基本　　　←いまここ
２．計算幾何のインフラ整備
 ３．直線、線分の計算幾何
４．円の計算幾何
５．多角形の計算幾何
６．AOJにある問題の解説

シリーズの構成
前言
対象者
シリーズの構成(再掲)
計算幾何とは&要点
誤差
誤差による計算幾何の実装上の影響
計算幾何に必要な数学
最後に
参考文献

前言

　外は新型コロナウイルスで騒がしく、様々なイベントが順延して、果てには大学も授業開始日を遅らせてきました。(2020/3/10現在)このまとまった期間を利用して、幾何ライブラリを整備していましたが、幾何特有の実装や手法などはクセものであり、界隈に落ちてる先人の足跡を一つ一つたどるという作業は、かなり大変でした。すでに完成されてるライブラリというのは多いが、まとまった解説をしてるような文献は少なく、さらに幾何特有の実装方法まで踏み込んだものは決して多くありませんでした。

　そこで、自分の勉強ついでに、競プロ界隈のため、将来の自分のために、本ブログ初のシリーズもの、シリーズ：「基礎的計算幾何ライブラリの作り方」を書きます。

　願わくば、本シリーズがこれから計算幾何を学ぶ人々にとって、より計算幾何の基礎を学ぶ敷居を下げ、魅力を伝えられるようなものになり、そして僕のレートをちゃんと上げてくれるものでありたいです。

　最後に、このシリーズを作るにあたって、多くの助言を与えてくれたshiroさん(@shiro537)、ライブラリの現物を大変参考にしていただいたこたつがめさん(@kotatsugame_t)、なにより信じられないほどのクオリティのPDFを残してくれた秀郁未さん、大変ありがとうございました！

対象者

　計算幾何について、何も知らない人でも読めるような内容です。必要な数学力としては、高校の数2Bぐらいまであれば理解できます。

　本シリーズの主な目標としては、2次元において、以下のようなものが計算できることです。

直線や線分の交点を求められる
円と直線、線分の交点を求められる。
凸包やその直径が求められる。
多角形の面積が求められる。

詳細はこちら
幾何ライブラリの構造.md - Google ドライブ

　以上の内容を、今回のシリーズでは「基礎的計算幾何」と定義して、これらについて解説をする記事になります。

シリーズの構成(再掲)

　トップの再掲になりますが、次のようにします。

１．計算幾何の基本　　　←いまここ
２．計算幾何のインフラ整備
３．直線、線分の計算幾何
４．円の計算幾何
５．凸包、多角形の計算幾何
６．AOJにある問題の解説

　AOJ(Aizu Online Judge)には、様々な分野の典型問題を重点的に学べるコースが存在しています。今回は、それらのコースによるACという検証のもとで、シリーズを構成していきます。

　今回、「１．計算幾何の基本」は

計算幾何とは何なのか
計算幾何にまつわる誤差の話
計算幾何特有の実装上の書き方
計算幾何に必要な数学の紹介

を紹介、解説します。

計算幾何とは&要点

計算幾何学は、大きく分けて2つの区分があります。[出典1]1つは、(古典的)幾何学の問題(中学生がやってるような幾何)を、大量のデータが存在してる時に、コンピューターを使い、できるだけ効率よく、そして正確に求めることです。2つは、コンピューターの力で現実世界にあるモノの形を、デザインするための計算というものです。
　
　本シリーズでは、1つ目の区分のもとの、基礎部分を解説します。

誤差

　幾何の計算の結果のほとんどが小数になりますが、コンピューターの中ではどのように小数を保持してるのでしょうか？
　
　C++では、doubleという型で保持しますが、これは「浮動小数点数」という型で、小数を次のように変換して、保持します。
　
　 $(符号つき)a*10^b$ 　この符号と $aとb$ を保持するわけです。
　
　これについては、以下のサイトの最初で解説されてます。(京都産業大学の教材)
www.cc.kyoto-su.ac.jp

　この中には次のような記述があります。

2.225074 10^-308 < double の絶対値 < 1.797693 10^+308

　こんだけ大きい値を持ててしまったのならば、誤差なんてなくないですか？

　しかし、残念ながらあります。

計算機イプシロン - Wikipedia

　この上の記事にもあるように、「1の次大きく表現できる小数と1の差」というのを計算機イプシロンと定めていて、doubleの計算機イプシロンは $2.220 \times 10^{-16}$ 。

　このことからもわかるように、doubleは大きい値を持つだけなら持てるけど、その一定桁より下の細かい桁の情報は消えてしまう、といいうわけです。

　ゆえに、計算幾何を扱う競技プログラミングの問題では、必ず許容誤差が与えられ、絶対誤差がその値以内ならば正解とみなされます。そして、いかに誤差の少ない計算で求めるか？という実装も必要です。

誤差による計算幾何の実装上の影響

EPS

　　計算幾何において、許容する誤差をEPSという小数の定数として、定めていくことがよくあります。競技プログラミングでは、

$10^{-7}$ 以下の誤差は許容する。

となっているのならば、念のためにEPSをこれよりも10~100倍小さい、 $10^{-9}$ などにするといいでしょう。しかし、上記の計算機イプシロンのことも踏まえて、 $10^{-16}$ に近いような値は危険です。よって、僕は $10^{-9}$ から $10^{-10}$ をお勧めします。

EPSを考慮した大小判定

EPSを考慮した大小判定は、簡単なようで、意外に面倒です。
僕は様々な先人の実装を見て、次のようにしています。[出典2]

//a > 0ならば+1, a == 0ならば0, a < 0ならば-1　を返す。　基本的にEPS込みの評価はこれで行う。
	//不等式は、加減算に直してこれに適用する。
	int sgn(const double a) {
		return (a < -EPS ? -1 : (a > EPS ? +1 : 0));
	}

　この書き方の利点としては、valのequalゾーンだけ広げて、それの正or負はそのままにしてる、というのがあります。

　これによって、次のように直感的にかけたりします。

//aが0より大きいなら
if(a > 0);//本来はこう
if(sgn(a) > 0);//sgn()関数で囲むだけ

//bが0以下なら
if(b <= 0);//本来はこう
if(sgn(b) <= 0);//やはり囲むだけ

//不等式の条件などは、次のように書く
if(a > b);
if(sgn(a - b) >= 0);

if(a <= b + c);
if(sgn(a - (b + c)) <= 0);

sqrt()の中身には気を付ける

　次にようなコードを見てください。

double a, b;
//aは計算の結果、[4.31, 7.0]で与えられる。
b = sqrt(a - 4.31);

　sqrt()の中身は1.0から0.0まで変動するので、これは人畜無害のコードのように見えます。

　しかし、実行をしてみると、環境によっては、「sqrt()に定義域外の(つまり、負の)引数が渡された」といったエラーメッセージが出るときがあります。

　理由としては、数学的には確かその4.31がaに与えられても、2進数に直す過程で、正確な値でなくなり、少し小さくなってる、という可能性があります。
現に
4.31
は僕の環境では
4.30999999999999960920150
と保存されているみたいです。

　その、微妙に小さくなってることによって、sqrt()の中身が負と判定されて撃墜、というようなことは普通に起きます。特に競技プログラミングの場合では大量の、悪意あるデータを処理するのでこういうスレスレ案件は割と発生します。

対策としては、次のように、sqrt()の中身を書き直すことです。

b = sqrt(max(0.0, a - 4.31));

計算幾何に必要な数学

直交座標平面

　2次元の場合は、垂直な線を2本引き、交わった点を0として、その上で等間隔に目盛りを振っていきます。点がどこにあるのかを、横の線(x軸と呼びます)ではどこの目盛りにあるのか、縦の線(y軸と呼びます)ではどこの目盛りにあるのかで表します。(さすがに説明不要だと思いますが・・・)

　余談ですが、この直交座標平面をデカルトがベットで寝てる時、飛んでいるハエの位置がどこにあるのかなーと考えた際にひらめいたそうです。

直線、半直線、線分

　無限に長い、まっすぐな線のことを直線といいます。

　ある点から始まる、無限に長い、まっすぐな線のことを半直線と言います。

　ある2つの点の間を結ぶまっすぐな線を線分と言います。
f:id:Sen_comp:20200311230628j:plain

ベクトル

　長さと向きを持った矢印です。始点がどこにあるのかは基本的に関係ありません。(0, 0)から(2, 1)へ伸びるベクトルも、(2, 2)から(4, 3)へ伸びるベクトルも同じものとして扱います。

　ベクトルは、始点からx軸にa、y軸にb伸びてる場合、(a,b)とも書きます。座標での点の表し方と一緒ですが。

　このベクトルという長さと向きを持った矢印を考えることで、直線、辺やそれらが作る角度というのを考えやすくなります。

　また、直交座標平面の上の点を、(0, 0)からそこへ伸びてるベクトルとみなすという考えもあり、それは位置ベクトルと言います。直交座標平面の上の位置の書き方も、ベクトルの上の表し方も同じなのは、こう見なせるでもあります。

内積

　ベクトル同士の演算(数字における掛算、割り算のようなもの)の1つです。

　(a, b)と(c, d)の内積は、「・」という記号を使い $(a, b) \cdot (c, d)$ と書きます。その値は、

$(a, b) \cdot (c, d)=ac+bd$ 　になります。

　また、内積は図にしてみると、次のような性質を持ちます。

f:id:Sen_comp:20200311230829j:plain

　演算の順番は入れ替えても同じ結果になります(式からそれはそう)。

　実をいうと内積の値自体は本シリーズではあまり重要ではないです、むしろ正か0か負かによく着目します。(でも内積の値が本質的な場所もあるので、意味はつかんだ方がいいです。)

外積

　内積と似た、別の演算です。

$(a, b) \times (c, d)=ad-bc$

　これは、(正負はさておき)2本のベクトルが作る平行四辺形の面積を表しています。つまり、2本のベクトルのなす角を $\theta$ とすると、

$\vec{A} \times \vec{B}=|A| |B| \sin \theta$ 　を満たします。

　これの正負は、次の図のようになります。

f:id:Sen_comp:20200311230923j:plain

　また、内積は順番を入れても値は変わりませんが、外積は正負が逆になります！演算の順番に気を付けよう！

　ちなみに、厳密にいうと、外積は3次元以上のベクトルでしか定義できませんし、内積と違って、かけた結果もベクトルになります。今回の2次元での定義は、3次元の時の外積の定義を無理やり2次元に当てはめたものです。

　このように、ベクトルという概念と、それに付随する演算として内積と外積を使います。内積と外積はいずれも、2つのベクトル $(a,b)$ と $(c,d)$ の四則演算で計算できるのがポイントです。

　四則演算で計算できる値は、ルートなどの演算に比べると誤差を産みづらいというメリットがあります。ですので、本シリーズの実装も、できるだけこのベクトル、内積、外積を使用したものになっています。

最後に

　次回は、2次元幾何における幾何の計算に必要な機能をどう実装する「計算幾何のインフラ整備」です。いよいよ本格的にコーディングしていきます。

参考文献

[出典1]　計算幾何学 - Wikipedia　2020/02/10 01:00アクセス
計算幾何学 - Wikipedia

[出典2]　geometry
https://www.ioi-jp.org/camp/2017/2017-sp_camp-hide.pdf

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

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