CPP-QUIZ Q5 のこたえ、マクローリン展開による sin/cos の高速化に潜んでいた大きな誤差を生じるバグ
概要
この記事を書くきっかけは先日書いた記事 CPP-QUIZ from Unreal 2017 ( part-1/2; 出題編 ) の Q5. について、こたえとして掲載したコードについて、 Twitter にて @nekketsuuu さんから次のようなご指摘を頂いた事でした。ありがとうございます😃
@USAGI_WRP 再びすいません。ふと思い立って CPP QUIZ Q5 回答の sin の std::sin に対する誤差をとってみたところ、大きい所で4桁程度の誤差が出ているようでした。どこが誤差を生んでいるのかまでは調べていませんが、誤差に周期性があるようなので低角度域へシフトするところでしょうかね……？ pic.twitter.com/cGSv4806yQ
— 熱血 (@nekketsuuu) 2017年12月6日
どう見ても誤差が大きすぎるし（この検証の元となった出題では std::sin との誤差が最大でも 1.0e-6f に収まるようにしたいというものだった）、明らかに癖の強い周期性が見える。
ここ数日、お仕事などでまとまったプライベート時間をなかなか確保できず、数日お待たせしてしまいましたが、この問題の原因を調査し整理できる時間を確保できたので記事として整理する事にしました。
問題を確認する
「にゃーん🐾」
とりあえず、現状の問題、誤差の発生状況を把握するため、[0 .. 2π] について視覚的にプロットを把握するに十二分と思われる適当な分解能で、
std::sinを計算し、SinCosBを計算し、 sin 成分を取り出し、std::sinとSinCosBの sin 成分の差の絶対値を計算して縦軸へ、入力の角度値を横軸へプロット
した絵です。こうしてみると、
π/2ごとに周期的に誤差が増減し、- 最大誤差が
1.0e-6fどころか1.0e-4fを超えてしまっている
ことが明らかです。
問題のコード（＝「こたえ編で紹介した現行の UE4 実装と等価なコード」）
#define PI 3.14159265358979f std::tuple< float, float > SinCosB( const float angle_in_radians ) { float quotient = angle_in_radians * 0.5f / PI; if ( angle_in_radians >= 0.0f ) quotient = (float)( (int)( quotient + 0.5f ) ); else quotient = (float)( (int)( quotient - 0.5f ) ); float a = angle_in_radians - quotient * 2.0 * PI; float s = 0.0f; if ( a > 0.5f * PI ) { a = PI - a; s = -1.0f; } else if ( a < -0.5f * PI ) { a = -PI - a; s = -1.0f; } else s = +1.0f; float p = a * a; // [ 1 / 2!, 1 / 3!, 1 / 4!, .., 1 / 11! ] constexpr float f2 = 1.0f / 2.0f; constexpr float f3 = f2 / 3.0f; constexpr float f4 = f3 / 4.0f; constexpr float f5 = f4 / 5.0f; constexpr float f6 = f5 / 6.0f; constexpr float f7 = f6 / 7.0f; constexpr float f8 = f7 / 8.0f; constexpr float f9 = f8 / 9.0f; constexpr float f10 = f9 / 10.0f; constexpr float f11 = f10 / 11.0f; return std::make_tuple ( a * ( +1.0f + p * ( -f3 + p * ( +f5 + p * ( -f7 + p * ( +f9 * p * ( -f11 ) ) ) ) ) ) , s * ( 1.0f + p * ( -f2 + p * ( +f4 + p * ( -f6 + p * ( +f8 + p * ( -f10 ) ) ) ) ) ) ); }
誤差が 1.0e-6f 未満となる十分な精度が得られるまでマクローリン展開を計算しています。この実装は std::sin と std::cos を使った同様の計算よりも wandbox 環境では 1.45 倍高速に動作した、めでたしめでたし、というものでした。
ところが、どうやらこの SinCos は一見すると正しくマクローリン展開を計算しているようで、残念なバグが潜んだコードだったのです。
ﾅ ﾅﾝﾀﾞｯﾃｰ!!
Ω ΩΩ
問題について調査する
一般的に浮動小数点数の取り扱いで生じる幾つかの可能性が思いつきましたが、これと直感で断言できるほどの確信は無く、簡単なものから可能性を排除していく事にしました。
1. 単純に浮動小数点数型の精度を上げてみる
float 固定だった浮動小数点数の型を template 化して同様にプロットを眺めてみます:
template < typename T > constexpr auto pi = (T)3.14159265358979323846264338327950288; template < typename T > std::tuple< T, T > SinCosC( const T angle_in_radians ) { T quotient = angle_in_radians * (T)0.5 / pi< T >; if ( angle_in_radians >= (T)0 ) quotient = (T)( (int)( quotient + (T)0.5 ) ); else quotient = (T)( (int)( quotient - (T)0.5 ) ); T a = angle_in_radians - quotient * (T)2.0 * pi< T >; T s = (T)0; if ( a > (T)0.5 * pi< T > ) { a = pi< T > - a; s = (T)-1; } else if ( a < (T)-0.5 * pi< T > ) { a = -pi< T > - a; s = (T)-1; } else s = (T)+1; T p = a * a; // [ 1 / 2!, 1 / 3!, 1 / 4!, .., 1 / 11! ] constexpr T f2 = (T)1 / (T)2; constexpr T f3 = f2 / (T)3; constexpr T f4 = f3 / (T)4; constexpr T f5 = f4 / (T)5; constexpr T f6 = f5 / (T)6; constexpr T f7 = f6 / (T)7; constexpr T f8 = f7 / (T)8; constexpr T f9 = f8 / (T)9; constexpr T f10 = f9 / (T)10; constexpr T f11 = f10 / (T)11; return std::make_tuple ( a * ( (T)1 + p * ( -f3 + p * ( +f5 + p * ( -f7 + p * ( +f9 * p * ( -f11 ) ) ) ) ) ) , s * ( (T)1 + p * ( -f2 + p * ( +f4 + p * ( -f6 + p * ( +f8 + p * ( -f10 ) ) ) ) ) ) ); }
「( ﾟ∀ﾟ)o彡ﾟおっぱい！おっぱい！」
SinCosC の template の実引数に float （プロット: 青）, double （プロット: 橙）, long double （プロット: 灰） を与え、それぞれの浮動小数点数型の精度ごとに std::sin との差の絶対値を絵にしてみましたが、この結果からは「どうやら問題は浮動小数点数型の精度に少なくとも直接的には起因した問題ではない」という事が明瞭となりました。
ついでの蛇足調査(1/2); そもそも float vs. long double の値の誤差ってどれくらいでるものなの？
[0 .. 2π] について図示に必要十二分な分解能で float と double それぞれの値の差の絶対値がどれくらい生じるかの図も用意してみました。
・・・こんなもんです・ｘ・
ついでの蛇足調査(2/2); SinCos の前半部分が無いとどうなるのか？
SinCos の実装は大きくわけて2段階の工程となります。
[-π .. +π]に射影しつつ2π * quotient + reminderを経由して、最終的に-π/2 .. +πと cos 用の符号を保存- マクローリン展開に基いて必要十分な精度まで正弦と余弦を計算する
この前半部分をやらずに、直接後半のマクローリン展開へ入力された値を放り込むとどんな事が起こるのか、もののついでにグラフを作っておきました。
みょーん・・・と、なってしまい 2π どころか πあたりで既に実用性がかなり怪しい精度になっています。この問題を隠蔽するために最も誤差が小さく他の象限とは対称となるだけの [-π/2 .. π/2] に角度を落としてマクローリン展開部分を計算していたわけです。
そもそもマクローリン展開ってそんな誤差がでるものなの？
でません。少なくとも数学的には。そんなわけで、疑う先はマクローリン展開の実装コードではないかと絞り込めてきました。
そもそも、入力値が π/2 の場合、数学的に正確な正弦は 1 ちょうどです。しかし、どうやら SinCosB やそれを template 化した SinCosC の実装では、入力値が π/2 のときに誤差が最大値となっている事がわかっています。なので、初期値が 0 付近の値の計算や 0 付近の値に対する誤差の計算の信頼性が云々・・・という事でも無さそうです。
ちなみに、入力値が π/2 の場合の場合の SinCosB または SinCosC の出力値は float, double, long double の何れで計算した場合も、出力値は 0.999843 となり、 vs. std::sin の差の絶対値は 0.000157 となります。
これは浮動小数点数型の精度、マクローリン展開の精度の何れからも不可解な結果です。
そこで、実装の分かり易く愚直なマクローリン展開の正弦を計算するコードを書いて試してみます。
#include <iostream> #include <iomanip> template < typename T > constexpr auto pi = (T)3.14159265358979323846264338327950288; /// @brief オーバーフローしない限りの任意の階数でマクローリン展開による正弦を計算 /// @tparam N マクローリン展開する階数 /// @tparam T 入力および計算と出力に用いる浮動小数点数型 /// @param a 入力角度 [rad] /// @return 正弦値 [-] template < std::uint8_t N, typename T > auto sin( const T a ) { static_assert( N % 2 == 1 ); static_assert( N <= 19, "factorial be overflow" ); // 展開された級数を加算/減算し結果を得る変数 T result = a; // 符号は項ごとに反転する bool is_positive = false; // マクローリン展開を項ごと愚直に計算 for ( decltype( N ) n = 3; n <= N; n += 2, is_positive = ! is_positive ) { // 項の分数の分子の側となる指数を逐次愚直に乗算して結果を得る変数 T power = 1; // 項の分数の分母の側となる階乗を逐次愚直に乗算して結果を得る変数 std::uint64_t factorial = 1; // 項の分子と分母それぞれを愚直に乗算を繰り返して計算 for ( auto m = n; m; --m ) { power *= a; factorial *= m; } // 項の符号に基いて計算中の結果へ加算/減算し計算精度を高める if ( is_positive ) result += power / (T)factorial; else result -= power / (T)factorial; } // マクローリン展開を要求された階数だけ計算した結果を返す return result; }
やや丁寧すぎるくらいにコメントをソースへ入れたので文章での解説は省略します。この sin に float, double, long double の π/2 を入力して出力を眺めてみましょう。
int main() { constexpr auto n = 11; std::cout << std::fixed << std::setprecision( 32 ) << sin< n >( pi< float > / 2 ) << '\n' << std::fixed << std::setprecision( 32 ) << sin< n >( pi< double > / 2 ) << '\n' << std::fixed << std::setprecision( 32 ) << sin< n >( pi< long double > / 2 ) << '\n' ; }
n = 11 （問題のある SinCosB などと同じマクローリン展開の階数）の結果:
0.99999988079071044921875000000000 0.99999994374105094507854119001422 0.99999994374105087037701150576297
あっさり、すっきり、精度6桁確保できていますね・ｘ・
ちなみに、 n = 13 にすると、
0.99999994039535522460937500000000 1.00000000066278027510691117640818 1.00000000066278009003360033313257
さらに桁単位で誤差が縮みます。
つまり、 SinCosB のマクローリン展開の実装に問題がある可能性がとても濃厚に疑われます🐰
SinCosB の何が問題だったのか？
次の2つのコードは同じでしょうか、また、違いがあるとすればどちらが正しいマクローリン展開でしょう。
// SinCosB のマクローリン展開の最終的な計算の実装 auto s0 = a * ( 1.0f + p * ( -f3 + p * ( +f5 + p * ( -f7 + p * ( +f9 * p * ( -f11 ) ) ) ) ) ); // クイズ出題にあたり参考とした元の UE4 FMath::SinCos での実装 auto s1 = ( ( ( ( ( ( -f11 ) * p + f9 ) * p - f7 ) * p + f5 ) * p - f3 ) * p + 1.0f ) * a;
s0: 0.99984318017959594726562500000000 s1: 0.99999988079071044921875000000000
CPP-QUIZ from Unreal 2017 では、出題の元ネタとして UE4 の FMath ライブラリーの実装に基きつつ、 UE4 についてはまったく知らない方でもほぼ純粋に C++ における実装問題として楽しめるよう、またコードも私なりに読みやすいよう、 "工夫" したつもりでした。一般に、数式でもマクローリン展開は低い階数から項を並べて書くことが多いので、ややこしい記述となる実装部位について読者にも少しでも読みやすいようにと思い、参考とした UE4 FMath ライブラリーの実装では右側が低い階数となっていた実装を「これは問題ない、同じ動作となる "はず" 」と書き換えていたのです。
よーく、よーく、眺めると、 s1 では +f9 と次の項との計算はもちろん + 演算子となっている ところが、 s0 では * 演算子となっています。 s0 と s1 を比較しやすいよう、 s1 の記述を s0 のように、つまりマクローリン展開についてよく数式で表す際に用いられるように低い階数から並べると、次のようになります。
|--- Mistake! auto s0 = a * ( 1.0f + p * ( -f3 + p * ( +f5 + p * ( -f7 + p * ( +f9 * p * ( -f11 ) ) ) ) ) ); auto s1 = a * ( 1.0f + p * ( -f3 + p * ( +f5 + p * ( -f7 + p * ( +f9 + p * ( -f11 ) ) ) ) ) ); |--- CORRECT!
やれやれ・・・なんという結末の「なあんだ、そんなこと！」感。調査に昨晩少々と今朝、記事にしながら併せて90分くらい掛けてしまいました。コードの記述がややこしい実装をやむを得ずする場合には、こうしたごく単純な書き間違えによるバグを埋め込んでしまう可能性が高まります。今回の例はやむを得ず高速化したいという要求から単純ではない実装を試みる必要がありましたが、気持ちが緩み過ぎて、そうしたバグが埋め込まれてしまいやすい実装を行う上で十二分に必要な措置を講じ確実で正確な検証を行うべきところを怠った事が実質的な "バグ" だったと言えます。
修正版 SinCosB の正弦と std::sin の差の絶対値:
最大でも 1.0e-6 未満に収まるようになりました😃 より具体的には入力値が 1.51163 [rad] （ ≃ 86.61 [deg] ）などで最大の誤差 1.79e-7 が観測されます。
幸い、今回は実害の無いクイズで、なおかつ異常に気づいた @nekketsuuu さんからの指摘を基に、早々にバグを発見、修正する記事をこのように書けましたが、これが実際の製品に埋め込まれていた可能性や、保守担当者が私ではなかった場合（特に若い子にこれを対処させる必要性が生じた場合など）を考えると、なかなか、危ないところでした。反省しつつ、もっとコードを、そしてメタ的（≃試験やアサーション）に安全に、検証を怠らずに、さらにメタメタ的な素養（≃意識、習慣）を強化したいと思います。
それでは、この後日談ネタも含めて CPP-QUIZ from Unreal 2017 の "楽しさ" として頂ければ幸いです。🎅メリクリ🎄
