метод Монте-Карло

//Трехмерное тело образовано объединением нескольких сфер произвольного размера 
//и взаимного расположения. Найти объем этого тела, используя метод Монте-Карло.
#include <iostream>
#include <cmath>
#include <vector>
#include <algorithm>
#include <functional>
#include <limits>
#include <numeric>
#include <ctime>
#include <cstdlib>
 
enum  T_coord_ind  
{
    X = 0,
    Y,
    Z,
    TOTAL_COORDS
};
 
T_coord_ind& operator++(T_coord_ind&  coord_ind)
{
    return coord_ind = T_coord_ind(coord_ind + 1);
}
 
typedef std::vector<double> T_coords;
 
class T_3_point
{
    T_coords  coords_;
public:
    typedef T_coords::iterator        iterator;
    typedef T_coords::const_iterator  const_iterator;
    typedef T_coords::value_type      value_type;   
 
    T_3_point
        (            
            value_type  x = 0,
            value_type  y = 0,
            value_type  z = 0
        ) : coords_(TOTAL_COORDS)
    {
        coords_[X] = x;
        coords_[Y] = y;
        coords_[Z] = z;
    }    
 
    iterator begin()
    {
        return coords_.begin();
    }
 
    iterator end()
    {
        return coords_.end();
    }
    
    T_3_point operator-(T_3_point  p) 
    {
        T_3_point  res_p;
        std::transform
            (
                begin(), 
                end(), 
                p.begin(),                
                res_p.begin(),                                
                std::minus<value_type>()
            );
        return  res_p;     
    }
 
    value_type& operator[](T_coord_ind  coord_ind)
    {
        return coords_[coord_ind];
    }
 
    value_type operator[](T_coord_ind  coord_ind) const
    {
        return coords_[coord_ind];
    }
 
    T_3_point::value_type  distance(T_3_point  p) const
    { 
        T_3_point  vector_from_this_in_p = p - *this;
 
        return sqrt
            (
                std::inner_product
                    (
                        vector_from_this_in_p.begin(), 
                        vector_from_this_in_p.end(),
                        vector_from_this_in_p.begin(),                        
                        0.0
                    )
            );
    }
 
    void input()
    {
        std::cout << std::endl
                  << "x = ";
        std::cin >> coords_[X];
        std::cout << "y = ";
        std::cin >> coords_[Y];
        std::cout << "z = ";
        std::cin >> coords_[Z];
    }
};//class T_3_point
 
typedef std::vector<T_3_point>  T_3_points;
 
class T_sfera
{
    T_3_point              center_;
    T_3_point::value_type  radius_;
public:
    T_sfera
        (
            T_3_point              center = T_3_point(), 
            T_3_point::value_type  radius = 0
        ) : center_(center), radius_(radius) 
    {}
 
    bool  is_vnutri_sfery(T_3_point  p) const
    {
        return  center_.distance(p) <= radius_;
    }
 
    void input()
    {
        std::cout << std::endl
                  << "Введите данные сферы:"
                  << std::endl
                  << "координаты центра:";
        center_.input();
        std::cout << "радиус: ";
        std::cin >> radius_;
    }
 
    T_3_point::value_type  get_min_coord(T_coord_ind  coord_ind)
    {       
        return center_[coord_ind] - radius_;
    }
 
    T_3_point::value_type  get_max_coord(T_coord_ind  coord_ind)
    {       
        return center_[coord_ind] + radius_;
    }
};//class T_sfera
 
typedef std::vector<T_sfera>  T_sfery;
 
//////////////////////////////////////////////////////////////////////////
class T_get_min_coord
{
    T_coord_ind            coord_ind_;
    T_3_point::value_type  min_coord_;
public:
    T_get_min_coord
        (
            T_coord_ind  coord_ind
        ) 
        : coord_ind_(coord_ind), min_coord_(std::numeric_limits<T_3_point::value_type>::max())
    {}
 
    void operator()(T_sfera sfera)
    {        
        if(sfera.get_min_coord(coord_ind_) < min_coord_)
        {
            min_coord_ = sfera.get_min_coord(coord_ind_);
        }
    }
 
    operator T_3_point::value_type()
    {
        return min_coord_;
    }
};//class T_get_min_coord
//////////////////////////////////////////////////////////////////////////
class T_get_max_coord
{
    T_coord_ind            coord_ind_;
    T_3_point::value_type  max_coord_;
public:
    T_get_max_coord
        (
            T_coord_ind  coord_ind
        ) 
        : coord_ind_(coord_ind), max_coord_(std::numeric_limits<T_3_point::value_type>::min())
    {}
 
    void operator()(T_sfera sfera)
    {        
        if(sfera.get_max_coord(coord_ind_) > max_coord_)
        {
            max_coord_ = sfera.get_max_coord(coord_ind_);
        }
    }
 
    operator T_3_point::value_type()
    {
        return max_coord_;
    }
};
//////////////////////////////////////////////////////////////////////////
class T_get_random_point
{
    T_3_point  min_XYZ_point_;    
    T_3_point  max_XYZ_point_;
public:
    T_get_random_point
        (
            T_3_point  min_XYZ_point,
            T_3_point  max_XYZ_point
        ) : min_XYZ_point_(min_XYZ_point), max_XYZ_point_(max_XYZ_point) 
    {}
    
    double  rand_between_0_and_1()
    {
        return static_cast<double>(rand()) / static_cast<double>(RAND_MAX);
    }
 
    T_3_point operator() ()
    {
        T_3_point::value_type  x = min_XYZ_point_[X] 
                       + rand_between_0_and_1() * (max_XYZ_point_[X] - min_XYZ_point_[X]);
 
        T_3_point::value_type  y = min_XYZ_point_[Y] 
                       + rand_between_0_and_1() * (max_XYZ_point_[Y] - min_XYZ_point_[Y]);
 
        T_3_point::value_type  z = min_XYZ_point_[Z] 
                       + rand_between_0_and_1() * (max_XYZ_point_[Z] - min_XYZ_point_[Z]);
 
        return T_3_point(x, y, z);
    }
};
//////////////////////////////////////////////////////////////////////////
class  T_is_v_objedinenii_sfer
{
    T_sfery  sfery_;    
public:
    T_is_v_objedinenii_sfer(T_sfery  sfery) : sfery_(sfery)
    {}
 
    bool operator() (T_3_point  point)
    {           
        for(T_sfery::const_iterator  cur_sfera_it = sfery_.begin();
            cur_sfera_it != sfery_.end(); ++cur_sfera_it)
        {
            if(cur_sfera_it->is_vnutri_sfery(point)) return true;
        }
        return false;
    }
};
//////////////////////////////////////////////////////////////////////////           
T_3_point::value_type get_objem_objedin_sfer
    (
        T_sfery  sfery,
        int      n
    )
{    
    //Находим минимальный параллелепипед, содержащий все заданные сферы.
    T_3_point  min_XYZ_point;
 
    for(T_coord_ind  cur_coord_ind = X; cur_coord_ind < TOTAL_COORDS; ++cur_coord_ind)
    {
        min_XYZ_point[cur_coord_ind] 
            = std::for_each(sfery.begin(), sfery.end(), T_get_min_coord(cur_coord_ind));    
    }
 
    T_3_point  max_XYZ_point;
 
    for(T_coord_ind  cur_coord_ind = X; cur_coord_ind < TOTAL_COORDS; ++cur_coord_ind)
    {
        max_XYZ_point[cur_coord_ind] 
            = std::for_each(sfery.begin(), sfery.end(), T_get_max_coord(cur_coord_ind));    
    }
 
    //Вектор разностей максимальных и минимальных координат.
    T_3_point  dliny_storon_parallelepipeda = max_XYZ_point - min_XYZ_point;
    
    T_3_point::value_type  objem_parallelepipeda = std::accumulate
                                        (
                                            dliny_storon_parallelepipeda.begin(),
                                            dliny_storon_parallelepipeda.end(),
                                            1.0,                                            
                                            std::multiplies<T_3_point::value_type>()
                                        );    
 
    //Генерируем вектор случайных точек.
    T_3_points  random_points(n);
 
    std::generate(random_points.begin(), random_points.end(), 
                  T_get_random_point(min_XYZ_point, max_XYZ_point));
 
    T_3_point::value_type  points_v_objedin_sfer_count 
        = std::count_if(random_points.begin(), random_points.end(), 
                        T_is_v_objedinenii_sfer(sfery));
 
    return objem_parallelepipeda * points_v_objedin_sfer_count / n;
}
//////////////////////////////////////////////////////////////////////////           
 
int main()
{
    std::locale::global(std::locale(""));
    srand(static_cast<unsigned>(time(0)));
 
    for(;;)
    {
        int sfer_total;
 
        std::cout << std::endl
                  << std::endl
                  << "Введите число сфер: ";
        std::cin >> sfer_total;
        if(sfer_total < 1) break;
        
        T_sfery  sfery(sfer_total);
        std::for_each(sfery.begin(), sfery.end(), 
                      std::mem_fun_ref(&T_sfera::input));
 
        for(;;)
        {
            std::cout << std::endl
                      << "Введите число испытаний : ";
            int n;
            std::cin >> n;
            if(n < 1) break;
            std::cout << "Ждите..."
                      << std::endl;
            std::cout << "Объем объединения "
                      << sfer_total
                      << " заданных сфер равен: "
                      << get_objem_objedin_sfer(sfery, n)
                      << std::endl;  
        }
    }  
    return 0;
}

#include "stdafx.h"
#include <stdlib.h>
#include <iostream>
#include <time.h>
#include <math.h>
using namespace std;
 
double f(double x);
 
int main()
{
    srand((unsigned)time(NULL));
    double x_rand,y_rand;
    int s=0;
    double PI=3.141592653589;
    int random_point=100000;
    for(int i=0;i<random_point;i++)
    {
        x_rand = (double)rand()/RAND_MAX*PI;
        y_rand = (double)rand()/RAND_MAX*(1+exp(PI/2));
        if(((x_rand<PI/2)&&(y_rand>0)&&(y_rand<f(x_rand)))||((x_rand>PI/2)&&(y_rand<0)&&(y_rand>f(x_rand)))) s+=1;
    }
    double integral=s/(double)random_point*(PI*(1+exp(PI/2)));
    cout << "S = "<< integral <<endl;
    system("Pause");
    return 0;
}
 
double f(double x)
{
    return exp(x)*pow(cos(x),3);
}