D - Line Crossing Editorial by rsk0315

座標計算に基づく解法

公式解説 における下記の記述が突飛に感じられる人向けです。

実は、\((A_i+B_i)\bmod N\) で分類することができます。

Twitter を見るに、浮動小数点数で座標計算をしようとして誤差に苦しんでいる人もしばしばいるようなので、そういう解法(がうまくいかない理由)についても触れます。

座標計算

単位円を考え、点 \(i\)\((\cos(\frac{-2i\pi}{N}), \sin(\frac{-2i\pi}{N}))\) にあるとします。

\(i\) と点 \(j\) を通る直線の傾きは次の式で表せます。

\[ \begin{cases} \displaystyle\frac{\sin(\frac{-2i\pi}{N})-\sin(\frac{-2j\pi}{N})}{\cos(\frac{-2i\pi}{N})-\cos(\frac{-2j\pi}{N})}, & \text{if}\: \cos(\frac{-2i\pi}{N})\ne\cos(\frac{-2j\pi}{N}); \\ \infty, & \text{if}\: \cos(\frac{-2i\pi}{N}) = \cos(\frac{-2j\pi}{N}). \end{cases} \]

ここで、\(\phi\)\(\theta\)\(\cos(\phi)\ne\cos(\theta)\) なる実数とすると、和積公式から下記が成り立ちます。

\[ \begin{aligned} &\phantom{{}={}} \frac{\sin(\phi)-\sin(\theta)}{\cos(\phi)-\cos(\theta)} \\ &= \frac{2\cos(\tfrac{\phi+\theta}2)\sin(\tfrac{\phi-\theta}2)}{-2\sin(\tfrac{\phi+\theta}2)\sin(\tfrac{\phi-\theta}2)} \\ &= -\cot(\tfrac{\phi+\theta}2). \end{aligned} \]

ここで、\(\cot(x) = \frac{\cos(x)}{\sin(x)}\) です。\(x\equiv 0\pmod{\pi}\) のときはよしなに定義されているとします。

\(\cos(\phi)\ne\cos(\theta)\) より \(\phi\not\equiv\theta\pmod{2\pi}\) かつ \(\phi\not\equiv-\theta\pmod{2\pi}\) なので、\(\sin(\tfrac{\phi+\theta}2)\ne0\) かつ \(\sin(\tfrac{\phi-\theta}2)\ne0\) です。

\(\phi = \frac{-2i\pi}{N}\) かつ \(\theta = \frac{-2j\pi}{N}\) とすると、

\[ \begin{aligned} -\cot(\tfrac{\phi+\theta}2) &= -\cot(\tfrac12\cdot(\tfrac{-2i\pi}{N}+\tfrac{-2j\pi}{N})) \\ &= \cot(\tfrac{(i+j)\pi}N) \end{aligned} \]

となります。また、\(\cot(0)=-\infty\) とすると、 \(\phi\equiv-\theta\pmod{2\pi}\) のとき \(-\cot(\tfrac{\phi+\theta}2)=\infty\) となります。 任意の \(1\le i\lt j\le N\) について \(-\tfrac{2i\pi}N\not\equiv-\tfrac{2j\pi}N \pmod{2\pi}\) ですから、\(\phi\equiv\theta\pmod{2\pi}\) のケースについて考える必要はありません。

平行性の判定

上記より、点 \(i\) と点 \(j\) を通る直線の傾きは \(-\cot(\tfrac{(i+j)\pi}N)\) と表せることになります。

\(\cot\) は周期 \(\pi\) の関数なので、\(\tfrac{(i+j)\pi}N \bmod \pi\) で分類すればよく、すなわち \((i+j)\bmod N\) で分類すればよいことがわかります。

誤差評価

\(\gdef\roundcirc#1{\operatorname{\circ}{({#1})}}\)

実数 \(x\) を浮動小数点数で表せるように丸めた値を \(\roundcirc{x}\) と書きます。型としては binary64、丸めモードとしては tiesToEven を考えます(いわゆる double のデフォルトの挙動)。 また、浮動小数点数の四則演算を \(\oplus\), \(\ominus\), \(\otimes\), \(\oslash\) と書きます。

\(\cot(x)\) が用意されていない言語が多そうなので、以降は \(\tan(x)\) を考えます。tan\(\roundcirc{\tan(x)}\) と異なる値を返す実装もしばしばありますが、一旦は \(\roundcirc{\tan(x)}\) を返す(いわゆる correctly-rounded な値を返す)とします。

note: 精度の都合で \(\tan(x)\) を正確に表せないという話ではなく、\(\tan(x)\) に最も近い浮動小数点数以外を返す場合があるという話です。

note: GNU C Library における誤差は Known Maximum Errors in Math Functions などに記載があります。

下記の例が示すように \(\tan(x) = \tan(x+\pi)\) に相当する性質を満たさないため、前項の考察に基づいて素朴に 1.0 / tan((i + j) * PI / n) と実装しても正解することはできないでしょう。

\[ \begin{aligned} \roundcirc{\pi} &= 7074237752028440\times 2^{-51}, \\ \roundcirc{\tan(\roundcirc{\pi}\oslash 3)} &= 7800463371553960 \times 2^{-52}, \\ \roundcirc{\tan(\roundcirc{\pi}\otimes 4\oslash 3)} &= 7800463371553953 \times 2^{-52}. \end{aligned} \]

余談

なお、IEEE 754 に記載のある tanPi と呼ばれる関数は \(x\) に対し \(\roundcirc{\tan(\pi x)}\) を返すものですが、\((i+N)\oslash N = (i\oslash N)+1\) とは限らないため、これを用いてもうまくいかないと思われます。これを提供している言語があるかはわかりません。

\[ \begin{aligned} 1\oslash 3 &= 6004799503160661\times 2^{-54}, \\ 4\oslash 3 &= 6004799503160661\times 2^{-52} \\ &= (1\oslash 3) + 1 - 2^{-54}, \\ \mathrm{tanPi}(1\oslash 3) &= 7800463371553961\times 2^{-52}, \\ \mathrm{tanPi}(4\oslash 3) &= 7800463371553958\times 2^{-52}. \end{aligned} \]

もちろん、和積公式を用いた変換をする前の式を用いる場合も、さまざまな誤差が生じるため、うまくいかないでしょう。一般に、誤差を含む値同士の減算は避けるべきです。

補足

\((A_i+B_i)\bmod N\) で分類できるのはなぜか?という問いへの答えにはなっていると思いますが、コンテスト中にどうして思いつけるのか?という問いへの答えにはなっていないと思います。実際、自分はコンテスト本番は座標計算を介さずに解きました。

座標計算で解こうとした人は、焦っていきなり「非自明な判定法を思いつけるようになりたい」と思うのではなく、式変形を通して導出できるようになると後々役に立つだろうと思います。

posted:
last update: