2021/11/30(火)NVR2021
今回で5回めで、コロナ禍で昨年に引き続き残念ながらオンラインです。今年は、大石CRESTが最終年度で成果報告会を開催してほしいという要請があり、時期も近くて参加メンバーも多くが被るので、共同開催することにしました。CREST報告会が土曜日、NVRが日曜日でした。いろいろ重なって大変でしたが、無事終わってホッとしています。
使用したスライドを張っておきます。1つ目はdouble-doubleに関する話題、2つ目はKrawczyk法の候補者集合に関する話題です。2つ目は軽い話題で発表時間も短め(20分)です。
2021/11/12(金)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_iかt_{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となるような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 |
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まで証明できました。数学としては何の意味もない問題でしょうけど、こういうのは楽しいですね。
2021/08/05(木)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関数の計算アルゴリズムをより高精度なものに変更。
- その他細かい修正。
また、digammaを高精度化したのは、ある問題を解くのに必要になったためです。忘れないうちにその話も記事にするつもりです。
2021/03/20(土)kv-0.4.51
変更点は以下の通りです。
- 区間演算を行うinterval.hppにfloorとceilを追加。
- 区間演算と領域の再帰的分割による最適化を行うoptimize.hppが場合によって正しく終了しなかったバグの修正。
- 高階微分を計算するhighderiv.hppが定数関数を渡されたとき動作しなかったバグの修正。
- 昨年11月のNVR2020の発表で作成した、複合Newton-Cotes公式を用いた数値積分のための関数群を追加。
highderiv.hppとoptimize.hppはともに数値積分の誤差項の計算に使われており、昨年11月の発表の際にいろいろ使ってみて不具合が発覚しています。やはり使われて鍛えられないとライブラリの完成度は上がらないので、ご利用下さると大変嬉しいです。