#include <stdlib.h>

#include <stdio.h>

#include <math.h>

// 各ベクトルアーキテクチャごとに...

#if defined( __VEC__ )

// AltiVec

    // 各ベクトルタイプに対応して float[4] 配列のベクトルタイプを準備

    typedef vector float vFloat;

    // 実際の各SIMD言語に仮想のSIMD言語をマップするマクロを定義する。

    // 行列乗算の場合、

    // 実行する必要がある作業は、

    // 2 つの命令セットアーキテクチャ(ISA)間で基本的に同じである。

    #define vSplat( v, i )  ({ vFloat z = vec_splat( v, i );

                                /* return */ z; })

    #define vMADD                   vec_madd

    #define vLoad( ptr )        vec_ld( 0, ptr )

    #define vStore( v, ptr )    vec_st( v, 0, ptr )

    #define vZero() (vector float) vec_splat_u32(0)

#elif defined( __SSE__ )

// SSE

    // ヘッダファイル「xmmintrin.h」では、インテルのCプログラミングインターフェイスに従って

    // SSEおよびSSE2を使用するためにC関数を定義する。

    #include <xmmintrin.h>

    // 各ベクトルタイプに対応して float[4] 配列のベクトルタイプを準備

    typedef __m128  vFloat;

    // 実際の各SIMD言語に仮想のSIMD言語をマップする

    // マクロを定義する。

    // iは即値でなければならないので、ここで

    // iをスタックベースのコピーにエイリアスし、それを4回複製するのは間違っている。

    #define vSplat( v, i )({ __m128 a = v; a = _mm_shuffle_ps( a,  a, \

                             _MM_SHUFFLE(i,i,i,i) ); /* return */  a; })

    inline __m128 vMADD( __m128 a, __m128 b, __m128 c )

    {

        return _mm_add_ps( c, _mm_mul_ps( a, b ) );

    }

    #define vLoad( ptr )    _mm_load_ps( (float*) (ptr) )

    #define vStore( v, ptr )    _mm_store_ps( (float*) (ptr),  v )

    #define vZero()             _mm_setzero_ps()

#else

// スカラー

    #warning To compile vector code, you must specify -faltivec,

                -msse, or both- faltivec and -msse

    #warning Compiling for scalar code.

    // 上記のベクトルバージョンで実現することを示す

    // 同等のスカラー

    // ベクトル、4つのスカラーを持つ構造体として宣言

    typedef struct

    {

        float               a;

        float               b;

        float               c;

        float               d;

    }vFloat;

    // 要素iをベクトル全体に展開し、それを返す

    #define vSplat( v, i )  ({ vFloat z; z.a = z.b = z.c = z.d =  ((float*)

                             &v)[i]; /* return */ z; })

 // 乗算と加算が融合した操作をサポートするアーキテクチャでその操作を実行する

    // 結果 = X * Y + Z

    inline vFloat vMADD( vFloat X, vFloat Y, vFloat Z )

    {

        vFloat result;

        result.a = X.a * Y.a + Z.a;

        result.b = X.b * Y.b + Z.b;

        result.c = X.c * Y.c + Z.c;

        result.d = X.d * Y.d + Z.d;

        return result;

    }

    // 所定のアドレスでスタートするベクトルを返す

    #define vLoad( ptr )  ( (vFloat*) ptr )[0]

    // ベクトルを所定のアドレスに書き込む

    #define vStore( v, ptr )    ( (vFloat*) ptr )[0] = v

    // ゼロでいっぱいのベクトルを返す

    #define vZero()     ({ vFloat z; z.a = z.b = z.c = z.

                            d = 0.0f; /* return */ z; })

#endif

// ベクトル行列乗算関数のプロトタイプ

void MyMatrixMultiply( vFloat A[4], vFloat B[4], vFloat C[4] );

int main( void )

{

    // vFloatタイプ（上で定義）は、4つの浮動小数を含む

    // ベクトルまたはスカラー配列である

    // そのため、これらはそれぞれが4×4の行列であり、C言語の格納順序で格納される

    vFloat          A[4];

    vFloat          B[4];

    vFloat          C1[4];

    vFloat          C2[4];

    int i, j, k;

    // A、B、C1、およびC2の要素へのポインタ

    float *a = (float*) &A;

    float *b = (float*) &B;

    float *c1 = (float*) &C1;

    float *c2 = (float*) &C2;

    // データを初期化する

    for( i = 0; i < 16; i++ )

    {

        a[i] = (double) (rand() - RAND_MAX/2) / (double) (RAND_MAX  );

        b[i] = (double) (rand() - RAND_MAX/2) / (double) (RAND_MAX  );

        c1[i] = c2[i] = 0.0;

    }

    // 総当たりの行列乗算を実行し、

    // これを使用して正確さをチェックする

    printf( "Doing simple matrix multiply...\n" );

    for( i = 0; i < 4; i++ )

        for( j = 0; j < 4; j++ )

        {

            float result = 0.0f;

            for( k = 0; k < 4; k++ )

                result += a[ i * 4 + k] * b[ k * 4 + j ];

            c1[ i * 4 + j ] = result;

        }

    // ベクトルバージョン

    printf( "Doing vector matrix multiply...\n" );

    MyMatrixMultiply( A, B, C2 );

    // 結果が正しいことを確認する

    // ここでは多少の丸め誤差を考慮する

    printf( "Verifying results...");

    for( i = 0 ; i < 16; i++ )

        if( fabs( c1[i] - c2[i] ) > 1e-6 )

            printf( "failed at %i,%i:%8.17g %8.17g\n",  i/4,

                         i&3, c1[i], c2[i] );

    printf( "done.\n" );

    return 0;

}