分享到plurk 分享到twitter 分享到facebook

版本 be00093eb43263a1c338bd81725955872a922d16

2015q3 Homework #1 Ext

目標

利用以下公式

\(\pi = 4 \int_{0}^{1} \frac{1}{1 + x^2} dx\)

採用離散積分的方法求圓周率,並著手透過 SIMD 指令作效能最佳化

Baseline

#include <stdint.h>
double compute_pi(size 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 編譯和運作,並且測試其效能改進

參考資訊