Intel MKLのSIMD exp関数を読んだ

Intel CPUで高速に動く、数学関数のSIMD実装に Intel Math Kernel Library というものがあります。

それに含まれる指数関数実装を読んでみました。

AVX-512（512bitレジスタがあるので、doubleを最大で八並列で演算できる）用の実装は、以下のようになっていました*1。

00000000004013d0 <__svml_exp8_z0>:
  4013d0:  f3 0f 1e fa                     repz nop     edx
  4013d4:  62 f1 7c 48 10 25 62 29 00 00   vmovups      zmm4 , [rip+0x2962]
  4013de:  62 f1 7c 48 10 0d 98 29 00 00   vmovups      zmm1 , [rip+0x2998]
  4013e8:  62 f1 7c 48 10 15 ce 29 00 00   vmovups      zmm2 , [rip+0x29ce]
  4013f2:  62 f1 7c 48 10 1d 04 2a 00 00   vmovups      zmm3 , [rip+0x2a04]
  4013fc:  62 f1 7c 48 10 2d 7a 2a 00 00   vmovups      zmm5 , [rip+0x2a7a]
  401406:  62 71 7c 48 10 2d b0 2a 00 00   vmovups      zmm13, [rip+0x2ab0]
  401410:  62 f2 fd 78 a8 e1               vfmadd213pd  zmm4 , zmm0 , zmm1 , {rz-sae}
  401416:  62 f1 7c 48 10 35 e0 2a 00 00   vmovups      zmm6 , [rip+0x2ae0]
  401420:  62 71 7c 48 10 0d 16 2b 00 00   vmovups      zmm9 , [rip+0x2b16]
  40142a:  62 71 7c 48 10 1d 4c 2b 00 00   vmovups      zmm11, [rip+0x2b4c]
  401434:  62 71 7c 48 10 25 02 28 00 00   vmovups      zmm12, [rip+0x2802]
  40143e:  62 71 dd 18 5c f1               vsubpd       zmm14, zmm4 , zmm1 , {rn-sae}
  401444:  62 72 dd 48 7f 25 32 28 00 00   vpermt2pd    zmm12, zmm4 , [rip+0x2832]
  40144e:  62 d2 ed 18 bc c6               vfnmadd231pd zmm0 , zmm2 , zmm14, {rn-sae}
  401454:  62 d2 e5 18 bc c6               vfnmadd231pd zmm0 , zmm3 , zmm14, {rn-sae}
  40145a:  62 f1 fd 48 54 3d dc 29 00 00   vandpd       zmm7 , zmm0 , [rip+0x29dc]
  401464:  62 71 c5 18 59 c7               vmulpd       zmm8 , zmm7 , zmm7 , {rn-sae}
  40146a:  62 72 d5 18 b8 ef               vfmadd231pd  zmm13, zmm5 , zmm7 , {rn-sae}
  401470:  62 72 cd 18 b8 cf               vfmadd231pd  zmm9 , zmm6 , zmm7 , {rn-sae}
  401476:  62 71 bd 18 59 d7               vmulpd       zmm10, zmm8 , zmm7 , {rn-sae}
  40147c:  62 72 bd 18 a8 df               vfmadd213pd  zmm11, zmm8 , zmm7 , {rn-sae}
  401482:  62 52 bd 18 a8 e9               vfmadd213pd  zmm13, zmm8 , zmm9 , {rn-sae}
  401488:  62 52 ad 18 a8 eb               vfmadd213pd  zmm13, zmm10, zmm11, {rn-sae}
  40148e:  62 52 9d 18 a8 ec               vfmadd213pd  zmm13, zmm12, zmm12, {rn-sae}
  401494:  62 d2 95 18 2c c6               vscalefpd    zmm0 , zmm13, zmm14, {rn-sae}
  40149a:  c3                              ret

この実装では、 $\exp(x)$ を求めるとき、まず $\frac{x}{\log 2}$ を1/16の位まで求めます。この値の整数部分は、最後のvscalefpdで使用されます。小数部分は、vpermt2pdを使った16エントリの表引きに使われます。

$x$ を $\frac{1}{\log 2}$ で割ったあまりを $t$ とすると、絶対値は $\frac{\log 2}{16}$ 未満のはずなので、 $\exp(t)$ を六次の多項式を使って近似します。

vandpdは、 $x$ の絶対値が大きすぎて $\frac{x}{\log 2}$ を1/16の位まで求める計算がおかしくなった時に $t$ がおかしな値となり、最終結果がおかしくなることを防ぐ役割があるようです。 vandpdで使われるマスク定数は0xbfffffffffffffffであり、絶対値が2以上の数値やNaN・±Infを絶対値が2未満の数に変更する効果があります。 $x$ の絶対値が大きすぎる場合はどうせvscalefpdで0か+Infになるので、適当に小さくすればよいということのようです。

こうして見てみると、vperfmt2命令やvscalef命令は指数関数を実装するためにあるという感じの命令ですね……。

*1:マスクレジスタを使わないことが確定している場合の実装です。