2015q3 Homework #1 Ext
目標
利用以下公式
\(\pi = 4 \int_{0}^{1} \frac{1}{1 + x^2} dx\)
採用離散積分的方法求圓周率，並著手透過 SIMD 指令作效能最佳化
Baseline
#include <stdint.h>
double compute_pi(size_t dt)
{
double pi = 0.0;
double delta = 1.0 / dt;
for (size_t i = 0; i < dt; i++) {
double x = (double) i / dt;
pi += delta / (1.0 + x * x);
}
return pi * 4.0;
}dt 為 [0, 1] 區間的大小，理論上 dt 切割地越細，計算結果越精確，但運算量自然越大。 請用 dt = 128M 作為輸入
AVX SIMD
#include <immintrin.h>
#include <stdint.h>
double compute_pi(size_t dt)
{
double pi = 0.0;
double delta = 1.0 / dt;
register __m256d ymm0, ymm1, ymm2, ymm3, ymm4;
ymm0 = _mm256_set1_pd(1.0);
ymm1 = _mm256_set1_pd(delta);
ymm2 = _mm256_set_pd(delta * 3, delta * 2, delta * 1, 0.0);
ymm4 = _mm256_setzero_pd();
for (int i = 0; i <= dt - 4; i += 4) {
ymm3 = _mm256_set1_pd(i * delta);
ymm3 = _mm256_add_pd(ymm3, ymm2);
ymm3 = _mm256_mul_pd(ymm3, ymm3);
ymm3 = _mm256_add_pd(ymm0, ymm3);
ymm3 = _mm256_div_pd(ymm1, ymm3);
ymm4 = _mm256_add_pd(ymm4, ymm3);
}
double tmp[4] __attribute__((aligned(32)));
_mm256_store_pd(tmp, ymm4);
pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
return pi * 4.0;
}參考 flops 去修改上述程式碼，使其得以在 GNU/Linux 編譯和運作，並且測試其效能改進
