浮動小数点演算の演算順序

たけおか

初出: 2010/NOV/20

0.はじめに

浮動小数点演算では、演算の順序が重要である。という、よく知られている話を、 一応、書いておく。
特に、変わったことは述べていない。

このページは、x86 80bit演算問題の補足のために存在するだけである。


1. まだ、はじめに

一般人は、数値の計算は子供の頃からやっていて、
「コンピュータは、数値計算なんかできて当たり前だ」と思っている。
しかし、計算機は、有限の精度で演算や記憶を行っているため、 数値計算はそれほど簡単ではない。

例えば、ある数 a に、bを乗じて、その後 c で除するとき、
C言語で、 
x= a*b/c;
と書きそうになる。
しかし、C言語でそれは大丈夫だろうか?

実は大丈夫ではない。
軟弱 Cプログラマには、この文書を最後まで読むことをおすすめする。



2. 評価順序



Lispであろうが、Fortranであろうが、C言語であろうが、
評価の順序は重要である。
数値の四則演算でさえ、評価の順序が非常に大事である。

2.1 乗除算
のは、結果が異なる場合がある。

次のような2つのコード断片で実験する。
  1. 小さな数(1.000000e-308)に対して、まず、3e16を掛ける乗算を行い、次に、3e16で割る除算を行う。 
     double a,b,t;

 a=1.0/1e308;
 printf("a= %e\n", a);
 t= a*3e16;
 b= t/3e16;
 printf("b= %e\n", b);
 --
a= 1.000000e-308
b= 1.000000e-308
    特に問題なく、もとの小さな数にもどる。


  2. 小さな数(1.000000e-308)に対して、まず、3e16で除する除算を行う。次に、3e16を乗ずる乗算を行う。 
     double a,b,t;

 a=1.0/1e308;
 printf("a= %e\n", a);
 t= a/3e16;
 b= t*3e16;
 printf("b= %e\n", b);
 -- 
a= 1.000000e-308
b= 0.000000e+00
    小さな数に対して、先に除算をする、小さくなりすぎて、0になってしまう。
    それに、いくら掛けても、0のままである。
ということで、乗算と除算の順序には気をつけなればならない。

ここで、Cプログラマは、b= a*3e16/3e16; と書きたいだろうが…(後述)



2.2 加算

ただの加算でさえ、注意が必要である。

変数aに入った大きな数 1e36 (10の36乗)に、変数cに入った小さな数(1e19 (10の19乗))を1000000000回加える。という、 極めて簡単な計算を行う。
(10の19乗は、普通の感覚では、大きな数だが)

次のような2つのコード断片で実験する。
  1. 先に、0に対して、変数cに入った小さな数(1e19)を1000000000回加える。その合計と、変数aに入った数 1e36を加える。 
    double a,c,s,t;
a=1e36;
c=1e19;

for(t=0., i=0;i<1000000000;i++)
        t+=c;
s=a+t;
-- 
 結果= 1.0000000100000000e+36
    a+(c+...+c) という気分である。
    そして、得られる結果は、人間が期待するものである。


  2. 数 1e36が入った変数sに対して、変数cに入った小さな数(1e19)を 1000000000回加える。 
    double a,c,s;
a=1e36;
c=1e19;
s=a;
for(i=0;i<1000000000;i++)
        s+=c;

-- 
 結果= 1.0000000000000000e+36
    ((((a+c)+c)+...+c) という気分である。
    この計算の結果、変数sの内容はまったく変化していない。
    これは、人間が期待していない結果であろう。
上の例2番では、初期値としてsに入っている値(1e36)は、 1e19に対して大きすぎる。
よって、1e36 + 1e19 の加算を行ったとき、桁落ちが生じて、結果として、1e19 は 0として加算されてしまう。

大きな数と、小さな数をとりまぜて、加減算する場合、小さな数同士、また、大きな数同士を先に計算し、 その後、大きな数と小さな数の計算を行わなければならない。
そうしなければ、ここで示したように、桁落ちによって、期待しない結果を得ることになってしまう。



3. C言語の評価順序


3.1 C言語文法における評価順序の規定

さて、浮動小数点演算の評価の順序が、非常に重要なことは分かった。
そこで、C言語の話をする。

C 言語の文法では、強さが同じ演算子については、ソースコード上の字面とは無関係に、 コンパイラは演算の順序を自由に入れ替えることができる、と、規定されている。

括弧"( )" を付けるのも無意味である。
コンパイラは、括弧を外した後に、演算子の強さを見て、最適化を行おうとする。
括弧による、演算の結合は保存されるが、括弧を付けたからといって、その演算が他より先に実行されるとは、 まったく保証されない。
というか、評価の順序の制御に、括弧はまったく意味を持たない。

さらに、加えて言えば、二項演算子の左オペランドと右オペランドのどちらを先に評価するかも、規定されていない。
一般人の気持ち的には、"+"演算子の左を先に評価して欲しいが、 C言語ではそんなものは保証されない。

つまり、
・a/b*c/d*e は 乗算と除算を自由に入れ替える事ができる
・a+b+c+d-e-f-g-hは加算と減算を自由に入れ替える事ができる
・プログラマは‘( )’などを使用してみても計算順序を制御できない


余談になるが…
二項演算子のオペランドのどちらを先に評価するかが規定されていない、ということは、
x = (*p++ << 8)| *p; というようなコーディングが誤っていることを意味する。
"|"演算子の左オペランドの *p++が、右オペランドの*pの参照前か参照後のどちらに実行されるか、C言語では規定されていないからである。
このコードは、ある時、たまたま、うまく意図どおりに動作したとしても、 コンパイラの最適化オプションを変えたり、違うコンパイラに持っていったら、まったく動作しないことが、よくある。


3.2 C言語における評価順序の制御

それでは、C言語では、どのように記述すれば、意図した順序で、演算を実行させることができるのか?

例えば、
a=1.0/1e308;
b= a/3e16*3e16;
C言語では、この式はどんな結果を出すかは、不定である。
除算と乗算のどちらが先に実行されるか不定であるから


これの解決は、演算の分解である。
すべての浮動小数点演算を、一演算子、一文に分解する。

例えば、先の式は、次のようにして、評価順序を明示する。
一時変数を導入して、演算順序を明示。
 a=1.0/1e308;
 tmp= a/3e16;
 b= tmp*3e16;
一演算子あたり一文とするために、一時変数を用意し、それへどんどん代入しながら、演算を進めるように記述する。
現在のコンパイル技術では、無駄な変数への格納、参照は発生しない。
(ただし、x86で、64bit もしくは、32bit 精度の演算を厳密に行わせるには、無駄と思われる格納が必要となるのだが…
「x86における浮動小数点演算の精度の制御と、80bit 浮動小数点演算問題」を参照 )


s=a+b+c+d; と書きたいが…
各演算をバラして、下のように書く。 
s=a;
s+=b; s+=c; s+=d;

s=a*b/c*d; と書きたいが… 各演算をバラして、下のように書く。 
t=a*b;
t=t/c;
s=t*d
ということで、C言語の浮動小数点演算のあるプログラムは、一演算、一文、 として記述するのが、定石である。

また、x86 は、内部で 80bit 演算を行っている、という特殊事情から、
世の中に多くある 64bit の計算機と互換性のある、ポータブルな演算のためには、
一時変数への格納が必須である。
「x86における浮動小数点演算の精度の制御と、80bit 浮動小数点演算問題」を参照


4. Fortran

数値計算の王道の言語は、やはり Fortranである。

Fortranは、浮動小数点演算については、 ソースコードの字面の左から右に(人間の習慣にあった)演算を行う。
コンパイラが勝手に演算の順序を入れ替えることはない。

整数演算についてはその限りではない。つまり、 コンパイラが最適化によって、演算の順序を変更することがある。
整数は、演算の順序を変更しても、計算の精度が変わることがないからである。
(整数でも、除算と乗算の順序が入れ替わると、困ることがよくある。 Fortranであっても、それは文を分けて制御しなければならない)



たけおか(竹岡尚三)のホームページ