SRM 400 DIV1 Hard CollectingBonuses

SRM 近似 誤差回避 普通な数式

問題 Editorial

問題

飲料メーカーがキャンペーン中。
ジュースのボトル1つと引き換えに、nつの異なるコードのうちランダムな1つが等確率で得られる。
kつの異なるコードを得たい。
それが達成されるのに必要なボトルの数の期待値を求めよ

  • 1 ≦ k ≦ n ≦ 10^18

解答

Editorialを参考にした。
まず漸近式を解くと、答えは"n * (H(n) - H(n-k)) = n * (1/(n-k+1) + 1/(n+k+2) + ... + 1/n)"であることがわかる。ここで、H(n)は調和数。
しかしnが大きすぎるので直接定義式を計算することはできない。そこで、近似公式を使う。
ただし、小さなnでは誤差が大きいので、n+k+1〜n+k+10^7程度の項は定義式通りに計算する。
求めた項(n+k+10^7+1)をmとする。次に、残りの項"H(n) - H(m-1)"を求めたい。
"H(n) = γ + ln(n + 1/2)"と近似できるので、残りの項は"n * (ln(n + 1/2) - ln(m - 1/2))"と近似できる。この近似式自体は精度が良いのでそれ自体は問題ない。
ここに第1の落とし穴がある。単に(ln(n+1/2) - ln(m-1/2))として求めると誤差が大きすぎる！なぜなら、ln(n+1/2),ln(m-1/2)はそれなりの大きさだが、その差分は非常に小さくなりうる。そのため、桁落ちが起こってしまう。さらにそれにnをかける必要があるので絶対誤差も増幅され、1e-9の精度が確保できない。
さて、どのようにして桁落ちを回避するのか。一般に、同程度の数同士の引き算を回避したい。ここで"ln(x/y) = ln(x)-ln(y)"を右から使う。そうすると桁落ちは回避できる。が、第2の落とし穴がある。
なぜなら、"(n+1/2)/(m-1/2)"を計算する時点で(1+ε)となり、εが非常に小さくなりうる。
1があるせいで指数が0になって丸め誤差(この用語は合ってるか不明。以下略)が大きく発生してしまう。
この問題は広く知られており、そのためにC99/C++11には"log1p"という関数が存在する。この関数"log1p(x)"はlog(1+x)を計算するが、xを0に近い値として渡せるために丸め誤差を小さくできる。
しかし単に"log1p( (n+1/2)/(m-1/2)-1)"としたのでは割り算の時点での誤差が回避できていないので、"log1p( (n+1/2)/(m-1/2)-1) = log1p( (n+1/2)/(m-1/2)-(m-1/2)/(m-1/2)) = log1p( (n-m+1)/(m-1/2))"と式変形をする。
ここまでやってやっと適当な精度が確保できる。

コメント

log1pは知らなかった。誤差回避の式変形は面白いと思った。

コード

#include <string>
#include <sstream>
#include <cmath>
using namespace std;
typedef long long ll;   

struct CollectingBonuses {
	ll parse(string s) {
		stringstream ss(s);
		ll x;
		ss >> x;
		return x;
	}
	double expectedBuy(string n_s, string k_s) {
		ll n = parse(n_s), k = parse(k_s);
		//f(k) = 0
		//f(i) = n / (n - i) + f(i+1)
		//f(k-i) = \Sigma_{j=i}^{k-1} (n / (n - j))
		//f(0) = \Sigma_{j=0}^{k-1} (n / (n - j))
		//     = n * \Sigma_{j=0}^{k-1} (1 / (n - j))
		//     = n * \Sigma_{j=n-k+1}^{n} (1 / j)
		//     = n * (H(n) - H(n-k))
		ll m = n - k + 1;
		double r = 0;
		while(m <= 10000000) {
			r += 1. / m;
			if(m == n) return n * r;
			m ++;
		}
		//H(n) - H(m-1) ~= log(n+.5)-log(m-.5) = log((n+.5)/(m-.5))
		//注意: log1p(x) = log(1+x)
		//log((n+.5)/(m-.5)) = log1p((n+.5)/(m-.5)-1) = log1p((n+.5)/(m-.5)-(m-.5)/(m-.5)) = log1p((n-m+1)/(m-.5))
		//・1に近い値より0に近い値のほうが精度は高くなる
		//・log(1+x)はそのために精度が低くなるので、log1p(x) = log(1+x)という関数がある
		//・0に近い値を引き算で求めようとしてもその前に精度が低くなってるから無意味。0に近い値を直接求めなければいけない
		r += log1p((n - m + 1) * 1. / (m - .5));
		return n * r;
	}
};