今やってる仕事では、マイコンを使って科学技術計算をやるのだが、科学技術計算では浮動小数を使って幅広い桁や有効数字の計算を行うことが要求される。

今回の処理は、20ミリ秒間隔でやるため、その科学技術計算も20ミリ秒以内に収めなければいけない。

紆余曲折あって、一連の演算の固定小数版と浮動小数版ができたので時間測定を行った。時間測定は処理の開始と終了時に出力ポートに0とか1とかを書き出してデジタルオシロスコープ(データロガー)で測定。すると固定小数版では1ミリ秒で終わるところが浮動小数版では8ミリ秒もかかることがわかった。今書きながら8ミリ秒だったら20ミリ秒以内だからいいじゃんと思ったが、まあ、これからいろんなエラー処理などを入れたりなど考えたら少しでも時間の余裕があるほうがよいということにしておこう。

で、固定小数で演算していろいろ試しているうちに、どうもオーバフローすることが頻繁に起こってしまう。ある程度固定小数に関しては勉強したつもりだが、まだまだ実装に考慮が足りないようだ。この辺に関する資料が自分の調べ方に問題があるのかもしれないが見つからなかった。(すごく難しそうなのは見つかりましたが今の私の能力では理解できなかった)

とにかく私の系では2の補数表現を用いた固定小数で、加減乗除算を行っていかに正しく計算を行うかについて興味があるわけで、それについて考察した内容を記す。(前置き長い)

軽く自分の理解をまとめるために小さい数で考える。

3ビットの整数があったら0から7までを表現できる。これは符号なし整数。(cでunsigned)
これを2の補数表現で正数と負数も表現すると-4から3までの数を表す。これは符号付整数。(cでsigned)

binary unsigned singed
000 0 0
001 1 1
010 2 2
011 3 3
100 4 -4
101 5 -3
110 6 -2
111 7 -1

こんなイメージ。じゃあ今度は符号付の固定小数を考える。符号1bit、整数部分2bit、小数部分0bitの数のことを1+2+0bit固定小数ということにする。そうすると、signedの整数は1+2+0bit固定小数といってよい。

こんな感じ。

binary 1+2+0 1+1+1 1+0+2
011 3.0 1.5 0.75
010 2.0 1.0 0.50
001 1.0 0.5 0.25
000 0.0 0.0 0.00
111 -1.0 -0.5 -0.25
110 -2.0 -1.0 -0.50
101 -3.0 -1.5 -0.75
100 -4.0 -2.0 -1.00

これで1+2+0なら幅1.0・最大値3.0・最小値-4.0の実数、1+1+1なら幅0.5・最大値1.5・最小値-2.0の実数、1+0+2なら幅0.25・最大値0.75・最小値-1.00の実数が表現できることになる。

じゃあ次に演算に関して考える。加減算は2の補数表現を使っているのですごく楽チンだ。
足してオーバフローを無視すればよい。この考え方をはじめて知ったときには目から鱗だった。

乗除算がやっかいで、2進数の乗除算は簡単な筆算でできるといった情報をなんとなくそのまま鵜呑みにしていたけど、2の補数表現ではそれが成り立たないのだ!例えば、signedで2*2=4を考える。結果はsigned 3bitでは表現できないので、signed 6bit(3+3bit)にします。ちなみに6bitのsingedは、

binary signed unsigned
011 111 31 31
... ... ...
000 100 4 4
000 011 3 3
000 010 2 2
000 001 1 1
000 000 0 0
111 111 -1 63
111 110 -2 62
111 101 -3 61
111 100 -4 60
... ... ...
100 000 -32 32

な訳で、010*010=000 100は普通に筆算して求めることができる。では、2*(-2)=-4はというと、010*110=001 100となって111 100にはならないのである。ちなみに001 100は12で、こいつはどういう意味かというとunsignedの2*6=12に対応しているということで、簡単な筆算で乗算ができるのはunsignedのときで2の補数表現を使ったsignedではそれが成り立たないことがわかった。

でも、科学技術計算で正の数だけで計算することは難しいわけで、signed intとか使うと内部的には2の補数表現なわけで、その辺の性質の理解をした上でプログラムをしないとオーバーフローとかしてこれがもう、人が装着する機械の制御プログラムだったりするわけですから、暴走すると非常に危ないわけです。

また、上の考察からもう一つ発見があって、これはsingedでもunsignedでも同じだが、符号付3bitで-4から3の数が表現できるのに対して、その乗算結果は-12から9。この数値は5bitでは表現できず、6bit必要なのだが、6bit整数全ての範囲を使っていない。

binary signed
011 111 31
... ...
001 001 9 これより大きくならない
... ...
000 001 1
000 000 0
111 111 -1
... ...
110 100 -12 これより小さくならない
... ...
100 000 -32

これだと、最上位の符号ビットの次のビットも意味がなくなってしまい、(1+2+0)*(1+2+0)=1+(1)+4+0という風に1ビットの無駄ができてしまうようだ。ちなみにこの無駄ビットを除いた1+4+0ビットでも-16から15までを表現できるので、それよりもビットに空きができてしまうようだ。

さて、私の系では、intもlongも32ビット長なので、多倍長演算を導入しない限り、それより長いビット長の演算はできない。この多倍長演算ってのもなかなか面白そうなのだが、さらっと調べた感じでは速度がでないようなので却下して、何が何でも32bit長のintを使って実数計算をやってやらなければならない。

では、intで安全な実数計算を行うにあたって、演算の性質と何がこいつで表現できるのか?ついて考察する。

安全な演算ということを考えると、まずは、演算対象そのものがオーバフローしてないことを前提として、加減算は演算後にオーバフローしてなければよし。乗算は上の性質を考慮して1+15+0bitから1+0+15bitの数としておけば、乗算結果は1+(1)+30+0bitから1+(1)+0+30bitに収まるのでオーバフローなしに乗算可能だと思う。除算は被除数を<<16して除数で割って結果がオーバフローしてないかの確認をすればよかろう。

という訳で、intで-32768以上32767以下の数が安全に計算できる範囲で各演算ロジックも上で定義できた。

こいつで表現できる数をまとめておく。

表現 最大 最小
1+0+15 3.051*e-5 0.9999 -1.0
1+1+14 6.103*e-5 1.9999 -2.0
1+2+13 1.220*e-4 3.9998 -4.0
1+3+12 2.441*e-4 7.9997 -8.0
1+4+11 4.882*e-4 15.999 -16.0
1+5+10 9.765*e-4 31.999 -32.0
1+6+9 1.953*e-3 63.998 -64.0
1+7+8 3.906*e-3 127.99 -128.0
1+8+7 7.812*e-3 255.99 -256.0
1+9+6 1.562*e-2 511.98 -512.0
1+10+5 3.125*e-2 1023.9 -1024.0
1+11+4 6.25*e-2 2047.9 -2048.0
1+12+3 1.25*e-1 4095.8 -4096.0
1+13+2 2.5*e-1 8191.7 -8192.0
1+14+1 5*e-1 16383.5 -16384.0
1+15+0 1 32767 -32768

これを見ながら実際の計算にどの精度の表現を使うかを決めて演算すればよいハズ。