C++で学ぶ精度保証付き数値計算法の初歩

東洋大学 情報連携学部 関根 晃太

1数値計算とは?

数値計算とは?
定義:
コンピュータを用いて，手で解くことが困難な数学の問題を近似的に解く方法

2

数値計算とは?
数値計算 :
コンピュータを用いて，手で解くことが困難な数学の問題を近似的に解く方法
知りたい情報モデル化(数学で記述する)
数学モデル化
コンピュータで解ける式に変形
数値モデル化プログラムと実行
数値結果
3

精度保証付き数値計算とは?

あなたのコンピュータは正しい結果を返していますか?
4精度保証付き数値計算とは?
精度保証付き数値計算 :与えられた問題(数学モデル)の解の存在範囲もしくは一意存在の範囲を丸め誤差の厳密評価を含めて特定する算法
知りたい情報モデル化(数学で記述する)
数学モデル化
コンピュータで解ける式に変形
数値モデル化プログラムと実行
数値結果
5

精度保証付き数値計算とは?
精度保証付き数値計算 :与えられた問題(数学モデル)の解の存在範囲もしくは一意存在の範囲を丸め誤差の厳密評価を含めて特定する算法
知りたい情報近似的に計算数学モデル化
数値モデル化誤差が入る!!

数値結果
6

精度保証付き数値計算とは?
・誤差についてプログラムのどこで誤差が発生するか?

・区間演算について四則演算の誤差はどのように把握するか?
・ライブラリについてもっと簡単に使えないか?

7

精度保証付き数値計算とは?
・誤差についてプログラムのどこで誤差が発生するか?

・区間演算について四則演算の誤差はどのように把握するか?
・ライブラリについてもっと簡単に使えないか?

8プログラムを少し書きます!
どこに誤差が発生するかわかりますか??
#include <iostream>
int main(void){
double a = 0.1;
double b = 0.01;
double c = a + b;
std::cout << c << std::endl;}

return 0;9精度保証付き数値計算とは?
コンピュータで扱える数値の一つ:
IEEE 754 Standardで定められている浮動小数点数
倍精度浮動小数点数の場合…
S, ei, mjを全て2進数決まった数しか表現できない…

S e10 e9 ・・・ e2 e1 e0 m51 m50 m49 m48 ・・・ m3 m2 m1 m0
符号1bit

指数部
11bit

値: (-1)S×(1+m)×2e-1023

仮数部
52bit

整数:e=e0×20 +・・・+e10×210
小数:m=m0×20-52 +・・・+m52×252-52
10精度保証付き数値計算とは?
Q. 10進数の0.1と0.01は倍精度浮動小数点数で表せますか?

11精度保証付き数値計算とは?
Q. 10進数の0.1と0.01は倍精度浮動小数点数で表せますか?
A. いいえ…
近似されるため，誤差発生!!

浮動小数点数

実数軸

0.1
特に，入力による誤差は非常に把握しにくい!!
(上に丸められたか，下に丸められたか判断できない…)
12

プログラムを少し書きます!
どこに誤差が発生するかわかりますか??
#include <iostream>int main(void){
double a = 0.1;
double b = 0.01;
double c = a + b;
std::cout << c << std::endl;}

誤差が発生するのでそのままではダメ…

return 0;13精度保証付き数値計算とは?
Q. 10進数の1と3は倍精度浮動小数点数で表せます.
その上で，1/3の結果は浮動小数点数で表せますか?

14精度保証付き数値計算とは?
Q. 10進数の1と3は倍精度浮動小数点数で表せます.
その上で，1/3の結果は浮動小数点数で表せますか?
A. いいえ…
但し，IEEE754 Standardで定められている丸めモードの指定をすることでどこの浮動小数点数にいったか把握可能!!
最近点丸め
下向き 上向き

浮動小数点数
実数軸

1/3

15

プログラムを少し書きます!
どこに誤差が発生するかわかりますか??
#include <iostream>
int main(void){
double a = 0.1;
double b = 0.01;
double c = a + b;
std::cout << c << std::endl;}

演算でも誤差が発生するのでそのままではダメ…

return 0;

16

プログラムを少し書きます!
どこに誤差が発生するかわかりますか??
#include <iostream>
int main(void){
double a = 0.1;
double b = 0.01;
double c = a + b;
std::cout << c << std::endl;}

return 0;

出力は??

17

プログラムを少し書きます!
どこに誤差が発生するかわかりますか??
#include <iostream>
int main(void){
double a = 0.1;
double b = 0.01;
double c = a + b;
std::cout << c << std::endl;}

return 0;

出力は??
そもそも決められておらず，何しているかわからない…
(浮動小数点演算が使われていたら，もちろんNG!)
18

プログラムを少し書きます!
どこに誤差が発生するかわかりますか??
#include <iostream>
int main(void){
double a = 0.1;
double b = 0.01;
double c = a + b;
std::cout << c << std::endl;}

return 0;

出力は??
そもそも決められておらず，何しているかわからない…
(浮動小数点演算が使われていたら，もちろんNG!)
19

精度保証付き数値計算とは?
・誤差についてプログラムのどこで誤差が発生するか?

・区間演算について四則演算の誤差はどのように把握するか?
・ライブラリについてもっと簡単に使えないか?

20

精度保証付き数値計算とは?
精度保証付き数値計算では…のように区間で表現し，常に真の値を含むようにする!!

浮動小数点数
実数軸0.121

精度保証付き数値計算とは?
もしかすると，精度保証付き数値計算を聞いたことがある人は…
丸め変更で簡単に何とかなるんじゃないの!?って思うかもしれません

22精度保証付き数値計算とは?
Q. 大抵のCPUがIEEE754 Standardに準拠しています.
では，丸めモードが標準で変更可能なプログラミング言語は??

23精度保証付き数値計算とは?
Q. 大抵のCPUがIEEE754 Standardに準拠しています.
では，丸めモードが標準で変更可能なプログラミング言語は??
A. C99(C/C++)
Fortran 2003
MATLAB
Octave・Pythonなどは標準では備えていない…
・C#は丸めモードの変更を非推奨

24精度保証付き数値計算とは?
Q. 大抵のCPUがIEEE754 Standardに準拠しています.
では，丸めモードが標準で変更可能なプログラミング言語は??
C99を使った丸めモードの変更

A. C99(C/C++)
Fortran 2003
MATLAB
Octave#include <fenv.h>
int main(void) {// 下向き丸めに変更
fesetround(FE_DOWNWARD );// 上向き丸めに変更
fesetround(FE_UPWARD );// 最近点丸めに変更
fesetround(FE_TONEAREST );}25

精度保証付き数値計算とは?
#include <iostream>#include <fenv.h>
int main(void) {// 下向き丸めに変更
fesetround(FE_DOWNWARD );double a_down = 0.1;}// 上向き丸めに変更
fesetround(FE_UPWARD );double a_up = 0.1;26

精度保証付き数値計算とは?
#include <iostream>#include <fenv.h>
int main(void) {// 下向き丸めに変更
fesetround(FE_DOWNWARD );double a_down = 0.1;}

丸め変更は四則演算(+, -, *, /)と平方根(sqrt)のみ対応!

// 上向き丸めに変更
fesetround(FE_UPWARD );double a_up = 0.1;

27

精度保証付き数値計算とは?
Q. なぜ下記のプログラムはダメなのでしょうか?
C99を使った1/3の包含#include <fenv.h>
int main(void) {double a, b;
fesetround(FE_DOWNWARD );a = 1.0/ 3.0;
fesetround(FE_UPWARD );b = 1.0/ 3.0;}28

精度保証付き数値計算とは?
Q. なぜ下記のプログラムはダメなのでしょうか?
C99を使った1/3の包含

最適化後の例#include <fenv.h>#include <fenv.h>

int main(void) {double a, b;
fesetround(FE_DOWNWARD );a = 1.0/ 3.0;
fesetround(FE_UPWARD );b = 1.0/ 3.0;}

int main(void) {double a, b;
fesetround(FE_DOWNWARD );a = 1.0/ 3.0;
fesetround(FE_UPWARD );b = a;}

コンパイラのオプションでコードが最適化されます…
(「-O0」で実行すると遅い…)

29

精度保証付き数値計算とは?
解決策(?)
C99にはFENV_ACCESSという丸めモードの変更に伴う最適化の抑制を行うコマンドが制定されている!!

30

精度保証付き数値計算とは?
解決策(?)
C99にはFENV_ACCESSという丸めモードの変更に伴う最適化の抑制を行うコマンドが制定されている!!
が，しかし，実装されているコンパイラは少ない…
(例えばgcc 9.1にはない…)

31

精度保証付き数値計算とは?
解決策(?)
C99にはFENV_ACCESSという丸めモードの変更に伴う最適化の抑制を行うコマンドが制定されている!!
が，しかし，実装されているコンパイラは少ない…
(例えばgcc 9.1にはない…)
現状は，volatile修飾子で部分的に最適化を抑制…
(が，しかし，最近制定されたc++20では，volatileは本来の目的以外で使用することを非推奨
将来はダメにするといわれ，現在はガクブル状態…)
32

精度保証付き数値計算とは?
今までのまとめ・演算ごとにも誤差が発生する・入出力でも誤差が発生する・精度保証では真の値を包含した区間で表現
・最適化によって演算順序が変わると精度保証できない

33

精度保証付き数値計算とは?
今までのまとめ・演算ごとにも誤差が発生する・入出力でも誤差が発生する・精度保証では真の値を包含した区間で表現
・最適化によって演算順序が変わると精度保証できない
0から入出力, 四則演算をつくり，そこから自分の解きたい問題まで…ってやっていたら非常に大変, やだなぁ…
34

精度保証付き数値計算とは?
今までのまとめ・演算ごとにも誤差が発生する・入出力でも誤差が発生する・精度保証では真の値を包含した区間で表現
・最適化によって演算順序が変わると精度保証できない
自分で区間演算を実装するのはやめて，既存のライブラリを使いましょう!
35

2.有名な精度保証付き数値計算ライブラリ

INTLAB:
言語:
MATLAB
特徴:
HP参照
数値線形代数に対する精度保証付き数値計算法の専門家であるS.M. Rump 氏が作成したライブラリ．MATLAB に実装されているBLAS とLAPACK を利用した高速・高精度な行列，数値線形代数に関わる精度保証 付き数値計算が強力．
S.M. Rump, INTLAB - INTerval LABoratory, In Tibor Csendes, editor, Developments
in Reliable Computing, pp.77104. Kluwer Academic Publishers, Dordrecht, 1999.
http://www.ti3.tuhh.de/rump/intlab/
36

2.有名な精度保証付き数値計算ライブラリkv Library:
言語:
C++
特徴:
HP参照
常微分方程式の解に対する精度保証付き数値計算法が専門家である柏木雅英氏が作成したライブラリ．常微分方程式の精度保証付き数値計算法が実装されている．また，Double-Double演算やMPFR を用いた高精度なスカラー計算に対する精度保証付き数値計算法が実装されている．
M. Kashiwagi: kv - C++ による精度保証付き数値計算ライブラリ,
http://veriedby.me/kv/37

kv ライブラリを使った精度保証プログラム
ここからは実際に，・インストール手順・プログラミングなどを実演します!
プログラミングのメモのみ残しておきます．

38

kv ライブラリを使った精度保証プログラム#include <iostream>
int main(void){
double a = 0.1;
double b = 0.01;
double c = a + b;
std::cout << c << std::endl;}

return 0;#include <iostream>#include <kv/interval.hpp>#include <kv/rdouble.hpp>
typedef kv::interval< double > idouble;
int main(void){
idouble a = idouble("0.1");
idouble b = idouble("0.01");
idouble c = a + b;
std::cout << c << std::endl;}

return 0;

39

VCPを使った近似行列プログラム(1)#include <iostream>#include <vcp/matrix.hpp>#include <vcp/matrix_assist.hpp>
typedef vcp::matrix< double, vcp::mats< double > > vmat;
int main(void){
int n = 5;
vmat A, B, C;
A.zeros(n, n);B.rand(n, n);
C = A + B;
std::cout << C << std::endl;}

return 0;40BLASとVCPを使った近似行列プログラム#include <iostream>#include <vcp/matrix.hpp>#include <vcp/matrix_assist.hpp>#include <vcp/pdblas.hpp>#include <vcp/vcp_timer.hpp>//typedef vcp::matrix< double, vcp::mats< double > > vmat;
typedef vcp::matrix< double, vcp::pdblas > vmat;
int main(void){
int n = 5000;
vmat A, B, C;
A.rand(n, n);B.rand(n, n);
vcp::time.tic();
C = A * B;
vcp::time.toc();
std::cout << C(0, 0) << std::endl;}

return 0;41

BLASとVCPを使った近似連立一次プログラム#include <iostream>#include <vcp/matrix.hpp>#include <vcp/matrix_assist.hpp>#include <vcp/pdblas.hpp>#include <vcp/vcp_timer.hpp>//typedef vcp::matrix< double, vcp::mats< double > > vmat;
typedef vcp::matrix< double, vcp::pdblas > vmat;
int main(void){
int n = 5000;
vmat A, x, b;A.rand(n, n);b.ones(n, 1);
b = A*b;
vcp::time.tic();
x = lss(A, b);
vcp::time.toc();
std::cout << x(0) << std::endl;}

return 0;42

mpfr,kv,VCPを使った保証付き連立一次方程式#include <iostream>

int main(void){
int n = 200;
vmat A, x, b;#include <kv/interval.hpp>#include <kv/mpfr.hpp>#include <kv/rmpfr.hpp>

A.rand(n, n);b.ones(n, 1);
b = A*b;#include <vcp/matrix.hpp>#include <vcp/imats.hpp>#include <vcp/matrix_assist.hpp>

vcp::time.tic();
x = lss(A, b);
vcp::time.toc();#include <vcp/vcp_timer.hpp>
typedef kv::mpfr< 150 > kvmpfr;
typedef kv::interval< kvmpfr > impfr;
typedef vcp::imats< kvmpfr > policy;
typedef vcp::matrix< impfr, policy > vmat;

std::cout << x(0) << std::endl;}

return 0;43BLAS,kv,VCPを使った保証付き連立一次方程式#include <iostream>

int main(void){#include <kv/interval.hpp>#include <kv/rdouble.hpp>

int n = 5000;
vmat A, x, b;#include <vcp/matrix.hpp>#include <vcp/pidblas.hpp>#include <vcp/matrix_assist.hpp>

A.rand(n, n);b.ones(n, 1);
b = A*b;#include <vcp/vcp_timer.hpp>

vcp::time.tic();
x = lss(A, b);
vcp::time.toc();

typedef kv::interval< double > idouble;
typedef vcp::pidblas policy;
typedef vcp::matrix< idouble, policy > vmat;

std::cout << x(0) << std::endl;}

return 0;44