線形補間の実装の比較 - 茂加部珈琲店

線形補完をいろんな方法で実装して、生成されるアセンブリを比較してみたいと思います

素直な実装

float lerp(float t, float v0, float v1) {
  return (1-t) * v0 + t * v1;
}

加減算が2, 乗算が2の実装です
t=0の時v0,t=1の時v1が保証されます（多分)

g++ 8.2 -O2

lerp(float, float, float):
  movss   xmm3, DWORD PTR .LC1[rip]
  subss   xmm3, xmm0
  mulss   xmm0, xmm2
  mulss   xmm1, xmm3
  addss   xmm0, xmm1
  ret

普通ですね
最初に1.fを読み込んでいます

clang++ 6.0 -O2

lerp(float, float, float):
  movss   xmm3, dword ptr [rip + .LCPI5_0] # xmm3 = mem[0],zero,zero,zero
  subss   xmm3, xmm0
  mulss   xmm3, xmm1
  mulss   xmm2, xmm0
  addss   xmm3, xmm2
  movaps  xmm0, xmm3
  ret

ほぼg++と同じ

次はFMA命令を有効化してみます

g++ 8.2 -O2 -march=native -ffp-contract=fast

lerp(float, float, float):
  vmovss  xmm3, DWORD PTR .LC0[rip]
  vsubss  xmm3, xmm3, xmm0
  vmulss  xmm0, xmm0, xmm2
  vfmadd231ss     xmm0, xmm3, xmm1
  ret

fmaが最後に入ってますね
精度的に若干有利なはず

clang++ 6.0 -O2 -march=native -ffp-contract=fast

lerp(float, float, float):
  vmovss  xmm3, dword ptr [rip + .LCPI5_0] # xmm3 = mem[0],zero,zero,zero
  vsubss  xmm3, xmm3, xmm0
  vmulss  xmm0, xmm0, xmm2
  vfmadd213ss     xmm3, xmm1, xmm0
  vmovaps xmm0, xmm3
  ret

ほぼ同じ

O3とかは危険な式変形やりだすので除外

精度を犠牲にした実装

float lerp(float t, float v0, float v1) {
  return v0 + (v1 - v0) * t;
}

乗算が減るけど、t=1のときにv0 + (v1 - v0)になるので結果がきっちりv1になる保証はありません

g++ 8.2 -O2

lerp(float, float, float):
  subss   xmm2, xmm1
  mulss   xmm2, xmm0
  movaps  xmm0, xmm2
  addss   xmm0, xmm1
  ret

普通ですね。clang++もほぼ同じなので省略します

次はFMAあり

g++ 8.2 -O2 -march=native -ffp-contract=fast

lerp(float, float, float):
  vsubss  xmm2, xmm2, xmm1
  vfmadd132ss     xmm0, xmm1, xmm2
  ret

除算とFMA命令のみ。はやそうです
ただt=1の時の誤差問題は解決できません
clang++も全く同じコードを吐きます

バランスの良い実装

float lerp(float t, float v0, float v1) {
  return t * v1 - (t * v0 - v0);
}

加減2,乗2と最初の実装と同じですが、順番が違います
t=0,t=1の結果がv0,v1になることが保証されます（多分）

g++ 8.2 -O2

lerp(float, float, float):
  mulss   xmm2, xmm0
  mulss   xmm0, xmm1
  subss   xmm0, xmm1
  subss   xmm2, xmm0
  movaps  xmm0, xmm2
  ret

普通ですね
clang++も同じコードを吐きます

FMAを入れてみると

g++ 8.2 -O2 -march=native -ffp-contract=fast

lerp(float, float, float):
  vfmsub132ss     xmm1, xmm1, xmm0
  vfmsub132ss     xmm0, xmm1, xmm2
  ret

FMA系命令2つだけになります
clang++でも同様

おまけ MSVCでやってみる

msvc 19.14 /O2 /fp:fast /arch:AVX2

lerp, COMDAT PROC
  vmulss   xmm3, xmm0, xmm2
  vmovaps xmm2, xmm1
  vfmsub213ss xmm2, xmm0, xmm1
  vsubss   xmm0, xmm3, xmm2
  ret      0

なんでやねん

まとめ

線形補完の実装は、精度と速度のトレードオフが存在する。
こんなに簡単なコードでも浮動小数点数が絡むと大変になる。浮動小数点数こわい。
最近のコンパイラはstd::fmaとか使わなくても勝手にやってくれる。
標準はstd::fmaだけですが、FMA系命令は他にも色々あります。