kashiの日記

2021/08/05(木)調和級数の部分和で遊んでみた

はじめに

次のようなツイートが話題になったことがありました。

調和級数

H_n = \sum_{k=1}^n \frac{1}{k} = 1 + \frac{1}{2} + \frac{1}{3} + \cdots + \frac{1}{n}

が、n \to \inftyのとき無限に発散することはよく知られています。この発散は非常に遅く、30を超える最小のnはかなり大きな数になります。これを計算で求めようと思うと、計算時間の問題もありますが、何よりも丸め誤差の蓄積の影響が馬鹿にならず、正しいnを求めるのは案外大変です。

以下、xに対してH_{n-1}<xだがH_n \ge xとなるようなnをF(x)と書くことにします。つまり、F(30)を計算したい。

精度保証せずに普通に計算

次の簡単なプログラムで、自然数なxに対してF(x)を計算できます。

#include <iostream>
#include <cmath>

int main()
{
	double s, tmp;
	long long int i;
	std::cout.precision(17);

	s = 0;
	i = 1;
	while (true) {
		tmp = s + 1.0 / i;
		if (floor(tmp) != floor(s)) {
			std::cout << i-1 << ":" << s << "\n";
			std::cout << i << ":" << tmp << "\n";
		}
		s = tmp;
		i++;
	}
}

計算結果は、

x	F(x)
2	4
3	11
4	31
5	83
6	227
7	616
8	1674
9	4550
10	12367
11	33617
12	91380
13	248397
14	675214
15	1835421
16	4989191
17	13562027
18	36865412
19	100210581
20	272400600
21	740461601
22	2012783315
23	5471312310
24	14872568832
25	40427833562
26	109894245237
27	298723533129
28	812014582525
29	2207275035877
30	6000125006292

となりました。果たしてこれは正しいでしょうか？

区間演算してみる

kvライブラリを使って、区間演算でこれが正しいかチェックしてみます。部分和を区間演算して、xを「またいで」いないかどうかチェックします。

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

typedef kv::interval<double> itv;
// typedef kv::interval< kv::dd > itv;

int main()
{
	itv s, tmp;
	long long int i;
	std::cout.precision(17);

	s = 0;
	i = 1;
	while (true) {
		tmp = s + 1.0 / itv(i);
		if (floor(tmp.lower()) != floor(tmp.upper())) {
			std::cout << "warning\n";
			std::cout << i-1 << ":" << s << "\n";
			std::cout << i << ":" << tmp << "\n";
			break;
		} else if (floor(tmp.lower()) != floor(s.lower())) {
			std::cout << i-1 << ":" << s << "\n";
			std::cout << i << ":" << tmp << "\n";
		}
		s = tmp;
		i++;
	}
}

これを実行すると、x=16のときは

\begin{eqnarray*} && H_{4989190} \in [15.999999890595899,15.999999899456953] \\ && H_{4989191} \in [16.000000091029196,16.000000099890251] \end{eqnarray*}

となって間違いなくF(16)=4989191と証明されました。しかし、x=17のときは、

\begin{eqnarray*} && H_{13562026} \in [16.999999921467317,16.999999960785197] \\ && H_{13562027} \in [16.999999995202607,17.00000003452049] \end{eqnarray*}

となってしまい、F(17)=13562027とは断定できないことが分かりました。

なお、このプログラムは単に和を左から足していますが、すべて正のたくさんの数を足すときは、小さい順に足した方が丸め誤差が小さくなり、区間演算の場合は区間の幅が小さくなることが知られています。大体のnが分かっているならば、逆順に足すことによってもう少し頑張れて、F(17)=13562027、F(18)=36865412まで証明することができます。

double-double区間演算で頑張る

倍精度ではF(16)まで、逆順に足す工夫をしてもF(18)までしか証明できませんでした。kvライブラリはdouble型を2つ束ねた疑似4倍精度で区間演算ができるので、それを使って頑張ってみました。上に挙げたプログラムの最初の方のtypedef行を変更するだけです。これで計算すると、x=30で

\begin{eqnarray*} && H_{6000022499692} \in [29.999999999999855,29.999999999999856] \\ && H_{6000022499693} \in [30.000000000000021,30.000000000000022] \end{eqnarray*}

となり、F(30)=6000022499693と判明しました。最初のtweetを見て計算してみた人は何人かいたようですが、やはり丸め誤差まみれになってしまってぴったりの値を計算するのは難しかったようです。ただ、この計算にはOpenMPを使って並列化しても55時間ほどかかりました(汗)。

計算結果をまとめておきます。

x	F(x)
2	4
3	11
4	31
5	83
6	227
7	616
8	1674
9	4550
10	12367
11	33617
12	91380
13	248397
14	675214
15	1835421
16	4989191
17	13562027
18	36865412
19	100210581
20	272400600
21	740461601
22	2012783315
23	5471312310
24	14872568831
25	40427833596
26	109894245429
27	298723530401
28	812014744422
29	2207284924203
30	6000022499693

こうして見ると、精度保証しない倍精度での計算結果は、F(23)までは(偶然)正しく、F(24)以降は間違っていたことが分かりました。下の赤字のところが間違っています。

x	F(x)
23	5471312310
24	14872568832
25	40427833562
26	109894245237
27	298723533129
28	812014582525
29	2207275035877
30	6000125006292

もうちょっと数学っぽく

有名なOEISを調べてみると、この数列はA002387に見つかりました。そこにはF(28)=812014744422まで載っていました。また、そのページからリンクされたここにはF(2303)まで(!)載っていますが、これが証明されたものかどうかは分かりません。自然数n>1に対して

F(n) = \lfloor \exp(n - \gamma) + \frac{1}{2} \rfloor

(\lfloor\rfloorはガウス記号、\gammaはEuler's gamma)という予想があるようなので、これで計算したものなのかも。

また、ガンマ関数の対数微分であるディガンマ関数

\mathrm{digamma}(x) = \frac{d}{dx}\log\Gamma(x) = \frac{\Gamma'(x)}{\Gamma(x)}

は、漸化式

\mathrm{digamma}(x+1) = \mathrm{digamma}(x) + \frac{1}{x}

を満たすことが知られており、これから直ちに、

H_n = \sum_{k=1}^n \frac{1}{k} = \mathrm{digamma}(n+1) - \mathrm{digamma}(1)

が従います。よって、与えられたxに対して、

n = \lfloor \exp(x - \gamma) + \frac{1}{2} \rfloor

で予測して、H_{n-1}とH_nを精度保証付きのディガンマ関数で計算し、

H_{n-1} < x \le H_n

が満たされていればF(x)=nが証明されたことになります。

幸いなことにkvライブラリには精度保証付きのディガンマ関数があるので、これを実装してみたところ、まったく精度が出ず話になりませんでした。原因はちょっと考えたらすぐ分かって、そもそもkvのディガンマ関数は

\mathrm{digamma}(x+1) = \mathrm{digamma}(x) + \frac{1}{x}

の漸化式を用いて1 \le x < 2の範囲に引き戻して計算する実装になっていて、何のことはない、ただ調和級数の部分和を計算するだけになっていたのでした。そこで、kv-0.4.52ではもうちょっと真面目にディガンマ関数を実装してみたという話でした。詳細はガンマ関数の精度保証付き計算メモの3章に書きましたが、簡単に言えば、xを引き戻した上で

\mathrm{digamma}(x) = \int_0^\infty \left( \frac{e^{-t}}{t} - \frac{e^{-xt}}{1-e^{-t}} \right) dt

を精度保証付き数値積分で計算していたのを、引き戻しを行わずに、

\mathrm{digamma}(x) = \log x - \frac{1}{2x} - 2 \int_0^\infty \frac{t}{(t^2+x^2)(e^{2 \pi t}-1)}dt

を使って計算するようにしました。

#include <iostream>
#include <kv/interval.hpp>
#ifdef DD
#include <kv/dd.hpp>
#include <kv/rdd.hpp>
#endif
#include <kv/gamma.hpp>

#ifdef DD
typedef kv::interval<kv::dd> itv;
#define MAX 66
#else
typedef kv::interval<double> itv;
#define MAX 32
#endif

int main() {
	#ifdef DD
	std::cout.precision(34);
	#else
	std::cout.precision(17);
	#endif

	itv gamma, x, y;
	int i;

	gamma = - kv::digamma(itv(1.));

	for (i=2; i<=MAX; i++) {
		std::cout << i << "\n";

		x = i;
		y = floor(mid(exp(x - gamma) + 0.5));
		std::cout << y.lower() << "\n";
		std::cout <<  kv::digamma(y) + gamma << "\n";;
		std::cout <<  kv::digamma(y + 1) + gamma << "\n";;
	}
}

新しいディガンマ関数を用いて計算すると、倍精度で

\begin{eqnarray*} && H_{16309752131261} \in [30.99999999999995,30.999999999999997] \\ && H_{16309752131262} \in [31.000000000000014,31.000000000000061] \end{eqnarray*}

となってF(31)=16309752131262が、double-doubleを使うと、

\begin{eqnarray*} && H_{9516116398696941841501814186} \in [64.99999999999999999999999999997317,64.99999999999999999999999999999567] \\ && H_{9516116398696941841501814187} \in [65.00000000000000000000000000007849,65.00000000000000000000000000010098] \end{eqnarray*}

となってF(65)=9516116398696941841501814187まで証明できました。

数学としては何の意味もない問題でしょうけど、こういうのは楽しいですね。

コメント（0件）

2021/08/05(木)kv-0.4.52

精度保証

kvライブラリを0.4.52にアップデートしました。

変更点は以下の通りです。

test-rounding.ccで、非正規化数を強制的に0にする"flush to zero"モードになっていないかどうかをチェックするようにした。
iccでコンパイルできるように修正。
optimize.hppで、「端にない」区間で傾きが単調ならばそこに最小値は存在しないというチェックを入れる改良を行った。
defint-newtoncotes.hpp, double-newtoncotes.hppで、最適な分割数の推定がうまくできていなかったバグを修正。
defint.hpp, defint-singular.hppで、急減少関数などの特定の被積分関数に対してstep sizeがうまく計算できていなかった問題を修正。
psa.hppの中のinv関数を、係数のオーバーフローが発生しづらいように改良。
gamma.hppの中のdigamma関数の計算アルゴリズムをより高精度なものに変更。
その他細かい修正。

iccはデフォルトで非正規化数を強制的に0にする"flush to zero"になっているので、そのままだと区間演算の結果が異常になってしまいます。例えば、上向き丸めモードで2^{-1022}/2を計算すると、0になってしまいます。すなわち、IEEE754に従っていません。オプション-no-ftz (windowsなら/Qftz-)を付けるとこの問題は回避できるようです。iccを使うときは、test-rounding.ccをコンパイル、実行してチェックしてみて下さい。

また、digammaを高精度化したのは、ある問題を解くのに必要になったためです。忘れないうちにその話も記事にするつもりです。

コメント（0件）

検索条件