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++;
	}
}
計算結果は、
xF(x)
24
311
431
583
6227
7616
81674
94550
1012367
1133617
1291380
13248397
14675214
151835421
164989191
1713562027
1836865412
19100210581
20272400600
21740461601
222012783315
235471312310
2414872568832
2540427833562
26109894245237
27298723533129
28812014582525
292207275035877
306000125006292
となりました。果たしてこれは正しいでしょうか？

区間演算してみる

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時間ほどかかりました(汗)。

計算結果をまとめておきます。
xF(x)
24
311
431
583
6227
7616
81674
94550
1012367
1133617
1291380
13248397
14675214
151835421
164989191
1713562027
1836865412
19100210581
20272400600
21740461601
222012783315
235471312310
2414872568831
2540427833596
26109894245429
27298723530401
28812014744422
292207284924203
306000022499693
こうして見ると、精度保証しない倍精度での計算結果は、F(23)までは(偶然)正しく、F(24)以降は間違っていたことが分かりました。下の赤字のところが間違っています。
xF(x)
235471312310
2414872568832
2540427833562
26109894245237
27298723533129
28812014582525
292207275035877
306000125006292

もうちょっと数学っぽく

有名な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まで証明できました。

数学としては何の意味もない問題でしょうけど、こういうのは楽しいですね。
OK キャンセル 確認 その他