いきなり余談

現在作成中のアレで、GPU Accelerationにに手を出すのはもう少し足場を固めてからと思ってたけど、CPUを使うコードをGPUを使うコードに変更するのは手間というか、設計レベルから別物になるかもしれないし、どうせ最終的にはGPU使うだろうしということで、自分の中でGPUの汎用計算利用を解禁した。
ただ、OpenGLを使ったGPGPUのみでCUDAはまだ禁止ということにしてる。
現在作成中のアレ、もしくはその派生物を配布する可能性がないとは言えないので。

現在、複雑な形の大きな面光源があっても結構正確で高速な直接照明の計算をGPUを使って行うコードを実装してる最中。
例の箱のシーンの天井全部が光源であるようなシーンでも1秒未満で描画可能なものがすでに出来ている。
前回のCPUのコードでもかなり早かったけど、あれは実質2個の半球点光源で近似しただけ（さらに言うとライティングの計算がやはり間違ってた）であって、光源が複雑な形である場合、近似の精度を上げるには半球点光源の数を増やさなければならず、例えば数百とか数千とかに増やすとかなり遅くなると思われる。
GPUを使うコードの方は、その数千とかの半球点光源を現に数秒で処理可能となっている。
ただ影についてはまだ実装しておらず、面光源を扱うということで半影、いわゆるソフトシャドウを表現出来なきゃいけないんだけど、今のところその実装方法を3パターンぐらい考えててどれがベストか悩んでるところ。
出来たらここに説明を書く予定。
（前回、「ベイカーはレンダラと違って半影領域の全体から見た比率はそう多くはならないだろう」と書いたけど、複数の光源が離れたところにあったり、デカイ面光源があるとほとんど半影になりますね。）

本題

そのソフトシャドウの実装方法を考えてる中でふと疑問に思った点があったので、それについて今回は書く。
内容はタイトルの通り。

浮動小数点数を例えば0.1ずつインクリメントするとすぐ誤差が表面化するのはよく知られてるが、整数として扱うと、具体的には1.0ずつインクリメントするとどうなるか調べてみた。

        float f = 0.0f;
        for (long i = 0; i < 123456789; i++) {
            f += 1.0f;
        }
        printf("f = %f\n", f);

こんなコードをCで書いてみたんだけど、出力はなんと 16777216.000000（約1億2千万回のループなのに約1千7百万弱） となる。
しかしこれは0.1ずつインクリメントした場合のように誤差が蓄積した結果ではなく、16777216以下までなら誤差なくインクリメント出来る。

なんでこんな現象が起きるかだけど、一言でいうと情報落ちが発生してるせい。
具体的な機序は次のような感じ。

float（IEEE 754 32bit）は仮数部23bit。ケチ表現を利用してるので有効桁数は2進数で23 + 1。
これを一応2 ^ (23 + 1)としてみると10進数で16777216となる。
ここで「仮数部が16777216までしか表現出来ないのが原因」とするのは微妙に間違っており、浮動小数点の仕組みについて誤解を招きやすいと思う。
正規化されることにより、16777216.0fという大きな数字がそのまま仮数部に16777216という値を持つことは無いからである。
（非正規化数はざっくり言うとゼロ付近の微小な値を表現するためのもので今回の話とは関係ない。）

で、floatを加算する場合、まず絶対値が小さいほうの指数部が大きいほうの指数部に合わせられる。
さらに大きい方の指数部に合わせた分だけ小さい方の仮数部が右シフトされる。
それから加算が行われる。

1.0fと16777216.0fの例でいうと、それぞれビット列は以下の通りになる。
左から符号部1ビット、指数部8ビット、仮数部23ビット。各部分の間に空白を入れた。
分かりやすさのため、ケチ表現で省略された1ビット（正規化で常に1になる）を括弧で仮数部に付け足した。

       1.0f = 0 01111111 (1)00000000000000000000000
16777216.0f = 0 10010111 (1)00000000000000000000000

ちょっと説明を端折るけど、それぞれ指数部は2のべき乗を表しており、1.0fの場合で0乗相当（数値にすると 2 ^ 0 = 1）、16777216.0fの場合で24乗相当（数値にすると 2 ^ 24 = 16777216）となる。

これらが加算される場合、1.0fの指数部は24 – 0 = 24増やされる。
それに合わせて1.0fの仮数部が24回右シフトされる。

ケチ表現で省略された1ビットが24回右に移動していくイメージだが、仮数部は23bitしかないので24回も移動すると右から漏れてしまう。
言い換えると仮数部の情報が無くなってしまい、1.0fが0.0f的な状態になってしまう。

なので、16777216.0fに何回1.0fを足しても何も足してない状態と同じになる。
これは指数部の差が24以上ある場合に起こるので、指数部の差が24未満の場合、例えば16777216.0f + 2.0fは情報落ちも誤差も発生せず正確に計算される。

また以上の事から、-16777216.0f〜16777216.0fまでは1.0f刻みで全ての整数を誤差なしで表現出来る事も解る。
（±16777216.0f〜±16777216.0f*2までは2.0f刻み、±16777216.0f*2〜±16777216.0f*4までは4.0f刻み、以降指数部が許す限り続く。となることも解る。）

個人的に今考えてる用途では16777216.0fまでの整数を誤差なしで表現出来るなら全然問題ないけど、一千七百万弱なんてのは最近のコンピュータにとっては大きくない値だし、intやunsigned intの範囲と比べてもかなり小さいので、何らかの事情が無い限りはfloatを整数的に使うのはやめといた方がいいだろう。