丸め誤差

突然ですが,C言語で

#define x1 0.5;
#define x2 0.1;

float sum1 = 0;
float sum2 = 0;

int i;

for(i = 0;i < 2;i++){sum1+=x1;}
for(i = 0;i < 10;i++){sum2+=x2;}

printf("%.16f\t%.16f\n" , sum1 , sum2);

を実行した結果を考えてみてほしい.

x1は0.5固定で,sum1に2回足される.つまり,sum1=0.5+0.5=1.0

x2は0.1固定で,sum2に10回足される.つまり,sum2=0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1=1.0

 

両方同じ結果だと思いますよね.それは,数学の話.

実行結果は次のようになります.

1.0000000000000000 1.0000001192092896

これは,コンピュータが本質的に孕んでいる問題で,C言語であろうと他の言語であろうと起こります.正確には,IEEE754という規格を採択しているとこういうことが起こります.

コンピュータで扱う実数は10進数ではなく2進数で,例えば0.1は

0.1_{(10)}=0.0001100110011001100\cdots _{(2)}

と無限に続く循環小数です.コンピュータのメモリは有限なので,どこかで打ち切らなければなりません.そこで,メモリになるだけ多くの小数を格納することを考えます.

さっきの例だと,最初の4つ続いている0はこれだけで4bitsも消費してしまい,bitの無駄です.指数は2^{-4}だから,-4_{(10)}=1100_{(2)}です(補数表現).この例では,減らせるbitはありませんでしたが,2^{-10}とかだと0が10個続く(10bits)のが,-10_{(10)}=1010_{(2)}なのでやっぱり4bitになり,6bitの短縮になりますね.つまり,2の指数を考えれば大幅にbitを減らせます.

すると,1.10011001100110011001100_{(10)}\cdots \times 2^{-4}となります.以降小数部分だけ考えます.

上の考え方を使えば,小数点より左は常に1なので,これにbitを割く必要はありません.これでもう1bit削れました.

さて,ここからがいよいよ打ち切りの問題.どこまで残すか.

IEEE754では単精度浮動小数点数(single, C言語ではfloat型)は小数部23bitsと定められていますので,次のように打ち切られます(|で打ち切られる)

0.1_{(10)} \to 1.\overbrace{10011001100110011001100}^{23{\rm bits}}|11\cdots _{(2)} \times 2^{ - 4}

そして,この入りきらなかった数字”11″は「丸め」られます.といっても磨くわけではありません.

IEEE754では,5種類の丸め方が定義されているのですが,最も用いられているのが(デフォルトで使われるのが),「最近接偶数丸め」です.すごい名前ですね.

「最近接偶数丸め」でググったらヒットするページのほとんどが10進数の四捨五入の話をしていて,コンピュータ内部の2進数処理とは異なっているので,本質を欠いている!!と思ったのが,このエントリを書こうと思った理由です.

「最近接偶数丸め」とは,基本的には零捨一入でちょうど中間の値の時には丸めた結果が偶数になる方をとる,という説明が2進数ではよくされます,が,よく分かりませんね.具体的には,

\cdots 00|00 \to 00

\cdots 00|01\to 00

\cdots {\bf 00|10 \to 00}

\cdots 00|11 \to 01

\cdots 01|00 \to 01

\cdots 01|01 \to 01

\cdots {\bf 01|10 \to 10}

\cdots 01|11 \to 10

という具合で丸めます.気をつけて欲しいのは,太字のところで,これが唯一,零捨一入と異なるところです.最近接偶数丸めでは,最下位ビットから2つ下のビットまで読み,それが”10″のときは,丸めた時の最下位ビットが0になるように丸めます.上の”10″が入りきらなかった場合,強制的に繰り上げるのが零捨一入です.

零捨一入では,最下位ビットのその1つ下の数しか見ないで繰り上げ,繰り下げを判断するのですが,これだと誤差が大きくなり,丸める前の和と丸めた後の和で誤差が発生します.(この誤差をバイアスと言います)

これは話題から逸れるのであまり詳しくは書きませんが,零捨一入では,例えば

1|00 \to 1

1|01 \to 1

0|10 \to 1

0|11 \to 1

ですが,これらの丸める前の和は11|10で,丸めた後の和は100です.

一方,最近接偶数丸めでは

1|00 \to 1

1|01 \to 1

0|11 \to 1

ですが,これらの丸める前の和は11|00で,丸めた後の和も11です.

 

冒頭のプログラムにおいて,0.1を10回足した方がうまくいかなかったのは,0.1を単精度浮動小数点数で表した時に丸め誤差が発生し,それが蓄積されていった,ということです.逆に,0.5を2回足した方がうまくいったのは,0.5_{(10)}=0.1_{(2)}なので単精度浮動小数点数で表しても正確に表現でき,丸める必要がなかったためだと言えます.

さて,0.1がうまくいかなかったのは丸め誤差の蓄積だということは説明しましたが,倍精度浮動小数点数(double型)ではどうでしょうか.

実行結果は,

1.0000000000000000 0.9999999999999999

となります.

double型というのは,小数部を52bits格納できます.大事なことは,23bitsみたいに奇数でなく,偶数ビット格納できるということです.これにより,0.1_{(10)} = 0.000\dot{1}10\dot{0}という4bitを周期とする循環小数がきちんと格納されます.

広告

論理回路クイズ

面白い問題を聞いた.

出典はこちら

●問題●
AND,OR,NOTゲートを使い、
入力X,Y,Zに対して各々の否定を出力する
三入力三出力の論理回路を作ってください。
ただし、NOTゲートは二つしか使えません。

最後の条件の使い所が超難しい.

自分なりの答えは以下.

“論理回路クイズ” の続きを読む

RSA暗号の数理

RSA暗号を最短で解説してみたい.

暗号学の習わしにならって,暗号はAliceとBobの間で交換されるとする.そして,第三者Eveが盗聴しているとする.

Aliceは

1. 2つの大きな素数p,qを選ぶ.

2. 合成数n=p\times qを計算する.

3. \phi(x)をオイラー関数として,\phi(n)=\phi(pq)=(p-1)(q-1)を計算する.

4. \phi(n)以下で,これと互いに素な数kを求める.

5. ku-\phi(n)v=1を満たすu,vを求める.

6. Bobにn,kを伝える.(これは公開通信路でEveにも傍受できる.)

Bobは

7. 伝えたい平文のメッセージMを用意する.

8. M^k \equiv b \ \ \ (mod\ n)を計算し,bを求める.

9. Aliceにbを伝える.(これも公開通信路でEveにも傍受できる.)

Aliceは

10. b^u\equiv M \ \ \ (mod\ n)を計算し,不思議なことにMを得る.

この摩訶不思議な操作の理論的な根拠は次の定理に依る.

gcd(b,n) = 1, gcd(k,\phi(n)) = 1 を満たすb,k,nが与えられたとき,

M^k \equiv b\ \ \ (mod\ n) の解を求める.

1. 一次方程式 ku - \phi(n)v = 1 を満たす正の整数解u,vを求める.

2. b^u \equiv M\ \ \ (mod\ n)を計算して得られた値がMとなる.

つまり,b^u \equiv (M^k)^u \equiv M^{ku} \equiv M^{1+\phi(n)v} \equiv M\times (M^{\phi(n)})^v \equiv M\ \ \ (mod\ n)

というだけの話.

ちなみに,M^{\phi(n)}\equiv 1\ \ \ (mod\ n)はオイラーの定理(数論)です.

これと,もう一つStep5の背景には次の定理

a,bを0でない任意の整数とするとき,

方程式 ax+by = gcd(a,b) は常に整数解を持つ.

が使われている.

そして,Eveはb,k,nを知っている.

しかし,上記定理のuを求める過程で,\phi(n)が必要だということが分かるが,\phi(n)nの素因数分解が分かっていないと計算できない.

だから,これだけではMを得ることはできない.

では,実際コンピュータはどうするのか.

Step1の巨大素数を見つける段階.ここでは,まず例えば5000以下ぐらいの雑魚な素数をエラトステネスの篩で取り除く.続いて,フェルマー素数判定法(フェルマーの小定理)の対偶を用いて,「素数っぽいヤツ」を見つける.しかしここではまだ分からない.カーマイケル数などがあるからだ.最終的には,ミラーラビン素数判定法で決定する.

Step4,Step5の互いに素な数を見つけるには,ユークリッドの互除法で最大公約数が求まるので大丈夫.

Step8,Step10の巨大数の巨大数乗は,繰り返し二乗法を使う.これは,冪の数を2進数で表して,底の2,4,8,16,,,2^n乗を各々計算して最後に掛ける.これをすると,計算が省略できる.

最後に,

k,nを公開鍵と言い,uを秘密鍵という.公開鍵は平文Mを暗号化するだけの鍵であり,秘密鍵は暗号bを復号化するだけの鍵である.

このような暗号のやり取りを公開鍵暗号方式という.(これに対して,1つの鍵で暗号化と復号化を行う方式を共通鍵暗号方式という.)