けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

TopCoder SRM 400 DIV1 Hard - CollectingBonuses (本番 2 人)

SRMDIV1Hard 期待値 回数の期待値 確率 誤差 式変形 調和級数 テク:スタートを0としてよい 期待値の線形性 TopCoder

ついこの間の AGC でもあった「〜が終了するまでの回数の期待値」に関する問題！！！

それにしても数値誤差しゃん...
150 人が提出してAC 2 人という恐ろしい問題。
　
実はこの回の writer であった ymatux さんから、直接話を聞いたのだった。Petr はこの回 resubmit をしていて、148 人が WA った中 AC を勝ち取ったらしい。「やっぱり Petr は違うなーと思った」と言っていた。

問題概要

n 種類の品目が当確率で排出されるガチャがある。このガチャで k 種類の異なる品目を揃えるまでに引く回数の期待値を、相対誤差または絶対誤差が 10^{-9} 以内の範囲で求めよ。

制約

  • 1 \le k \le n \le 10^{18}

考えたこと

答えの数式を得ることはそれほど難しくない

  • 1 個目を得るまでの期待値: 1 回
  • 2 個目を得るまでの期待値: \frac{n}{n-1} 回
  • 3 個目を得るまでの期待値: \frac{n}{n-2} 回
  • ...
  • k 個目を得るまでの期待値: \frac{n}{n-k+1} 回

というわけなのでこれを合計すればよい。もっというと

  • H(n) = 1 + \frac{1}{2} + \dots + \frac{1}{n}

として

  • n(H(n) - H(n-k))

が答えである。難しいのは「 n \le 10^{18}」という制約と「数値誤差」。数値誤差の怖さに気づかずに大量の Failed System Error が出たようだ。

H(n) の計算: n <= 1018 なので

n がでかいので愚直に全部足すことはできない。なので解析的になんとかしよう。

f(n) はいわゆる調和級数とよばれるものなので、解析的な取り扱いができそうである。数IIIで習った区分求積法によって、以下のことがいえる。

  • f(n) 〜 \log{n} + \gamma

ここで \gamma はオイラー定数とよばれる値で、0.57... とかなんとか。これの収束性が気になるが、以下の記述がある！

これによるとなんと、

  • f(n) 〜 \log(n + \frac{1}{2}) + \gamma

とした方が近似精度がよい！！！そして二乗の精度だとわかる。なので 10^{-9} の精度で求めたかったら、

  • n \le 10^{5} くらいまでは自力で計算
  • n \gt 10^{5} では \log{n + \frac{1}{2}} + \gamma を計算

という感じで H(n) を求めることができる

まだ罠があるっぽい！！！

ところがこれでもまだ罠があって「そんなん気づくかーーー！！！」という感じなのだけれども、 　

n(H(n) - H(n-k)) 　

が答えだけど、 H(n) - H(n-k) を安易に計算してはいけない！！！一般に数値計算において値が近い浮動小数同士の引き算は超危険で、たとえば 　

1.3244324434... - 1.3244324431... = 0.00000000003... 　

となって、相対精度がどっかーんと落ちてしまう。これを回避するために、 k が十分小さいときは 　

H(n) - H(n-k) = \log(\frac{2n + 1}{2n - 2k + 1}) 　

とすることが考えられる。 n-k がある程度大きければこう考えて良いと思われるので作戦としては

  • \frac{1}{n-k+1} + \frac{1}{n-k+2} + \dots + \frac{1}{n} の計算の先頭部分について、分母が十分大きくなるまで (安全のため 10^{7} くらい) やる
  • 残りは、 H(n) - H(n-k) = \log(\frac{2n + 1}{2n - 2k + 1}) を活用

という感じでやろう。

まだまだ罠があるっぽい！！！

ここまでケアしてもまだ罠があって、 \log(\frac{2n + 1}{2n - 2k + 1}) の \log の中身は 1 に近いので、このままだと精度が低い。そこで、log1p を用いるとよい。これは \log(1+x) を計算してくれる関数。 というわけで 　 H(n) - H(m) = \log(1 + \frac{2(n-m)}{2m+1}) = {\rm log1p}(\frac{2(n-m)}{2m+1}) 　

ということになる。

#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <cmath>
using namespace std;
#define COUT(x) cout << #x << " = " << (x) << " (L" << __LINE__ << ")" << endl

// n/n + n/(n-1) + n/(n-2) + ... + n/(n-k+1)
// f(n) = 1 + 1/2 + ... + 1/n として
// n(f(n) - f(n-k)) が答え

class CollectingBonuses {
public:
    double expectedBuy(string sn, string sk) {
        long double res = 0.0;
        const long long M = 100000000;

        istringstream isn(sn); long long n; isn >> n;
        istringstream isk(sk); long long k; isk >> k;

        long long m = n - k + 1;
        while (m < M) {
            res += (long double)(1.0) / m;
            if (m == n) return res * n;
            ++m;
        }
        --m;

        // f(n) - f(m) 〜 log((2n+1)/(2m+1)) = log(1 + 2(n-m)/(2m+1))
        res += log1pl((n - m) * 2.0 / (m * 2.0 + 1));
        return res * n;
    }
};