gistfile1.h

// Benchmark divpd and approximate division using cvt/rcp + Newton method.
#include <stdio.h>
#include <sys/time.h>
#include <xmmintrin.h>

#define NLOOP 100000
#define ITER 256

__m128d outputs[ITER];
__m128d values[ITER];

static inline __m128d rcpm(__m128d a)
{
    __m128d b = _mm_cvtps_pd(_mm_rcp_ps(_mm_cvtpd_ps(a)));
    b = _mm_sub_pd(_mm_add_pd(b, b), _mm_mul_pd(_mm_mul_pd(a, b), b)); 
    b = _mm_sub_pd(_mm_add_pd(b, b), _mm_mul_pd(_mm_mul_pd(a, b), b)); 
    return b; 
}

void
test_div(int n, double *outs, double *ins)
{
    int i;

    for (i = 0; i < n; i++) {
        outs[i] = 1.0 / ins[i];     // it will be mapped to divsd.
    }

}

void
test_sse(int n, __m128d *outs, __m128d *ins)
{
    int i;

    for (i = 0; i < n; i++) {
        outs[i] = rcpm(ins[i]);
    }

}

int
main(int argc, char **argv)
{
    int i;

    srand48(123L);
    for (i = 0; i < ITER * 2; i++) {
        ((double *)values)[i] = drand48();
    }

    struct timeval s, e;
    double divelap, ssedivelap;

    // Dry run.
    for (i = 0; i <  NLOOP; i++) {
        test_sse(ITER, outputs, values);
    }

    gettimeofday(&s, NULL);

    for (i = 0; i <  NLOOP; i++) {
        test_sse(ITER, outputs, values);
    }

    gettimeofday(&e, NULL);

    ssedivelap = (double)(e.tv_sec - s.tv_sec) +
          (double)(e.tv_usec - s.tv_usec) / (double)1.0e6;

    printf("[SSE approx div] %f secs\n", ssedivelap);

    // Dry run
    for (i = 0; i <  NLOOP; i++) {
        test_div(ITER * 2, (double *)outputs, (double *)values);
    }

    gettimeofday(&s, NULL);

    for (i = 0; i <  NLOOP; i++) {
        test_div(ITER * 2, (double *)outputs, (double *)values);
    }

    gettimeofday(&e, NULL);

    divelap = (double)(e.tv_sec - s.tv_sec) +
          (double)(e.tv_usec - s.tv_usec) / (double)1.0e6;

    printf("[Scalar div] %f secs\n", divelap);
    printf("SSE ver is %f times faster.\n", divelap / ssedivelap);


    return 0;

}