float g_Const_ODE_Func[4] = { 1.0f, 2.0f, 3.0f, 4.0f };

float ODE_Func_C( float t, float x, float constant )
{
        return ( t * ( constant - x ) ); // (1)
}


void EulerKernel_C( float* p_t, float x[4], float StepSize )
{
        x[0] += ( StepSize * ODE_Func_C( *p_t, x[0], g_Const_ODE_Func[0] ) ); //
        x[1] += ( StepSize * ODE_Func_C( *p_t, x[1], g_Const_ODE_Func[1] ) ); // 
        x[2] += ( StepSize * ODE_Func_C( *p_t, x[2], g_Const_ODE_Func[2] ) ); //
        x[3] += ( StepSize * ODE_Func_C( *p_t, x[3], g_Const_ODE_Func[3] ) ); // (2)
        *p_t += StepSize; // (3)
}

#define ALIGNED         __declspec( align(16) ) // 16 byteに境界する  e.g. ALIGNED float x[4];
#define NAKED           __declspec( naked )

// Name         ODE_Func_SSE_v1()   
// Function     dx/dt = t * ( Const - x )
// Input
//              XMM0 : [temp] t
//              XMM1 : [temp] x
//              XMM2 : [temp] Const
// Output
//              XMM0 : t * ( Const - x )
NAKED void ODE_Func_SSE_v1()
{
        __asm {
                subps xmm2, xmm1
                mulps xmm0, xmm2
                ret
        }
}

ALIGNED float g_StepSize[4]         = { 0.01f, 0.01f, 0.01f, 0.01f }; // t[]のキザミ値
ALIGNED float g_Const_ODE_Func[4]   = { 1.0f,  2.0f,  3.0f,  4.0f  }; // 常微分方程式の特性を決定する定数


// Name         EulerKernel_SSE_v1()        
// Function     オイラー法で常微分方程式を４つ同時に計算する
// Input
//              XMM0 : [temp] ODE_Func_v1()のパラメータt
//              XMM1 : [temp] ODE_Func_v1()のパラメータx
//              XMM2 : [temp] ODE_Func_v1()のパラメータConst
//
//              XMM6 : [static] tの初期値がロード済みのレジスタ
//              XMM7 : [static] xの初期値がロード済みのレジスタ
// Output
//              XMM6 : 1ステップ進んだ時のt[]の値
//              XMM7 : 1ステップ進んだ時のx[]の値
NAKED void EulerKernel_SSE_v1()
{
        __asm {
                movaps xmm0, xmm6 // t[]をロード
                movaps xmm1, xmm7 // x[]をロード
                movaps xmm2, [g_Const_ODE_Func]

                call ODE_Func_SSE_v1

                mulps xmm0, [g_StepSize] //
                addps xmm7, xmm0         // x[] += ( g_StepSize[] * result[] )

                addps xmm6, [g_StepSize] // t[] += g_StepSize[]
                ret
        }
}

// Name         ODE_Func_SSE_v2()   
// Function     Dy = t * ( Const - x )
// Input
//              XMM0 : [static, read only] t
//              XMM1 : [static, read only] x
//              XMM2 : [static, read only] Const
// Outout
//              XMM7 : t * ( Const - x )
// Other
//              XMM4 : [temp] Constのコピーと復元に使用
NAKED void ODE_Func_SSE_v2()
{
        __asm {
                movaps xmm4, xmm2 // Constをxmm4にコピー
                subps  xmm2, xmm1
                mulps  xmm2, xmm0
                movaps xmm7, xmm2
                movaps xmm2, xmm4 // xmm2を復元
                ret
        }
}

// Name         EulerKernel_SSE_v2()        
// Function     オイラー法で常微分方程式を４つ同時に計算する
// Input
//              XMM0 : [static] ODE_Func_v2()のパラメータt
//              XMM1 : [static] ODE_Func_v2()のパラメータx
//              XMM2 : [static] ODE_Func_v2()のパラメータConst
//              XMM3 : [static] t[]のキザミ値StepSize
// Output
//              XMM0 : 1ステップ進んだ時のt[]の値
//              XMM1 : 1ステップ進んだ時のx[]の値
// Other
//              XMM7 : [temp] ODE_Func_SSE_v2()の計算結果
NAKED void EulerKernel_SSE_v2()
{
        __asm {
                call ODE_Func_SSE_v2

                mulps xmm7, xmm3 //
                addps xmm1, xmm7 // x[] += ( g_StepSize[] * result[] )

                addps xmm0, xmm3 // t[] += g_StepSize[]
                ret
        }
}

#include "CycleMeasure.h"
#pragma comment( lib, "CycleMeasure.lib" )

#define SAMPLE          100 // サンプル数
#define PRE_CALCULATE   100 // 1回のループで常微分方程式を計算する回数

#define PRINT_RESULT    \
                        printf("\t t = %6.2f, x = (%f,%f,%f,%f)\n", t[0], x[0], x[1], x[2], x[3] );

void PerformanceTest_EulerKernel_C()
{
        int i, j;
        CycleMeasure cm;
        UINT64 cycle[SAMPLE];
        ALIGNED float x[4] = { 0.5f, 0.5f, 0.5f, 0.5f };
        ALIGNED float t[4] = { 0.1f, 0.1f, 0.1f, 0.1f };

        for( i = 0; i < SAMPLE; i++ ) {
                cm.Begin();
                for( j = 0; j < PRE_CALCULATE; j++ )
                        EulerKernel_C( t, x, 0.01f );
                cm.End();
                cycle[i] = cm.GetCycle();
                // PRINT_RESULT
        }

        // printf()が上のループで影響を与えるため, ループを分離している
        printf("EulerKernel_C()のパフォーマンステストの結果\n");
        for( i = 0; i < SAMPLE; i++ )
                printf("%d [cycles]\n", cycle[i]);
}


void PerformanceTest_EulerKernel_SSE_v1()
{
        int i, j;
        CycleMeasure cm;
        UINT64 cycle[SAMPLE];
        ALIGNED float x[4] = { 0.5f, 0.5f, 0.5f, 0.5f };
        ALIGNED float t[4] = { 0.1f, 0.1f, 0.1f, 0.1f };

        __asm {
                movaps xmm6, [t] // t[]をXMM6へロード
                movaps xmm7, [x] // x[]をXMM7へロード
        }

        for( i = 0; i < SAMPLE; i++ ) {
                cm.Begin();
                for( j = 0; j < PRE_CALCULATE; j++ ) {
                        __asm {
                                call EulerKernel_SSE_v1
                                movaps [t], xmm6 // XMM6をt[]にストア
                                movaps [x], xmm7 // XMM7をt[]にストア
                        }
                }
                cm.End();
                cycle[i] = cm.GetCycle();
                // PRINT_RESULT
        }

        // printf()が上のループで影響を与えるため, ループを分離している
        printf("EulerKernel_SSE_v1()のパフォーマンステストの結果\n");
        for( i = 0; i < SAMPLE; i++ )
                printf("%d [cycles]\n", cycle[i]);
}


void PerformanceTest_EulerKernel_SSE_v2()
{
        int i, j;
        CycleMeasure cm;
        UINT64 cycle[SAMPLE];
        ALIGNED float x[4] = { 0.5f, 0.5f, 0.5f, 0.5f };
        ALIGNED float t[4] = { 0.1f, 0.1f, 0.1f, 0.1f };

        __asm {
                movaps xmm0, [t] // t[]をXMM1へロード
                movaps xmm1, [x] // x[]をXMM6へロード
                movaps xmm2, [g_Const_ODE_Func]
                movaps xmm3, [g_StepSize]
        }

        for( i = 0; i < SAMPLE; i++ ) {
                cm.Begin();
                for( j = 0; j < PRE_CALCULATE; j++ ) {
                        __asm {
                                call EulerKernel_SSE_v2
                                movaps [t], xmm0 // XMM0をt[]にストア
                                movaps [x], xmm1 // XMM1をt[]にストア
                        }
                }
                cm.End();
                cycle[i] = cm.GetCycle();
                // PRINT_RESULT
        }

        // printf()が上のループで影響を与えるため, ループを分離している
        printf("EulerKernel_SSE_v2()のパフォーマンステストの結果\n");
        for( i = 0; i < SAMPLE; i++ )
                printf("%d [cycles]\n", cycle[i]);
}