浮動小数点の誤差と対策

IEEE 754浮動小数点の基本

浮動小数点を理解するためには、データの内部構造を知ることが有効です。IEEE 754（2019年改訂）が浮動小数点演算の規格です。

単精度（binary32）・倍精度（binary64）の内部形式

はじめに用語を整理します。値の大きさを表す部分を指数部（exponent）、値の細かさを表す部分を仮数部（fraction）と呼びます。精度（precision）は先頭の隠れビットを含む有効ビット数、有効桁は10進での目安です。

単精度（IEEE 754 binary32、C言語のfloat相当）は符号1ビット、指数部8ビット、仮数（格納部）は23ビットです。先頭の1は正規化により暗黙的に扱われるため、精度は24ビットになります。倍精度（IEEE 754 binary64、C言語のdouble相当）は符号1ビット、指数部11ビット、仮数（格納部）は52ビットです。こちらも同様に、精度は53ビットです。

10進の有効桁は約7桁と約16桁で、扱う値の範囲により使い分けられます。センサ値の平均、フィルタ係数、制御ゲインなど、値が示す意味と必要な桁数を照らし合わせて選ぶことが重要です。

非正規化数・無限大・NaNの扱い

指数部がすべて0の非正規化数や、指数部が全1の±Inf・NaNは、実装によっては通常の演算より多くの時間を要することがあります。指数部が全1のとき、仮数部が0の場合は±Inf、非0の場合はNaNとなります。組み込みFPUでは、処理系によっては非正規化数の演算をソフトウェア側にフォールバックする実装もあり、演算時間が伸びるおそれがあります。

非正規化数を0として扱うflush-to-zero（FTZ）の設定を備えたFPUもあります。ただし有効にするとゼロ近傍の値が失われるため、精度への影響を確かめた上で選びます。ゼロ近傍のセンサ入力を長時間扱う処理では、値の一部だけ処理時間が長くなる可能性もあります。

この章の要点

• IEEE 754（2019年改訂）は、単精度と倍精度の内部形式を定めています
• 10進の有効桁は約7桁と約16桁で、扱う値の範囲で使い分けます
• 非正規化数やNaNは、実装によって演算時間が伸びることがあります

丸め誤差が生じる仕組み

丸め誤差は「2進数では表せない値がある」ことにより生じます。精度を上げれば解決するわけではないことに注意が必要です。

10進数の小数と2進数での近似（0.1問題）

10進数の小数は、2進数では近似値でだけ表現できることに注意が必要です。たとえば、0.1は2進数では循環小数となるため、桁数は無限に続きます。そのため、有限桁の浮動小数点では正確に表現できません。多くの環境では、0.1を10回加算しても厳密な1.0にはなりません。精度を上げれば誤差そのものは小さくなりますが、10進数の小数が2進数で有限桁に収まらないという原因は残ります。

正規化と丸めモード

浮動小数点の値は、仮数の先頭の桁が0にならない形へそろえて保持します。この形へそろえる処理を正規化と呼びます。0.1は0.01×10¹とも1.00×10⁻¹とも表記できますが（ここでは考え方を示すために10進数で例示します。IEEE 754が扱うのは2進数の指数表現です）、先頭の桁が0でないほうを正規化された表現と定めています。表現を1つに決めることで、仮数の桁を余さず使えます。

正規化しても表しきれない端数は丸めます。IEEE 754の既定はroundTiesToEven（最近接。ちょうど中間のときは下位ビットが偶数になる側へ）です。常に切り上げると誤差が一方向へ偏りますが、偶数側へ寄せると偏りが打ち消しあいます。

丸めモードはC言語からも切り換えられます。C99の標準ライブラリの<fenv.h>で、fesetround()により丸めモードを設定できます。ただし、切り換えると、既定のまま動く他の環境と結果が一致しなくなります。他環境と突きあわせて検証する場面では、既定のままが無難です。

この章の要点

• 10進数の小数は2進数で有限桁に収まらず、近似値でしか表せません
• 精度を上げても誤差が小さくなるだけで、表現できない原因は残ります
• 既定はroundTiesToEven（最近接偶数丸め）で、丸めモードを変えると他環境と結果がずれます

• テスト支援システムの商品情報はこちら
• PDFで詳しく見る

等価比較が不適切な理由と対処

浮動小数点同士を==で比べると、期待どおりに「一致した」結果を得られない場合があります。==は値を比較しますが、10進で同じ値のつもりでも、2進への変換と丸めで実際に格納される値がわずかに違ってしまうためです。値が異なれば==は偽を返します。

値の一致を問うかわりに、差が許容範囲に収まっているかを問う形に変えることで解決できます。許容範囲の決め方には、比率で見る相対誤差と、差そのもので見る絶対誤差の2通りがあります。

相対誤差ベースのepsilon判定

相対誤差とは、2つの値の差を値の大きさで割った比率です。この定義は、1991年にACM Computing Surveysへ掲載された論文が示しています。

David Goldberg, "What Every Computer Scientist Should Know About Floating-Point Arithmetic", ACM Computing Surveys, Vol.23, No.1, March 1991 です。同論文は相対誤差を「2つの数の差を実数値で割ったもの」と定義しています。

比率で見るのは、浮動小数点の誤差が値の大きさに比例するためです。仮数の桁数は値によらず一定なので、値が8倍になれば表せる刻みも8倍になります。同じ論文では、ある値を8倍したとき末尾の桁で数えた誤差は0.5から4.0へ8倍になる一方、相対誤差は0.8eps（ε：マシンエプシロン）のまま変わらない例をあげています。

この性質が相対誤差の利点です。センサから取得した生の値でも物理量に変換した後でも、同じ考え方を用いて算出できます。

判定はfabs(a - b) <= eps * fmax(fabs(a), fabs(b))と書きます。maxはC標準の関数ではないため、実装ではmath.hのfmaxか、自前のマクロを用います。左辺が2つの値の差、右辺が「大きいほうの値の何割まで許すか」です。epsには、1に足すと値が変わる最小の刻み（マシンエプシロン。C言語ではfloat.hのFLT_EPSILON／DBL_EPSILON）の数倍を置きます。

注意点は、値が0に近づくと右辺も一緒に小さくなり、判定が事実上の完全一致に近づくことです。実装では、比較本体の前に前処理を置きます。

• isnan(a) || isnan(b) のときはfalseとし、用途によっては別扱いにします。
• isinf(a) || isinf(b) のときは a == b で処理し、同符号のInfだけを一致とします。
• a == b のときはtrueとし、+0と-0も一致として扱います。

その上で本体を fabs(a - b) <= fmax(abs_eps, rel_eps * fmax(fabs(a), fabs(b))) とします。しきい値は例として、rel_eps = 4*FLT_EPSILON（float）、4*DBL_EPSILON（double）などを用います。

絶対誤差との併用パターン

絶対誤差とは、2つの値の差そのものの大きさです。値の大きさで割らないため、ゼロ近傍でも実用的な判定をおこなえます。しきい値を量の単位で決められる点も利点です。一例として、温度を「0.1℃未満の相違は同じ」とみなすことがあげられます。

弱点は相対誤差の裏返しです。扱う値の桁が大きく動くと、しきい値を決め直す必要が出ます。0.1という許容幅は、値が1のときは十分です。一方で値が10000の場合は、完全一致に近い値しか許容しなくなり、必要以上に厳しい許容幅となるおそれがあります。

そこで、ゼロ近傍は絶対誤差、それ以外は相対誤差を用いる形で両者を併用する手法が有効です。fabs(a - b) <= fmax(abs_eps, rel_eps * fmax(fabs(a), fabs(b)))と書けば、右辺は2つのしきい値の大きいほうを採ります。値が小さいうちはabs_epsによる絶対誤差の判定が働きます。値が大きくなると相対誤差を用いて判定します。

この章の要点

• ==は値を比較するため、2進変換と丸めで値がずれると期待どおりの結果になりません
• 値の大きさに比例する誤差には相対誤差、ゼロ近傍には絶対誤差が向きます
• 両者を併用し、NaNは先に別扱いで判定するのが実務的な形です

桁落ち（catastrophic cancellation）の実例

近い値同士の引き算の結果、有効桁が一気に失われる現象を桁落ちと呼びます。ここまで解説した誤差の蓄積とは異なり、単一の演算で起き得ることに注意が必要です。

二次方程式の解の公式での桁落ち

1991年に Computing Surveys へ掲載された論文のSection 1.4「Cancellation」が、この現象を示しています。(-b ± sqrt(b² - 4ac)) / (2a) の分子で、bとsqrt項の値が近いとき、片方の解の有効桁がほぼ消えます。教科書どおりの公式が、そのままでは実装向きではない典型例です。

この桁落ちを解決するには、有理化で式を書き換える方法があります。桁落ちしないほうの解をx₁とすると、もう一方はx₂ = c / (a x₁)で求められます。この形にすれば有効桁を保てます。

加算並べ替えでの回避

加算では、絶対値の小さい順に足すという単純な工夫により、誤差の蓄積を抑えられます。浮動小数点では、演算の順序が結果を変えることも押さえておきたいポイントです。

この章の要点

• 近い値同士の引き算では、単一の演算で有効桁が一気に失われます
• 二次方程式の解の公式は、有理化で式を書き換えれば有効桁を保てます
• 加算は絶対値の小さい順に足すと、蓄積する誤差を抑えられます

• テスト支援システムの商品情報はこちら
• PDFで詳しく見る

補正和による誤差蓄積の抑制

長時間の平均や数値積分をおこなうと、微小な誤差が積み重なり無視できない値に達するおそれがあります。W. Kahan氏の1965年論文「Further remarks on reducing truncation errors」（Communications of the ACM, Vol.8）で提案された補正和アルゴリズムが、この蓄積を抑える工夫です。

補正和アルゴリズムの原理

Kahan氏の原論文が示すとおり、各加算で失われる下位ビットを補正項として保持し次の加算に持ち越します。追加の変数一つとわずかな演算コストで、誤差の増加を抑えられます。

平均化・積分演算での適用

補正和が効くのは、加算の回数が多く、1回当たりの誤差が積み上がる場面です。同論文の定理8では、単純に足し上げると各項が最大nepsで揺らぐのに対し、補正和では揺らぎが2epsに収まると示しています。nは加算の回数です。単純和の誤差は回数に比例しますが、補正和は回数に依存しないことも特徴の一つです。

センサ値の長時間平均や制御ループの積分項など、数万回を超えるループほど差が出ます。数十回の加算では演算コストのほうが目立ちます。

この章の要点

• 補正和は、各加算で失われる下位ビットを次の加算へ持ち越します
• 単純和の誤差が回数に比例するのに対し、補正和は回数に依存しません
• 数万回を超えるループほど差が出て、数十回では演算コストが目立ちます

組み込み向けFPUの実装と単精度の限界

組み込みシステムで浮動小数点を扱う多くの場面は、マイコンに内蔵されたFPUで対応できますが、性能には違いがあります。制約を知っておくと選択が現実的になります。

単精度専用FPUと倍精度対応FPUの実装差

単精度専用FPUを搭載したマイコンでは、単精度演算だけをハードウェアで処理します。倍精度対応FPUを搭載したマイコンでは、単精度・倍精度の両方をハードウェアで処理できます。倍精度に対応するかどうかは製品構成により異なるため、データシートでの確認が必要です。マイコンベンダーが提供するDSP最適化ライブラリでも、環境により単精度FPU向けと倍精度FPU向けでビルドが分かれています。

倍精度エミュレーションのコスト

単精度FPUで倍精度演算を書くと、ソフトウェアエミュレーションに落ちて桁違いに遅くなります。倍精度FPUを搭載した構成を選ぶことで、このリスクを下げられます。その際は、DSP関数側も倍精度APIを選ぶ必要があります。

この章の要点

• 単精度専用FPUと倍精度対応FPUがあり、対応は製品構成で異なります
• 単精度FPUで倍精度を書くと、ソフトウェアエミュレーションに落ちます
• DSPライブラリも単精度向けと倍精度向けでビルドが分かれています

• テスト支援システムの商品情報はこちら
• PDFで詳しく見る

固定小数点との使い分け判断

浮動小数点はすべての用途に向くわけではなく、用途に応じた固定小数点との使い分けが重要です。

FPU非搭載マイコンでの現実解

FPU非搭載のマイコンでは、float型の四則演算がハードウェア命令に置き換わりません。コンパイラが用意する浮動小数点エミュレーションのライブラリ関数が呼ばれ、多数の積和演算に展開されます。演算コストは1桁以上跳ね上がります。この場合は、固定小数点を用いて組み直すほうが高速となる場合も多いです。

別記事「FPUなしで高精度な固定小数点演算の設計と実装」との連携

代表的な組み込み向けDSPライブラリが提供するq7・q15・q31固定小数点型のように、精度・値域・演算コストの3軸で比較すれば切り換えの判断が容易になります。詳細は別記事「FPUなしで高精度な固定小数点演算の設計と実装」で解説しています。浮動小数点は、内部形式と誤差の性質を把握し使い分けることで、有効な道具として使えます。

この章の要点

• FPU非搭載のマイコンでは、浮動小数点演算がライブラリ呼び出しになります
• 演算コストは1桁以上上がり、固定小数点のほうが高速なこともあります
• 精度・値域・演算コストの3軸で比較すると、切り換えの判断が容易です