Объем тела методом Монте-Карло в С++

В прикладных расчетах методы Монте-Карло применяются достаточно часто. Причем нередко соответствующие методы исследования являются единственно возможными. В данном случае не будем останавливаться на особенностях методов Монте-Карло, ограничимся лишь некоторыми практическими рецептами их использования, в частности, при вычислении объемов (площадей) тел. Предположим, необходимо вычислить объем тела (области пространства), форма которой достаточно сложна для того, чтобы это сделать аналитически. В этом случае тело, объем которого нужно вычислить, условно и целиком помещается внутрь другого тела, объем которого известен или легко вычисляется - например, квадрат или сферу. Далее случайным образом внутри этого внешнего тела выбираются точки и вычисляется отношение' точек, которые попали внутрь меньшего тела (объем которого вычисляется) к общему числу точек. Если общее число точек стремится к бесконечности, то указанное отношение стремится к отношению объемов внутреннего и внешнего тел. На практике вместо случайного выбора точек используют сетку из равномерно распределенных точек. В этом случае говорят о квази-методах Монте-Карло. В примере ниже вычисляется объем усеченного конуса, накрытого полусферой. Объем такого тела вычисляется по формуле: $$V=\frac{1}{3}\pi R^{2}H+\frac{2}{3}\pi R^{3}$$ где R-радиус, а H-высота усеченного конуса.

#include 
#include 
using namespace std;
int main(){
//Число pi:
const double pi=3.1415;
//Рабочие переменные программы:
double R,H,V,V0,x,y,z;
//Число точек N (по каждой из координат) 
//и число внутренних точек n:
int N=1500,n=0;
//Ввод параметров R и H:
cout << "Enter R = ";
cin >> R;
cout << "Enter H = ";
cin >> H;
//Объем параллелепипеда:
V0=4*R*R*(H+R);
//Перебор всех точек:
for(int i=0;i <= N;i++){
 x=2*i*R/N-R;
 for(int j=0;j <= N;j++){
 y=2*j*R/N-R;
 for(int k=0;k <= N;k++){
 z=k*(H+R)/N;
 //Подсчет внутренних точек:
 if(((sqrt(x*x+y*y)/R <= z/H)&&(z <= H))||((x*x+y*y+(z-H)
*(z-H) <= R*R)&&(z > H)))
 n++;
 }
 }
}
//Объем тела:
V=V0*n/pow(N+1,3);
//Вывод вычисленного и точного значений:
cout << "V = " << V << " : " << pi*R*R*H/3+2*pi*pow(R,3)/3 << endl;
return 0;
}