#include <stdlib.h>

double Rand_0to1()
{
        return ( (double)rand() / (double)RAND_MAX );
}

double Rand_i0toi1( double i0, double i1 )
{
        double rand_0to1;

        rand_0to1 = (double)rand() / (double)RAND_MAX;

        return i0 + ( i1 - i0 ) * rand_0to1;
}

typedef struct {
        double x0, x1;
        double y0, y1;
        double z0, z1;
} INTERVAL;

double KamabokoFunc( double x, double y )
{
        return sin(x);
}

double CalcVolume_MC_HitOrMiss( INTERVAL I, int sample )
{
        double x, y, z, V0, V;
        int in = 0;

        // 直積空間[I.x0, I.x1] x [I.y0, Iy1] x [I.z0, I.z1]の体積
        V0 = ( I.x1 - I.x0 ) * ( I.y1 - I.y0 ) * ( I.z1 - I.z0 );

        for( int i = 0; i < sample; i++ ) {
                x = Rand_i0toi1( I.x0, I.x1 );
                y = Rand_i0toi1( I.y0, I.y1 ); // (これは不要)
                z = Rand_i0toi1( I.z0, I.z1 );
                if( z <= KamabokoFunc(x, y) )
                        in++;
        }

        // [I.x0, I.x1] x [I.y0, Iy1]の領域のf1の体積
        V = V0 * (double) in / (double) sample;

        return V;
}

double CalcVolume_BasicMC( INTERVAL I, int sample )
{
        double x, y, BaseArea, V = 0;

        BaseArea = ( I.x1 - I.x0 ) * ( I.y1 - I.y0 ) / (double)sample;

        for( int i = 0; i < sample; i++ ) {
                x = Rand_i0toi1( I.x0, I.x1 );
                y = Rand_i0toi1( I.y0, I.y1 );
                V += KamabokoFunc( x, y );
        }

        V *= BaseArea;

        return V;
}

double CalcVolume_Midpoint( INTERVAL* pI, int div_count )
{
        double width_x, width_y, x = 0, y = 0, sum_x, sum_y;

        width_x = ( pI->x1 - pI->x0 ) / ( double )div_count;
        width_y = ( pI->y1 - pI->y0 ) / ( double )div_count;

        sum_x = 0;
        sum_y = 0;

        for( int j = 0; j < div_count; j++ ) {
                for( int i = 0; i < div_count; i++ ) {
                        x = pI->x0 + ( i + 0.5 ) * width_x;
                        sum_x += fabs( KamabokoFunc( x, y ) );
                }
                sum_x *= width_x;
                sum_y += sum_x;
                sum_x = 0;
                y = pI->y0 + ( j + 0.5 ) * width_y;
        }

        sum_y *= width_y;

        return sum_y;
}