2021/11/30(火)NVR2021

11月27,28日に、NVR2021という研究集会を開催しました。

今回で5回めで、コロナ禍で昨年に引き続き残念ながらオンラインです。今年は、大石CRESTが最終年度で成果報告会を開催してほしいという要請があり、時期も近くて参加メンバーも多くが被るので、共同開催することにしました。CREST報告会が土曜日、NVRが日曜日でした。いろいろ重なって大変でしたが、無事終わってホッとしています。

使用したスライドを張っておきます。1つ目はdouble-doubleに関する話題、2つ目はKrawczyk法の候補者集合に関する話題です。2つ目は軽い話題で発表時間も短め(20分)です。

2021/11/12(金)kv-0.4.53

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

今回は2つのbug-fixです。

1つ目は、非常に特殊な条件ですが、xがdouble-double数で、その上端が(DBL_MAX+正数)というoverflow寸前の巨大数だった場合に、sqrt(x)の計算結果の幅が広くなってしまう(double-doubleとしてふさわしい精度が出ずにdouble程度の区間幅になってしまう)というバグの修正です。実際、
#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

int main(){
        kv::interval<kv::dd> x;
        std::cout.precision(34);

        x = std::numeric_limits<kv::dd>::max();

        std::cout << x << "\n";
        std::cout << sqrt(x) << "\n";
}
を実行すると、
[1.340780792994259561100831764080293e+154,1.340780792994259672743259815885529e+154]
と計算されてしまっていました。修正後は、
[1.340780792994259672743259815885478e+154,1.340780792994259672743259815885529e+154]
ときちんと計算されるようになりました。中田真秀先生の、
のtweetで気づくことができました(kvではなくQDライブラリへの言及ではありますが)。中田先生ありがとうございました。

2つめは、自動step幅調節機能付き定積分を行う、defint_autostepで、定積分
\int_a^b f(x) dx
を計算するとき、終端値bが区間(例えば\pi)だと無限ループ(あるいは非常に長時間のループ)に陥ってしまうことがあるというバグです。この関数は、(尺取り虫のように自動的に計算精度を見ながら)積分区間を
a=t_0, t_1, t_2, \ldots, t_{n-1}, t_n=b
のように分割して積分を行っていますが、そのプロセスは2passになっていて、t_iでTaylor展開の係数を計算して適切なstep sizeを推定し、[t_i,t_{i+1}]での精度保証付きの積分値を計算し、その区間幅を見て広すぎるか狭すぎるかしたら再度step sizeを修正してもう一度精度保証付きの積分値を計算する、ということをしています。ここで、[t_i,t_{i+1}]での積分値を計算するとき、t_it_{i+1}のどちらかが(幅が広めの)区間だった場合は積分値の計算結果の区間幅が極端に大きくなり、精度を出そうとしてt_{i+1}を引き戻してstep sizeを小さくしようとします(そんなことをしても精度は上がらないのですが)。従って、bが区間だった場合には[t_i,t_{i+1}=b]と最後の一歩を計算しようとすると、精度不足でt_{i+1}が引き戻されてしまい、いつまで経ってもbに到達しない、という現象が発生していました。最後の一歩のときはstep sizeの見直しをしないように修正しました。ずっと存在していたバグですが、0.4.51以前は最後の一歩の区間幅の再計算が広めになるバグがあってたまたま発生しづらくなっていて、0.4.52で直ったことにより顕在化したようです。

この現象は、浅井 大晴さん、多田 秀介さん、宮内 洋明さんの3名の、Bessel関数J_0(x)の計算が終わらないという指摘で気付きました。Bessel関数は、\displaystyle J_n(x) = \frac{1}{\pi} \int_0^{\pi} \cos(x \sin t - nt) dt で計算しており、終端値が区間である積分をしていたので、この現象が発生しました。みなさま、ありがとうございました。

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となるようなnF(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)=13562027F(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まで証明できました。

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

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を高精度化したのは、ある問題を解くのに必要になったためです。忘れないうちにその話も記事にするつもりです。

2021/03/20(土)kv-0.4.51

ものすごく久しぶりに、kvライブラリを0.4.51にアップデートしました。

変更点は以下の通りです。
  • 区間演算を行うinterval.hppにfloorとceilを追加。
  • 区間演算と領域の再帰的分割による最適化を行うoptimize.hppが場合によって正しく終了しなかったバグの修正。
  • 高階微分を計算するhighderiv.hppが定数関数を渡されたとき動作しなかったバグの修正。
  • 昨年11月のNVR2020の発表で作成した、複合Newton-Cotes公式を用いた数値積分のための関数群を追加。
大した変更や追加はないのですが、コロナ対応での多忙で時間が取れず、3月になってようやく少し時間が出来たので思い出しながらファイルを整理し、公開できる状態までもっていきました。数値積分のところのドキュメントに手間取りました。

highderiv.hppとoptimize.hppはともに数値積分の誤差項の計算に使われており、昨年11月の発表の際にいろいろ使ってみて不具合が発覚しています。やはり使われて鍛えられないとライブラリの完成度は上がらないので、ご利用下さると大変嬉しいです。
OK キャンセル 確認 その他