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

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

ABC191 - D - Circle Lattice Points

ABC 精度 要復習 良問

https://atcoder.jp/contests/abc191/tasks/abc191_d

概要

座標平面上で、$ (X, Y) $中心で半径が$ R $の円の内部や辺上に、何個の格子点(座標がすべて整数の点)がある？

制約：$ 0 \leq |X|, |Y|, R \leq 10 ^ 5 $。$ X, Y $は小数点以下4桁まで与えられる

考え方

  1. 円のうち、y座標で見れば最大でも$ O(Y) $通りのy座標の点しか答えの候補にならない。y座標の値さえ決まれば、円内に格子点となれる点の数は簡単に二次方程式を解くことで数えられる。
  2. doubleで実装。->WA　doubleでは精度が足りない。これを補う工夫が必要。
  3. long long なら$ 10 ^ 18 $の精度で計算できる。今回はこれならギリギリいけそう。
  4. 入力をlong longに直して(ここ罠あり)、自前でsqrtを書く。
  5. AC!

考察

この問題の解法の大筋は1通りであるが、精度を解決するためには2通りのやり方がある。

  1. $ 10 ^ 4 $倍して、long longで計算する。
  2. long double で精度の問題を解決し、辺上ギリギリの点を数えるために$ R $にEPSを足すことで考える。

アルゴリズム

この問題を解くアルゴリズム自体は平易である。格子点の数を数えるやり方として、「円の内部や辺上に含まれている格子点のy座標はたかだか$ O(Y) $通りしかない。」を利用して、y座標を固定した時に、円の内部や辺上にある格子点は、円の式$ (x - X) ^ 2 + (y - Y) ^ 2 = R ^ 2 $から二次方程式を解けば求まるとわかる。

誤差の問題を考えなければ、√の計算に$ O(\log |Y|) $かかると考えられるので、全体で$ O(|Y| \log |Y|) $程度の計算量で収まる。

精度の問題

しかし、ここで考えなければいけないのはdoubleの精度で足りるか？である。 doubleは倍精度浮動小数点であり、これの計算機イプシロン(doubleで表せる1以上で最小の数と1の差)は$ 2.220 \times 10 ^ {-16} $である。

今回の場合、ルートで最大$ (10 ^ 5) ^ 2 - (10 ^ {-4} ^ 2) $の計算をすることになる。この時、要求される精度は$ 10 ^ 18 $レベルであり、明らかにdoubleの計算機イプシロンを超えている。　つまり、doubleによる計算では今回の場合精度が足りないのである(実際多くの人はここで落ちたがゆえに青上位の難易度になっている)

精度問題の回避手段(方法1)

では、どのように精度問題を回避するのか？ここで着目したいのは、doubleの分解能?は$ 2.220 \times 10 ^ {-16} $レベルであるが、long long は$ 9.0 \times 10 ^ {18} $レベルの精度の演算ができる所である。つまり、64bit整数として計算してあげれば精度面ではなんとか(ギリギリだが)なりそうということだ。

というわけで、まず入力の$ X, Y, R $をすべて$ 10 ^ 4 $倍して、整数としよう。この時、実は以下のように書くとWAになる。

long long mX = X * 10000.0;

ここでも多くの人がつまずいた。そのまま代入による切り捨てや、ceil()は、負の数ではうまく働かなかったり、$ X \times 10000.0 $の値自体を床関数で計算しても綺麗に求まらないことになる。

実際に、以下の場合でエラーが生じる。

99999.9999 * 10000.0
99999.9999は倍精度浮動小数点で表すと、1.525878904724121 * 2^{16}であり、
(Wolfram Alphaで計算すると)、999999998.999...となってしまう。

このような状況では、roundll()で四捨五入するべきである。 そもそも浮動小数点では多くの値は循環小数となり、綺麗に表せない。個の誤差は普段とても小さく、(coutで出力すれば)見た目上は綺麗に表せているものの、計算を重ねることによって誤差が蓄積して最終的な値はずれてしまう現象が起きる。この時、誤差が小さいうちに毎回roundなどの処置を考えるべきである

long long mX = round(X * 10000.00);

二乗根の計算(方法1)

STLのsqrtは、doubleなどの浮動小数点仕様である。この問題では、doubleは精度不足から、std::sqrtを用いることができない。

そこで、$ \sqrt{x} $は単調増加であることを利用して、二分探索で求めることができる。具体的には以下の通りである。

ll my_sqrt(ll x) {
    ll ok = 0, ng = 2 * 1e9;
    while (ng - ok > 1) {
        ll mid = (ok + ng) / 2;
        if (mid * mid <= x)
            ok = mid;
        else
            ng = mid;
    }
    return ok;
}

long doubleで精度問題を対処(方法2)

long double(4倍精度浮動小数点)の計算機イプシロンは$ 1.926 \times 10 ^ {-34} $であり、今回の計算の場合ではこれで精度が事足りるとわかる。

この時、計算幾何によるあるテクニックとして、辺上にあるものの数え上げは、わずかに大きい外形の内部で考える、がある。 例えば、今回の場合は半径$ R $ではなく、ある十分小さな値$ \epsilon $を用意して、$ R + \epsilon $が半径の円の内部を代わりに考える。 この操作を施せば、辺上にあるものは安全に数え上げることができる。

今回の場合、一番精度が危ういようなケースは以下の通り。

$$ X = 0, Y = 100000, R = 100000, x = 0.0001, y = 0 $$

この場合は、ギリギリ入ってないが、$ R $の代わりに$ (R + \epsilon) $を用いて、ギリギリ入ってしまう場合は、 $ 10 ^ {10} + 10 ^ {-8} = (10 ^ 5 + \epsilon) ^ 2 $であり、解くと$ \epsilon = 5.0 \times 10 ^ -{14} $となる。これを超えるとこのようなエラーが生じるので、 $ \epsilon < 4.0 \times 10 ^ {-14} $とするのが無難である。

このように、計算幾何での$ \epsilon $はこのように最悪ケースでエラーする限界を考えて決めるとわかりやすい。

ACソースコード

解法1

C++　自前で10000倍にceilやroundを実装した。

解法2

C++　EPSを$ 10 ^ {-14} $にした。$ 4.0 \times 10 ^ {-14} $でも落ちるみたいなので、安全を取ろう。

最後に

コンピューターアーキテクチャや開発よりの所を問われた。こういうところは実務でも生きてくるので落としちゃダメだ。