数値計算講義 第１ ２ 回 多倍長演算・ 精度保証計算

数値計算講義 第１ ２ 回 多倍長演算・ 精度保証計算

カ ーネン コ

金 子

晃

kanenko@mbk.nifty.com alexei.kanenko@docomo.ne.jp

http://www.kanenko.com/

多倍長演算

通常のプロ グラ ミ ン グ言語でサポート さ れて いる 変数の桁数以上の 精度で計算する こ と

無限精度演算と か任意精度演算と も いう

原理： ソ ロ バン を 繋げれば

部屋に納ま る 限り

いく ら でも 大き な桁数で 計算でき る

i.e.

変数を 繋げれば

メ モリ ーの許す限り

いく ら でも 大き な桁数で 計算でき る

KETA(1) KETA(2) KETA(3) KETA(4) KETA(5)

任意精度浮動小数の最も 簡単なフ ォ ーマッ ト

正規化さ れた浮動小数の仮数部を 十進４ 桁ずつ整数配列に格納する ． 小数点の位置を 冪指数と し て 記憶する

こ の他に

全体の符号と 仮数部の桁数を 記憶する 必要がある ので

結局

多倍長浮動小数

と は次のよ う な内容の整数配列

と なる ：

^{全体の符号}

A(2) · · ·

^指数部

十進

符号付き

A(3) · · ·

^{仮数部の現在の長さ}

使用し て いる ユニッ ト 数

A(4) · · ·

^{仮数部の最初の４ 桁}

十進

〜

※

言語の場合は

適当な構造体を 定義すれば

も っ と プロ グラ ムし 易く なる

多倍長浮動小数の演算

加法 小数点の位置を 大き い方の数に合わせて 加え る 例：

のと き

= 0.1234×10² + 0.04323×10²

桁揃え 立場によ り 二つに分かれる ：

(1)

小数点以下４ 桁に精度を 固定し て いる と き

こ のと き は

第２ の数を 丸めて 加え る ：

= 0.1234×10² + 0.0432×10² = 0.1666×10² (2)

精度はも う 少し と っ て あり

二つの数は正確

(i.e.

上記の部分の後ろ に

が精度の分だけ続いて いる

こ のと き は

仮数部の長さ を １ 増やし て 続け る ：

= 0.1234×10²+ 0.04323×10²= 0.16663000×10²

A(5)

には

を 入れる のが正し い

間違え て

を 入れる と

ができ て し ま う

繰り 上がり 処理 加法の結果ある ユニッ ト で５ 桁になっ たら 必要と なる ． 繰り 上がり は次々 波及する かも し れな いので， 汎用に書く のは結構大変．

cf.

全加算器の設計．

最上位のユニッ ト で繰り 上がり が起こ る と ， 指数の調整が必要と なる ． 例：

４ 桁ずつの格納法だと ， こ の変更は大変

(

すべて のユニッ ト で１ 桁のずら し が必要と な る

なので， 冪は

の倍数と し

= 0.00011110×10⁸

でアレ ン ジさ れたも のと し て し ま う のが簡単．

(

ただし こ の種の処理が続く と 有効桁数が損な われる と いう 欠点がある

減法 同様だが， 引けない場合は上から 借り て 来る 必要が有る ．

便法と し て ， 負になっ て も よ いから 対応する 部位から 引いて し ま い，

後で一括アレ ン ジする と いう 流儀も 有る ．

加減算だけが続く と き は

こ の方が高速

かつ１ 万回程度の演算数なら

各ユニッ ト にも 途中結果を 格納する 余裕がある

実際の計算では

データ は符号付き なので

加法・ 減法のど ち ら を 実行する かは

第２ オペラ ン ド の符号と 合わせて 決ま る

答の符号は最上位ユニッ ト が正規化後にも 正になら ないと き 反転する

こ の場合は正規化を やり 直さ ねばなら ない

乗法 指数部は単なる 和

仮数部は基本的に整数の乗算と 同様 例：

右の計算よ り

仮数部は

と なる

簡易アレ ン ジ法の下では

(1)

精度が４ 桁に指定さ れて いる と き は

答

精度が８ 桁以上のと き は

答

と なる

1234 x 5678 ---

6170 7404

8638 9872 ---

7006652

一般に

と し

の仮数の第

ユニッ ト を

と 記せば

C_k=

k

X

j=1

(Aj×B_k−j+1

の上位４ 桁

の下位４ 桁

と なる

ただし

右辺で添字が範囲外と なっ た因子は

と みなす

特に， 中央付近のユニッ ト では総ユニッ ト 数程度の加算が生ずる

よ っ て 乗算結果は一回毎にアレ ン ジする のがよ い

oooo oooo oooo oooo oooo oooo oooo oooo oooo oooo +

oooo oooo 0.

0.

0. oooo oooo

oooo oooo

除法 指数部は引き 算

仮数部の除法は仮商を 立て る 小学生の筆算式にやればでき る が

低速

反復法

特に

法の援用など

種々 の工夫が行われて いる が

いずれにし て も 他の演算に比し て 非常に時間がかかる

組み込み函数の実装 単精度や倍精度で

など の 組み込み函数がど のよ う に実装さ れて いる かを 調べ

真似を する

チュ ーニン グが進みすぎて いる も のは真似ができ な いので

基本に帰る

例：

sqrt(x) Newton

法を 適用し た通常のラ イ ブラ リ 函数の実装法が そのま ま 使え る

し かも

割り 算を 避けた工夫が多倍長では非常に有効

exp(x)

概算で

の形にする ：

log 10⌉, y=x−nlog 10, exp(y)

は

展開で求める

通常の倍精度浮動小数用ラ イ ブラ リ では

同様に

冪を 括り 出し た後

の計算には

倍精度が保証さ れた複雑な係数の近似多項式な ど が 使われる が

それは倍精度を 越え たら 使え な い：

Q^o

数値計算の中には

函数近似法と いう ト ピッ ク がある

C++

によ る 多倍長演算の実装

多倍長浮動小数のク ラ ス

を 定義し

その上の演算を 実装する と 演算子を オーバーロ ード でき る

多倍長数に対する 組み込み函数も 定義し て 追加でき る

等の

の特徴が使え る

こ の結果

倍精度浮動小数に対する のと ほと んど 同じ 表記で 多倍長数での計算が可能と なる ：

infprec A, B, C;

...

C=A*B-3*A+2*B/A;

....

高速に実装する には

多倍長浮動小数型の引数のコ ピーが 頻発し ないよ う に

サブルーチン 間のデータ の受渡し の書き 方を 工夫する 必要がある

ただの配列

なら

普通に引数に書けば

のアド レ スが渡る

参照渡し

し かし

が

構造体宣言さ れて いる と

デフ ォ ールト では

の内容が コ ピーさ れて サブルーチン に渡さ れる

こ れを 防ぐ には

void sub(infprec &A){

...

のよ う な書き 方で

参照渡し を 宣言する と よ い

(

ポイ ン タ ーではないので

指すアド レ スは変化せず

安全

更に

の内容を

の中で変更さ れたく 無い場合は

void sub(const infprec &A){

...

と する と

アド レ ス参照なのに

の中で

の要素の値が修正不可能と なる

Q^o C

言語流に

宣言を ポイ ン タ で行う こ と も 可能

こ のと き は

生成・ 破壊演算子でメ モリ を 確保

開放する よ う に書く

ま た

構造体の内部を

にし て

のセキュ リ ティ 機能を 使う べき である ．

既存の多倍長演算シ ステム

UBASIC

老舗だが

特に昔の

のアセン ブラ で 書かれて おり

著者に

への移植の意志が無かっ たので

今は使われない

PARI

ボルド ー大で開発さ れ

言語で書かれた

主に整数論の計算用 ツ ールだが

初等函数の多倍長浮動小数演算も サポート

libpari.a

の形でラ イ ブラ リ と し て も 提供さ れて おり

こ の中の函数を 自分のプロ グラ ムにリ ン ク し て 使え る 他

ソ ースが公開さ れて いる ので

必要な部分を 取っ て 来て 自分のプロ グラ ム に取り 込んだり も でき る

も 多倍長演算にこ れを 利用し て いる

GMP GNU MultiPrecision library. gcc

用のラ イ ブラ リ と し て

年頃から 出て いる

今後はフ リ ーの多倍長演算ラ イ ブラ リ の主流になる と 思われる

そ の他個人の作品 小生を 含めて 多く の人が作っ て 個人的に使っ て いる

最近のお勧め：

http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib/

数式処理ソ フ ト に組み込ま れた多倍長演算機能 主な数式処理ソ フ ト

など は

いずれも その中で精度を 任意に指定し た多倍長浮動小数の演算が可能

し かし

こ れら を

自分で書いた

のプロ グラ ムと リ ン ク する のは困難

GMP

の使用例

言語編

#include<stdio.h>

#include<stdlib.h>

#include<gmp.h> /*

ヘッ ダフ ァ イ ルの挿入

int main(void) {

int i;

mpf_t s,t; /*

型宣言

mpf_set_default_prec(1000);

mpf_init(s); /*

変数を ０ で初期化

mpf_init_set_str(t,"1e0",10); /*

変数を １ で初期化

は

for (i=1;i<=10000;i++){

mpf_add(s,s,t); /* s <- s+t */

mpf_div_ui(t,t,i); /* t <- t*i */

}

mpf_out_str(stdout,10,1000,s); /*

結果の出力

は

mpf_clear(t); /*

必須

mpf_clear(s);

return 0;

}

コ ン パイ ルは５ 階に

が標準でイ ン スト ールさ れたので次のよ う にする ：

ラ イ ブラ リ

は

言語用で

用は

こ のラ イ ブラ リ では， 他に

整数型

有理数型

を サポート

指令の詳細はディ レ ク ト リ 内のマ ニュ アル

を 参照せよ

8 π

の計算

古代エジプト  

⁴

=·· 3.16 Archimedes 223

71 < π < 22

7 (

内接正

角形の辺長

劉徽

世紀）

7 = 3.1428571· · ·, 157

50 = 3.14, 3927

1250 = 3.1416

祖沖之（

世紀） 約率

7 ,

密率

113 = 3.14159292· · ·, 3.1415926< π <3.1415927 Aryabhat.a (¯

イ ン ド ， ５ 世紀

Adriaan Anthoniszoon (1527–1607

）

年に

フ ラ ン ス

代数学の父

公式

∞

Y

n=2

cos π 2ⁿ =

r1 2

s 1 2 +1

2 r1

2 v u u t1

2 +1 2

s 1 2 +1

2 r1

2· · ·

を 用いて

年に

を 小数点以下

桁計算

Ludolph van Ceulen (

オラ ン ダ，

正

角形の周長を 用いて 小数点以下

桁計算

関孝和 小数点以下

桁

正

角形の周長を

桁ほど 計算し

補外によ り

桁を 正し く 求めた

単位円に内接する 正

角形の半周長は

2ⁿ, L2² = 2√ 2.

半角公式よ り

2ⁿ⁻¹ = 2 sin π 2ⁿcos π

2ⁿ

∴

cos π 2ⁿ

=· · ·= L2² n

Y

k=3

cos π 2^k

= 2√ 2

n

Y

k=2

cos π 2^k

√1 2

= 2

n

Y

k=2

cos π 2^k

n→∞

−→ π.

cos π

2ⁿ

は半角公式で漸化式が作れる ：

2ⁿ = r1

2 +1

2cos π

2ⁿ⁻¹, cos π 2² = 1

√2

正

角形の半辺長：

正

角形の半辺長：

正

角形の半辺長：

今

正

角形の半辺長が

N + c2

N² +· · ·

^{と いう} 漸近展開を 持つと 仮定し

c0=π

が半円周の長さ と すれば

2¹⁵+ c2

2³⁰ +· · · 3.14159265238659134580· · ·=c0+ c1

2¹⁶+ c2

2³² +· · · 3.14159265328899276527· · ·=c0+ c1

2¹⁷+ c2

2³⁴ +· · ·

こ の三式から

を 消去する と

6c0+ c1

2¹⁵ +· · ·= 8×3.14159265328899276527· · ·

−2×3.14159265238659134580· · · , 3c0+ c1

2¹⁵ +· · ·= 4×3.14159265238659134580· · ·

−3.14159264877698566948· · ·,

∴

3 ×3.14159265328899276527· · ·

−2×3.14159265238659134580· · · +¹₃ ×3.14159264877698566948· · ·

= 3.14159265358979323894· · ·

左辺の

^は

2⁴⁵

^{のオーダー}

cf:

正

角形の半辺長：

正

角形の半辺長：

A. Sharp (1653-1742) 1699

下記公式によ り 小数点以下

桁：

π

6 = Arctan√¹ 3 = √¹

3

1− 3¹·3 +₅¹

·3² −+· · ·

J. Machin (1680–1751) 1706

下記公式を 導き

桁計算：

π

4 = 4Arctan¹₅ −Arctan 1 239

= 4 1

5 − 1

3·5³ + 1

5·5⁵ −+· · ·

− 1

239 − 1

3·239³ + 1

5·239⁵ −+· · ·

※

松永良弼

小数点以下

桁

下記公式を 導き 利用：

π = 3

1 +₄¹_·²₆ +_4·¹₆²_·^·₈³_·²₁₀+_4·6¹_·8²_·^·₁₀³²_·^·₁₂⁵²_·14+· · ·

W.Shanks 1873 Machin

の公式によ り 小数点以下

桁計算

こ の計算があま り に圧巻だっ たので

以後電子計算機の発明

ま で誰も 挑戦し よ う と する 者は居なかっ た

Ludolph

の計算について

正

角形の半辺長：

3.14159265358979323846264338327950288395418471219027 cf. π

の真値：

3.14159265358979323846264338327950288419716939937510

ち ょ う ど

桁ま で真値なので

ど のよ う に答え たか興味深い！

松永の計算について

松永の級数の値：

70

項ま での和：

3.141592653589793238462643383279502884197169398014488669322 78

項ま での和：

3.141592653589793238462643383279502884197169399375088142023 79

項ま での和：

3.141592653589793238462643383279502884197169399375101484119 80

項ま での和：

3.141592653589793238462643383279502884197169399375104756842

松永は

微強 と 記し て いる ので

実際には

桁計算し たと みて よ いであろ う

本当に

項計算し たのか

それと も

関孝和の補外法のよ う な加速法を 援用し たのか ？

コ ン ビュ ータ によ る

の計算

1949

世界最初

の電子計算機

で

桁

時間

．

の値の

桁以降が誤っ て いる のを 発見

1961 D.Shanks-J.W.wrench

によ り

万桁

約９ 時間

によ り

万桁

約

時間）

1973 Guilloud-Bouyer

によ り

万桁

約

時間）

1981

三好

金田 

によ り

万桁

約

時間）

1983

後

金田 

で

万桁（ 約

時間）

1985 Gosper

によ り

万桁

約

時間？ ）

によ り

万桁

約

時間）

1988

金田

田村 

で２ 億桁（ 約

時間）

1989 G.V.Chudnovsky

で

億桁（ ２ ヶ 月？ ）

 自家製計算機で

億桁（ ？ ）

1995

高橋

金田 

で

億桁（ 約

時間）

1996 G.V.Chudnovsky

 自家製計算機で

億桁（ 一週間？ ）

高橋

金田 

で

億桁（ 約

時間）

1999

高橋

金田 

で

億桁（ 約

時間）

2002.12.6

金田

日立チーム

で

兆

億桁

約

時間

の計算は

の項式を 用いて いたが

現在では

算術幾何平均の改良型を 用いる ：

2

を 初項と する 二つの数列を

2 , b_n=p

a_n−1b_n−1

で定める と き

n→∞

2a_nb_n 1−

n

X

j=0

2^j(a²_j−b²_j)

精度保証計算

精度保証計算と は

実数

を 有限浮動小数で近似する と き

x0< x < x1

と 上下から 挟むこ と によ り

最終的な答の誤差の大き さ を 保証付き で与え る も の

昔は区間演算と 呼んだ

実現する には

非常に計算量がかかり

非実用的だっ た

粗い計算だと

演算の度に区間がど んど ん大き く な り

すぐ に挟んでも 意味の無いよ う な 大き さ に区間が広がっ て し ま う

最近の

は切り 捨て

切り 上げのモード 変更ができ る ので

それぞれのモード で１ 回ずつ同じ 計算を する だけで

普通の計算の２ 倍の計算量で精度保証が可能と なる

う ま く 使う と

数学の定理の実用的証明の道具と な る

の丸めモード ：

通常

は正し い計算値を 維持する ため

少し 多めの桁で計算し

結果を 返すと き に対応する 変数のフ ォ ーマ ッ ト に合わせて 丸める

デフ ォ ールト は四捨五入

指定のフ ォ ーマ ッ ト で表現可能な最も 近い値に丸める

こ の他に

切り 上げ

切り 捨て

０ に最も 近い値

の 計４ 種が選択可能

こ れら は実数の集合

から 表現可能な浮動小数の集合

_{への写像と し て}

でき ち んと 規定さ れて いる

Near(-x)=-Near(x), 15

Up(-x)=-Down(x),Down(-x)=-Up(x), x

が表現可能浮動小数なら

任意の丸め方

について

演算

の全て について

0

Down UpNear Tozero Down Up

Tozero Near

1 2

−1

sqrt

ま では

に規定がある が，

，

など ， 他の組み込み函数について は規定が無いので，

CPU

の丸めモード を

や

に変え ただけでは，

上下から 挟める かど う かは現状では期待でき な い．

=⇒

^{自分で書く か}

など のソ フ ト に組み込ま れたも のを 利用する し かない

例： 計算機では

を

回足し て も

になら ない

藤代先生の基礎講義

と いう のを 種々 のモード で実験し て みる

最近接丸め

四捨五入

方式

標準

のと き

s=1.000000000000000e-01 = 9A 99 99 99 99 99 B9 3F P10

i=1s=9.999999999999999e-01 = FF FF FF FF FF FF EF 3F

上への丸め

切り 上げ

方式を 指定し たと き

s=1.000000000000000e-01 = 9A 99 99 99 99 99 B9 3F P10

i=1s=1.000000000000001e+00 = 03 00 00 00 00 00 F0 3F

下への丸め

切り 捨て

方式を 指定し たと き

s=9.999999999999999e-02 = 99 99 99 99 99 99 B9 3F P10

i=1s=9.999999999999998e-01 = FE FF FF FF FF FF EF 3F

ケプラ ーの問題

数値計算を 証明に用いた最初の例

球を 空間にでき る だけ密に詰める 方法は

果物屋さ んがオレ ン ジを 積んでいる やり かただろ う

不思議なこ と に

最下層の一段の並べ方が図

でも 図

でも

その後に無駄無く 積み上げて でき た結果は３ 次元的には同じ にな る ！ 図

は

によ る その説明図

図

図

図

こ れよ り 詰めら れないこ と を

が予想し て 以来約

年振り に

が精度保証付き 数値計算によ り 証明し た

が有限個のパラ メ ータ の最小問題に変換し て いた

その最小解は

問題の解の空間充填率を 主要項と する 漸近式を 持つ

はその誤差項を 改良し

それを 厳密に評価する こ と によ り

理論的最小解と

予想の解が一致する こ と を 示し た

現実の

での精度保証計算

INTEL x86

の場合：

プロ グラ ム見本

CPU

のモード を 保持する レ ジスタ ー

で表さ れる

の第

ビッ ト が 丸めのモード を 制御し て いる ．

第

ビッ ト は精度を 制御．

下記の函数

イ ン ラ イ ン アセン ブラ

は

の名前で

に定義さ れて いる ので， こ れを イ ン ク ルード する だけでも よ い．

void setround(int mode){

__asm __volatile ("fldcw %0" : : "m"(mode)) }

int getround(int *mode)

__asm("fnstcw %0" : "=m" (*(mode))) return *mode;

}

#define FE_TONEAREST 0x0000

#define FE_DOWNWARD 0x0400

#define FE_UPWARD 0x0800

#define FE_TOWARDZERO 0x0c00

AMD64

互換の

ビッ ト マ シン では， 上のデータ は共通だが，

10

バイ ト 長の

で実験し ないと 差が見え ない．

PowerPC

の場合：

プロ グラ ム見本

言語用のラ イ ブラ リ 函数が用意さ れて いる ．

void fesetenv(const fenv_t *); /*

丸めモード 設定

void fegetenv(fenv_t *); /*

丸めモード 問い合わせ

#define FE_TONEAREST 0x00000000

#define FE_TOWARDZERO 0x00000001

#define FE_UPWARD 0x00000002

#define FE_DOWNWARD 0x00000003

roundppc.c

は現在のイ ン テルマ ッ ク では， コ ン パイ ルし て も 正常に 動かない． 昔の

でコ ン パイ ルさ れたバイ ナリ

が，

エミ ュ レ ーショ ン で正し く 動く ので， こ れで

を 味わっ て く ださ い．

本日の講義内容の自習課題

1 Kerosoft

社の多倍長演算ラ イ ブラ リ を 使っ て みる ．

inflib.f,infdec.f

がラ イ ブラ リ で

など が それを 利用し たプロ グラ ムである

コ ン パイ ル例：

g77 -c infdec.f inflib.f

g77 log.f inflib.o infdec.o -o log ./log

2 GMP

を 使っ た

プロ グラ ム例を 実行し て みる

コ ン パイ ルの仕方は

gcc gmptest.c -lgmp -o gmptest ./gmptest

3 CPU

の丸めモード を 変更する 実験を し て みる

は

で普通にコ ン パイ ルでき る

roundppc.c

の方はコ ン パイ ルでき ないので，

を 実行する だ けにせよ ．

PowerMac

を 持っ て いる 人は

を 持ち 帰っ て コ ン パイ ル・ 実行し て みよ

同様に，

ビッ ト の

用の

は５ 階ではコ ン パイ ルはでき る が

実行し て も 何も 変わら ない． こ れは，

が丸めを

ビッ ト でやっ て いる ため，

64

ビッ ト に落と し たと き に， みな同じ 値になっ て し ま う ためのよ う で ある ．

32

ビッ ト

機を 持っ て いる 人は

を 持ち 帰っ て

コ ン パイ ル・ 実行し て みよ

本日の範囲のプロ グラ ム練習問題

問題

社の多倍長演算ラ イ ブラ リ を 使っ て

の解を

桁近似計算せよ ．

問題

言語で次の

の級数を 実装し

を

桁計算せよ ：

5− 1

3·5³ +· · ·+ (−1)ⁿ

(2n+ 1)·5²ⁿ⁺¹ +· · ·

−4 1 239 −1

3 1

239³ +· · ·+ (−1)ⁿ 2n+ 1

1

239²ⁿ⁺¹ +· · · .

［ こ れは次問の多倍長演算ラ イ ブラ リ を 構築し なく て も 単独で解ける

］ 問題

を 利用し て 何か興味を 持っ た数を

桁計算し て みよ ．

(

例え ば

第５ 回の課題に有っ た

の正の解など

問題

多倍長数の加減乗算を 実行でき る ラ イ ブラ リ を

言語で作れ

それを 用いて

2

を

桁計算せよ

［ 第４ 回に紹介し た

法と 割り 算を 避ける 工夫を 組み合わせた

方法によ れば

多倍長数同士の除算を 使わずに計算でき る

］