PageSourceSearch
Opt out API access

https://verifiedby.me/kv/nv-boost/

html verifiedby.me collected 2026-10-02 18:32:19 UTC 16,211 bytes, 438 lines download raw bytes 

1<HTML>
2<HEAD>
3<meta http-equiv="content-type" content="text/html; charset=UTF-8">
4<TITLE> Numerical Verification using Boost C++ Libraries </TITLE>
5</HEAD>
6<BODY>
7
8<div align=right>
9最終更新: 2013/7/4
10</div>
11
12<H1>
13<div align=center>
14boostで精度保証しよう
15</div>
16</H1>
17
18<div align=right>
19柏木 雅英
20</div>
21
22<HR>
23学生の卒論準備のために研究室内向けに作成したページですが、
24一般公開することにしました。
25<HR>
26
27<H2>
281. はじめに
29</H2>
30
31精度保証付き数値計算を行う方法はいろいろあるが、
32特に非線形計算でC++を使った本格的なプログラムを書くことを
33考えたとき、boostに含まれているintervalクラスとublasクラスを
34使うのは一つの有力な方法だと思われる。
35
36<P>
37
38ここに、そのboost.intervalとboost.ublasを使って作った
39比較的単純な非線形方程式の全解探索のプログラムがある。
40これを動かしてみて、boostに慣れてみて欲しい。
41
42<H2>
432. boostとは
44</H2>
45
46boostとは、C++の巨大なテンプレートクラスライブラリ集である。
47C++はオブジェクト指向言語に分類されるが、それ以上にテンプレートによる
48記述力が他の言語の追随を許さぬレベルにある点が特色である。
49
50<P>
51
52テンプレートとは、簡単に言えば、変数の型を未定にしたまま
53クラスや関数を定義する機能である。C言語のマクロ機能を発達させたものという
54見方も出来る。
55<UL>
56<LI>
57関数テンプレートの例:<A HREF="sample/template-func1.cc">template-func1.cc</A>,
58<A HREF="sample/template-func2.cc">template-func2.cc</A>
59<LI>
60クラステンプレートの例:<A HREF="sample/template-class1.cc">template-class1.cc</A>,
61<A HREF="sample/template-class2.cc">template-class2.cc</A>
62</UL>
63
64<P>
65
66C++にはSTLと呼ばれる標準テンプレートクラスライブラリがあるが、
67boostはそれに飽き足らぬ人達が標準化を待たずに勝手に開発した
68テンプレートの集合体である。
69
70<P>
71boostの和書として、
72<UL>
73<LI>
74<A HREF="http://www.amazon.co.jp/dp/4798017043">
75Boost C++ Libraries プログラミング
76</A>
77<LI>
78<A HREF="http://www.amazon.co.jp/dp/4894714191">
79Boost C++をチューンアップする最先端ライブラリ
80</A>
81</UL>
82が出版されている。
83
84<H2>
853. boostのインストール
86</H2>
87
88boostは巨大なクラスライブラリの集合体だが、
89<UL>
90<LI> ヘッダファイルをincludeするだけで使えるクラス
91<LI> ヘッダファイルだけでなく、ライブラリをリンクする必要のあるクラス
92(正規表現、threadなど)
93</UL>
94
95の二種類に分けられる。
96boostの全機能を使おうと思えば、リンクするライブラリをmakeする
97必要がある。その方法は適当に検索するとか本を見るとかして下さい。
98
99<P>
100
101intervalとublasは両方ともヘッダファイルだけで使えるクラスであり、
102そのようなクラスしか使わないならば、boostパッケージのヘッダファイルだけを
103抜き出してどこかに置いておくだけで問題ない。
104boostパッケージを展開して出来る「boost」というディレクトリをどこかに
105置いておき、コンパイル時に「-I(その場所)」とすればよい。
106
107<P>
108
109例えば/usr/local/includeに置いたなら、
110<P>
111cc -I/usr/local/include hoge.cc
112<P>
113カレントディレクトリに置いたなら
114<P>
115cc -I. hoge.cc
116<P>
117
118という感じ。
119
120
121<H3>
1223.1 intervalへのpatch
123</H3>
124
125ublasの方は頻繁に更新されているが、実はintervalの方は
126ほとんど更新されておらず、やや荒削りな部分が見受けられる。
127プログラムを作成しながら、どうしてもここは直して欲しいという
128部分がいくつか見付かったので、本体に直接手を入れることにした。
129
130<P>
131
132以下にその変更点をpatchの形で示す。
133
134<UL>
135<LI> <A HREF="patch/boost-1.50.0.diff">boost 1.50以降用patch</A>
136<FONT SIZE=-1>
137<LI> <A HREF="patch/boost-1.46.1.diff">boost 1.46.1-1.49.0用patch</A>
138<LI> <A HREF="patch/boost-1.44.0.diff">boost 1.44用patch</A>
139<LI> <A HREF="patch/boost-1.43.0.diff">boost 1.43以前用patch</A>
140</FONT>
141</UL>
142
143boostディレクトリの見えている場所で、
144<TT>patch -p0 &lt; boost-1.xx.x.diff</TT>とすればパッチを当てられる。<BR>
145このpatchは、以下の内容を含む。
146
147<UL>
148<LI> 丸めの向きの変更について、
149 <UL>
150 <LI> Windows/Linux/Mac OS X
151 <LI> 32bit/64bit
152 <LI> SSE2あり/無し
153 </UL>
154のいずれの組合せでも正しく高速に動くようにした。
155(特にVisual C++で遅いのがやや改善された。)
156<LI>

156 intervalのコンストラクタに関して、is_convertibleとenable_ifを用いて
157安全化した。具体的には、interval&lt;T&gt;をT1で
158初期化するようなコードを書いたとき、元のコードではは無条件に
159そのコンストラクタが生成されてしまっていたが、
160T1からTに変換可能な場合のみ、そのコンストラクタを生成するようにした。
161複雑に入れ子になったクラスを扱うと元のままでは重大な問題を引き起こす。
162<LI> interval&lt;double&gt;+intのような、
163「interval&lt;T&gt;とTに変換可能な定数の四則演算」を可能にした。
164<LI> -O3など強力に最適化をかけても大丈夫なように、最適化抑制を丁寧にした。
165<LI> intervalの出力に関する軽微な修正。
166<LI> VC++ 64bitでのrint問題の修正。
167<LI> gcc 4.7以降でのto_int問題の修正
168</UL>
169
170柏木の開発しているプログラムは、もはやpatchの当たったintervalを
171前提に書いている。
172
173<BLOCKQUOTE>
174<FONT COLOR="gray" SIZE=-1>
175手抜きするならここに
176<UL>
177<LI>
178<A HREF="headers/boost-1.54.0-headers-patched.tar.gz">
179patchを当てたboostの全header (1.54.0)
180</A>
181<LI> 
182<A HREF="headers/boost-1.54.0-interval-headers.tar.gz">
183patchを当てた<TT>boost/numeric/interval/</TT>以下 (1.54.0)
184</A>
185</UL>
186を置いておくのでこれ使ってもいいよ。
187</FONT>
188</BLOCKQUOTE>
189
190<H2>
1914. 区間演算(interval)
192</H2>
193
194<H3>
1954.1 概要
196</H3>
197
198数値を、[下端, 上端] のような2つの数の組で表される閉区間で表し、
199「結果として有り得る数値の範囲を包含するように」演算を行う方法である。
200区間の上端と下端を表す数の型はテンプレートになっており、
201doubleだけで無く他の型を使うことも出来る。
202
203<H3>
2044.2 使い方
205</H3>
206
207使い方のサンプル
208(<A HREF="sample/test-interval.cc">test-interval.cc</A>)
209を見て欲しい。基本的な使い方はこれで十分と思う。
210詳細は
211<A HREF="http://www.boost.org/doc/libs/">Boost Library Documentation</A>
212の中の
213<A HREF="http://www.boost.org/doc/libs/release/libs/numeric/interval/doc/interval.htm">Intervalの項</A>を参照のこと。
214<A HREF="http://www.kmonos.net/alang/boost/classes/interval.html">日本語の解説</A>も参考になるかも。
215<P>
216<TT>boost::numeric::</TT>とか<TT>std::</TT>とか長くて面倒と
217思うだろうが、<TT>using namespace</TT>を使うのが常套手段だろうか
218(<A HREF="sample/test-interval2.cc">test-interval2.cc</A>)。
219<TT>std</TT>ならともかく<TT>interval</TT>なんて得体の知れないものを
220名前空間に引っ張り出すのは気持ち悪いと思うなら、
221<TT>typedef</TT>を使う手もある
222(<A HREF="sample/test-interval3.cc">test-interval3.cc</A>)。
223<TT>namespace</TT>の別名を使って短縮する方法もある。
224(<A HREF="sample/test-interval4.cc">test-interval4.cc</A>)。
225
226<H3>
2274.3 crlibmと組み合わせる
228</H3>
229
230区間演算を行う上でsinなどの数学関数をどうするかは重大な問題である。
231boost::intervalではsin_down, sin_upなどのsinを下向き丸め、上向き丸めで
232計算するような関数を「自分で用意すれば」sinの区間演算を行うことが出来る。
233そこで、
234<A HREF="http://lipforge.ens-lyon.fr/www/crlibm/">
235CRlibm(Correctly Rounded mathematical library)
236</A>
237と組み合わせて、boostでの数学関数の区間演算を行うことにする。
238<P>
239crlibmのmakeは苦労するかもだが頑張って下さい。
240現状ではたぶんUNIX系OSのみで、windows不可。
241<P>
242boostでcrlibmを使うために、
243<A HREF="sample/use-crlibm.hpp">use-crlibm.hpp</A>
244を作成した。
245<A HREF="sample/test-crlibm.cc">test-crlibm.cc</A>
246のようにuse-crlibm.hppをincludeし、interval型はpolicyを変更して
247crlibmを使うようにする。
248コンパイルするときは
249<PRE>
250g++ -I/usr/local/include -L/usr/local/lib test-crlibm.cc -lcrlibm
251</PRE>
252のように、crlibmのheaderとlibraryのありかをきちんと教えて-lcrlibmを付ける。
253
254<H3>
2554.4 入出力に注意
256</H3>
257
258boost.intervalは、残念なことに10進数と2進数の相互変換に関して

259ほとんどサポートが無い。
260<P>
261数値の入力に関して言えば、
262<A HREF="sample/test-interval5.cc">test-interval5.cc</A>のように
263例えば10進数の0.1を区間演算で扱うには注意が必要である。
264<P>
265数値の出力については、<TT>boost/numeric/interval/io.hpp</TT>を
266includeすれば区間を表示できるが、これは単にC++標準の方法で(最近点への丸めで)
267下端、上端を表示しているだけであり、そのつもりで見る必要がある。
268正しくは、下端は下向き丸めで、上端は上向き丸めで10進数に変換すべきで
269あるが、そうはなっていない。
270
271<H3>
2724.5 丸めの向きの変更
273</H3>
274
275区間演算を行う際の丸めの向きの変更はboost.intervalの内部で
276きちんと行われるが、自作プログラムの一部で丸めの向きを変更したいことも
277あるであろう。その場合は、
278<A HREF="sample/test-rounding.cc">test-rounding.cc</A>のように
279行えばよい。
280ただし、丸めの向きの変更はコンパイラの最適化が大敵で、
281普通にコンパイルしたら正しく動いているように見えたものが <TT>-O3</TT>
282だとうまく動かないようなことがよく起こる。
283<TT>force_rounding</TT>(与えられた引数を<TT>volatile</TT>変数に代入することに
284よって、最適化による計算式の除去を防ぐ)を使うなどして、
285慎重に記述する必要がある。
286
287<H2>
2885. 行列計算(ublas)
289</H2>
290
291<H3>
2925.1 概要
293</H3>
294
295ublasは、行列計算を行うテンプレートライブラリである。
296行列やベクトルの各要素はテンプレートになっており、doubleや
297floatばかりでなくintervalを入れたり、自作の型を入れることも出来る。
298自動微分などの特殊な型を使うことが多い非線形計算の精度保証には
299うってつけのライブラリであると言える。
300
301<P>
302
303なお、blasという名前から高速に動作するようなイメージを持つかも
304知れないが、blasと同等のインターフェースを持つという意味で
305ublasと名付けられており、特に高速なわけではない
306(自分でforループを回して書いたプログラムと同程度)。
307
308<H3>
3095.2 使い方
310</H3>
311
312<UL>
313<LI>
314使い方のサンプル
315(<A HREF="sample/test-ublas.cc">test-ublas.cc</A>)
316<LI>
317LU分解を利用した逆行列計算のサンプル
318(<A HREF="sample/test-invert.cc">test-invert.cc</A>,
319<A HREF="sample/matrix-inversion.hpp">matrix-inversion.hpp</A>)
320<LI>
321要素を区間にしたサンプル
322(<A HREF="sample/test-ublas-interval.cc">test-ublas-interval.cc</A>)
323</UL>
324
325詳細は
326<A HREF="http://www.boost.org/doc/libs/">Boost Library Documentation</A>
327の中の
328<A HREF="http://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm">ublasの項</A>を参照のこと。
329<A HREF="http://www.kmonos.net/alang/boost/classes/ublas.html">日本語の解説</A>も参考になるかも。
330
331<A HREF="http://www.page.sannet.ne.jp/d_takahashi/boost/ublas/index.html">
332boost::numeric::ublas 線形代数ライブラリの使い方
333</A>
334も大変参考になる。
335
336
337<P>
338なお、STLの一部として<TT>std::vector</TT>が存在するため、
339<TT>using namespace std;</TT>している場合などに
340<TT>boost::numeric::ublas::vector</TT>と衝突してエラーになりやすい
341(<A HREF="sample/test-ublas2.cc">test-ublas2.cc</A>)。
342<A HREF="sample/test-ublas3.cc">test-ublas3.cc</A>のように
343ublasのvectorは明示的に指定することにしてusing namespaceに頼らないのが
344安全かも知れない。
345
346<H3>
3475.3 マクロ<TT>NDEBUG</TT>について
348</H3>
349ublasではExpression Templateと呼ばれる技法を用い、
350余分な一時変数を作成しないことで高速化を図っている。
351但し、この機能はコンパイル時に <TT>-DNDEBUG</TT>を付けて
352<TT>NDEBUG</TT>マクロを定義した場合のみ有効になる。
353<TT>NDEBUG</TT>はC++においてrelease版(debug版で無い)を
354makeする場合に付加するマクロとして一般的であり、実際Visual C++で
355release版としてbuildすると自動的にこのマクロが付加される。
356<P>
357なお、ublasでは、<TT>NDEBUG</TT>マクロが定義されていない場合は
358配列の添字溢れのチェックが行われるが、<TT>NDEBUG</TT>を定義すると
359チェックが省かれて高速化されるようになっている。
360
361<H2>
3626. 全解探索プログラムを動かす
363</H2>
364
365大きめのプログラムとして、非線形方程式の全解探索を行う

366プログラムのサンプルを示す。高速化のための凝った工夫を取り除いた
367シンプルなものである。
368
369<UL>
370<LI>
371<A HREF="sample/allsol-simple.hpp">allsol-simple.hpp</A> : 全解探索プログラム本体
372<LI>
373<A HREF="sample/test-allsol-simple.cc">test-allsol-simple.cc</A> : 全解探索を行うsample (main)
374<LI>
375<A HREF="sample/autodif.hpp">autodif.hpp</A> : 自動微分クラス
376<LI>
377<A HREF="sample/ivector-uty.hpp">ivector-uty.hpp</A> : 区間行列・ベクトルに関するいくつかの関数
378<LI>
379<A HREF="sample/matrix-inversion.hpp">matrix-inversion.hpp</A> : 逆行列
380</UL>
381
382必要なのはこれらのファイル。
383<A HREF="sample/test-allsol-simple.cc">test-allsol-simple.cc</A>
384をコンパイル、実行すれば、
385その中で定義された非線形関数の全解探索が行われる。
386
387
388<P>
389
390<A HREF="sample/test-allsol-simple.cc">test-allsol-simple.cc</A>
391での解きたい関数の定義の仕方が分かりにくいかも知れない。
392解きたい関数ごとにクラスを定義し、そのクラスでは
393<TT>operator()</TT>をオーバーロードする。
394そのクラスのインスタンス<TT>f</TT>に対して<TT>f(x)</TT>と書くと
395<TT>operator()</TT>が引数<TT>x</TT>で呼び出されるので、
396そこに関数の定義を書き、インスタンス<TT>f</TT>を関数とみなす。
397これは関数オブジェクトと呼ばれる技法で、関数ポインタを使うよりも
398一般的にオーバーヘッドが少ない。
399
400<H2>
4017. 数値型のプログラミング
402</H2>
403
404ここまでは、既存のintervalとublasの使い方を説明してきたが、
405精度保証付き数値計算では、自動微分などそれ以外の「数値っぽい振る舞いをする型」
406を自作する必要がある。
407特に自作型を入れ子にする場合など、注意すべき点も多いので、
408その方法を <A HREF="numbertype/index.html">数値型のプログラミング</A>
409にまとめた。
410
411<H2>
4128. Affine Arithmetic
413</H2>
414
415<A HREF="affine/index.html">Affine Arithmetic</A>
416にまとめた。
417
418<H2>
4199. おわりに
420</H2>
421
422まだ書いてない。自動微分とかpsaとかについて書く予定。
423アプリケーションとして、全解探索とかODEとか数値積分とかも。
424<P>
425ublasで内部型がdoubleの場合にlapack/blasを使う方法、
426intervalの内部型としてmpfr型を使い任意精度区間演算を行う方法、
427なども書いていきたい。
428
429<H2>
430変更履歴
431</H2>
432
433<A HREF="history.html">変更履歴</A>参照。
434
435<HR>
436
437</BODY>
438</HTML>

Line numbers count LF bytes from the start of the resource, as the search results do. Vendor segments are library code the classifier recognised; they are stored but not indexed. Bytes are shown as Latin1 characters, one per byte.