近年(もう随分前になるのかな?)のプロセッサは除算をニュートン・ラフソン法で行っているらしいので、以前作った備忘録：除算命令を使わない除算を参考に divps の代替を作ってみます。
以下の環境
macOS Big Sur (11.5.2)
Intel(R) Core(TM) i7-8700B CPU @ 3.20GHz
gcc version 11.2.0 (GCC)
でコンパイルして実行します。
#include <math.h>
#include <inttypes.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <time.h>
#include <xmmintrin.h>
/*
* DIVPS
*/
inline static __v4sf divps1(__v4sf a, __v4sf b)
{
return _mm_div_ps(a, b);
}
/*
* DIVPS の代替
*
* x[1] = RCPPS(b)
* x[n+1] = x[n] * (2 - b * x[n])
*
*/
inline static __v4sf divps2(__v4sf a, __v4sf b)
{
__v4sf two = { 2, 2, 2, 2 };
__v4sf x1 = _mm_rcp_ps(b);
__v4sf x2 = x1 * (two - b * x1);
return a * x2;
}
/*
* 以下、動作確認用
*/
inline static uint64_t RDTSC()
{
#if defined(_MSC_VER) || \
(defined(__clang__) && \
(defined(__amd64) || \
defined(__i386) || \
defined(__x86) || \
defined(__x86_64)))
return __rdtsc();
#elif defined(__GNUC__) && \
(defined(__amd64) || \
defined(__i386) || \
defined(__x86) || \
defined(__x86_64))
uint32_t eax, edx;
__asm__ volatile ("RDTSC" : "=a" (eax), "=d" (edx));
return ((uint64_t)edx << 32) | eax;
#else
timespec c;
clock_gettime(CLOCK_REALTIME, &c);
return ((uint64_t)c.tv_sec * 1000 * 1000 * 1000 + c.tv_nsec);
#endif
}
static void print_tsc(uint64_t tsc, uint64_t cnt)
{
uint64_t i = tsc / cnt;
uint64_t f = (100 * (tsc % cnt) + (cnt / 2)) / cnt;
printf("%" PRIu64 ".%02" PRIu64, i, f);
}
int main()
{
__v4sf a = { 7, 7, 7, 7 };
size_t i, fcount = (1 << 20);
size_t vcount = fcount / 4;
float *buffer = (float *)_mm_malloc(2 * fcount * sizeof(float), 16);
__v4sf *p = (__v4sf *)buffer;
__v4sf *q = (__v4sf *)(buffer + fcount);
uint64_t tsc1, tsc2, tsc3;
if (!buffer)
{
fprintf(stderr, "_mm_malloc 失敗\n");
return 2;
}
for (i = 1; i <= fcount; i++)
{
buffer[i - 1] = (float)i;
buffer[i - 1 + fcount] = (float)i;
}
tsc1 = RDTSC();
for (i = 0; i < vcount; i += 8)
{
p[i + 0] = divps1(a, p[i + 0]);
p[i + 1] = divps1(a, p[i + 1]);
p[i + 2] = divps1(a, p[i + 2]);
p[i + 3] = divps1(a, p[i + 3]);
p[i + 4] = divps1(a, p[i + 4]);
p[i + 5] = divps1(a, p[i + 5]);
p[i + 6] = divps1(a, p[i + 6]);
p[i + 7] = divps1(a, p[i + 7]);
}
tsc2 = RDTSC();
for (i = 0; i < vcount; i += 8)
{
q[i + 0] = divps2(a, q[i + 0]);
q[i + 1] = divps2(a, q[i + 1]);
q[i + 2] = divps2(a, q[i + 2]);
q[i + 3] = divps2(a, q[i + 3]);
q[i + 4] = divps2(a, q[i + 4]);
q[i + 5] = divps2(a, q[i + 5]);
q[i + 6] = divps2(a, q[i + 6]);
q[i + 7] = divps2(a, q[i + 7]);
}
tsc3 = RDTSC();
_mm_free(buffer);
printf("divps1=");
print_tsc(tsc2 - tsc1, vcount);
printf(", divps2=");
print_tsc(tsc3 - tsc2, vcount);
printf("\n");
for (i = 0; i < 10; i++)
{
printf("%2d: (%f,%f,%f,%f) (%f,%f,%f,%f)\n",
(int)i,
p[i][0], p[i][1], p[i][2], p[i][3],
q[i][0], q[i][1], q[i][2], q[i][3]);
}
return 0;
}
実行結果
$ gcc -O3 -W -Wall sample.c && ./a.out
divps1=2.99, divps2=3.61
0: (7.000000,3.500000,2.333333,1.750000) (7.000000,3.500000,2.333333,1.750000)
1: (1.400000,1.166667,1.000000,0.875000) (1.400000,1.166667,1.000000,0.875000)
2: (0.777778,0.700000,0.636364,0.583333) (0.777778,0.700000,0.636364,0.583333)
3: (0.538462,0.500000,0.466667,0.437500) (0.538462,0.500000,0.466667,0.437500)
4: (0.411765,0.388889,0.368421,0.350000) (0.411765,0.388889,0.368421,0.350000)
5: (0.333333,0.318182,0.304348,0.291667) (0.333333,0.318182,0.304348,0.291667)
6: (0.280000,0.269231,0.259259,0.250000) (0.280000,0.269231,0.259259,0.250000)
7: (0.241379,0.233333,0.225806,0.218750) (0.241379,0.233333,0.225806,0.218750)
8: (0.212121,0.205882,0.200000,0.194444) (0.212121,0.205882,0.200000,0.194444)
9: (0.189189,0.184211,0.179487,0.175000) (0.189189,0.184211,0.179487,0.175000)
divps1 と divps2 の数値は、概ね __v4sf での割り算一回分（メモリ・アクセスを含む）の CPU サイクル数のハズです。アセンブリ出力を見ると、divps2 の呼び出し部は、divps1 の呼び出し部より２倍くらい長いのですが速度的には大差がありません。
アセンブリ出力（抜粋）
L5:
movaps %xmm0, %xmm1
subq $-128, %rax
divps -128(%rax), %xmm1
movaps %xmm1, -128(%rax)
movaps %xmm0, %xmm1
divps -112(%rax), %xmm1
movaps %xmm1, -112(%rax)
movaps %xmm0, %xmm1
divps -96(%rax), %xmm1
movaps %xmm1, -96(%rax)
movaps %xmm0, %xmm1
divps -80(%rax), %xmm1
movaps %xmm1, -80(%rax)
movaps %xmm0, %xmm1
divps -64(%rax), %xmm1
movaps %xmm1, -64(%rax)
movaps %xmm0, %xmm1
divps -48(%rax), %xmm1
movaps %xmm1, -48(%rax)
movaps %xmm0, %xmm1
divps -32(%rax), %xmm1
movaps %xmm1, -32(%rax)
movaps %xmm0, %xmm1
divps -16(%rax), %xmm1
movaps %xmm1, -16(%rax)
cmpq %rax, %rbx
jne L5
L6:
movaps (%rax), %xmm2
subq $-128, %rax
rcpps %xmm2, %xmm4
movaps %xmm2, %xmm3
movaps %xmm1, %xmm2
mulps %xmm4, %xmm3
subps %xmm3, %xmm2
mulps %xmm4, %xmm2
mulps %xmm0, %xmm2
movaps %xmm2, -128(%rax)
movaps -112(%rax), %xmm2
rcpps %xmm2, %xmm4
movaps %xmm2, %xmm3
movaps %xmm1, %xmm2
mulps %xmm4, %xmm3
subps %xmm3, %xmm2
mulps %xmm4, %xmm2
mulps %xmm0, %xmm2
movaps %xmm2, -112(%rax)
movaps -96(%rax), %xmm2
rcpps %xmm2, %xmm4
movaps %xmm2, %xmm3
movaps %xmm1, %xmm2
mulps %xmm4, %xmm3
subps %xmm3, %xmm2
mulps %xmm4, %xmm2
mulps %xmm0, %xmm2
movaps %xmm2, -96(%rax)
movaps -80(%rax), %xmm2
rcpps %xmm2, %xmm4
movaps %xmm2, %xmm3
movaps %xmm1, %xmm2
mulps %xmm4, %xmm3
subps %xmm3, %xmm2
mulps %xmm4, %xmm2
mulps %xmm0, %xmm2
movaps %xmm2, -80(%rax)
movaps -64(%rax), %xmm2
rcpps %xmm2, %xmm4
movaps %xmm2, %xmm3
movaps %xmm1, %xmm2
mulps %xmm4, %xmm3
subps %xmm3, %xmm2
mulps %xmm4, %xmm2
mulps %xmm0, %xmm2
movaps %xmm2, -64(%rax)
movaps -48(%rax), %xmm2
rcpps %xmm2, %xmm4
movaps %xmm2, %xmm3
movaps %xmm1, %xmm2
mulps %xmm4, %xmm3
subps %xmm3, %xmm2
mulps %xmm4, %xmm2
mulps %xmm0, %xmm2
movaps %xmm2, -48(%rax)
movaps -32(%rax), %xmm2
rcpps %xmm2, %xmm4
movaps %xmm2, %xmm3
movaps %xmm1, %xmm2
mulps %xmm4, %xmm3
subps %xmm3, %xmm2
mulps %xmm4, %xmm2
mulps %xmm0, %xmm2
movaps %xmm2, -32(%rax)
movaps -16(%rax), %xmm2
rcpps %xmm2, %xmm4
movaps %xmm2, %xmm3
movaps %xmm1, %xmm2
mulps %xmm4, %xmm3
subps %xmm3, %xmm2
mulps %xmm4, %xmm2
mulps %xmm0, %xmm2
movaps %xmm2, -16(%rax)
cmpq %rax, %rdx
jne L6
逆数なら、a を掛ける必要がなくなったり、演算精度を落としてもいいのなら n2 の計算を省ぶことも考えられます。すると、divps を使うより良くなる余地が出てきます。ということで rcpps/rcpss が用意されているようです。(rsqrtps/rsqrtss も同様)