Metoda Monte Carlo całki podwójne

Witam,
Mam nadzieję, że zakładam topic w dobrym dziale ; )

Otóż mój problem jest taki :

Mam taki oto algorytm liczenia całek podwójnych metodą monte carlo :

#include <iostream>
#include <math.h>
using namespace std;


double f(double x, double y)
{

    return 1/((y+x+1)*(y+x+1));

}

int main()
{
	int N = 5000000;
	int i;
	double s,xp,xk,dx,dy,yp,yk;
	cout << "Podaj poczatek przedzialu x: " << endl;
	cin >> xp;
	cout << "Podaj koniec przedzialu x: " << endl;
	cin >> xk;
    cout << "Podaj poczatek przedzialu y: " << endl;
	cin >> yp;
	cout << "Podaj koniec przedzialu y: " << endl;
	cin >> yk;
	srand((unsigned)time(NULL));
	dx = (xk-xp);
	dy = (yk-yp);
	for (int j = 0; j<200; j++){
	s = 0;
	for (i=1;i<=N;i++) {
		s += f(xp+((double)rand()/(double)(RAND_MAX)*dx),yp+((double)rand()/(double)(RAND_MAX)*dy));
	}
	s = (dx*dy / N) *s;
	cout << "Wartosc calki wynosi : " << s << endl;
	}
	system("pause");
}

To o dziwo działa, ale teraz za nic nie mogę wymyśleć, jak to zaadaptować, kiedy obszar nie jest prostokątem tzn :
np.
0 < x < 2 x/2 < y < 2x

Generalnie z tego co rozumiem z tej metody, to w bloku gdzie liczę s+=... itd muszę dodać jakiś if, który sprawdza czy wylosowany punkt (x,y) należy do obszaru, jeśli nie, to niczego nie dodaje, jeśli tak, algorytm się wykonuje w tej postaci s+=... itd.

Miał już ktoś taki problem kiedyś ? : D

pozdro ; ) i z góry dzięki

Po pierwsze musisz zdecydować, czy N to jest ilość losowań, czy ilość wylosowanych punktów spełniających warunek. Poniżej zakładam, że N jest ilością losowań.
int ileDobrych=0;

        for (i=1;i<=N;i++) 
       {
              double kandydatX=((double)rand()/(double)(RAND_MAX)*dx);
              double kandydatY=((double)rand()/(double)(RAND_MAX)*dy);
              i//f(kandydatX i KandydatY leżą gdzie trzeba)
              if(kandydatY>=kandydatX/2 && kandydatY<=2*kandydatX)
              {
                    ileDobrych++;
                    s += f(xp+kandydatX,yp+kandydatY);
              }
        }
        s = poleObszaruPoKtorymCałkujesz*(s/ ileDobrych) ; //pole * średnia wartość funkcji

Jeżeli obszar, po którym całkujesz jest skomplikowany, to do obliczenia jego pola można też użyć metody Monte Carlo: pole = (dxdy)(ileDobrych/N).

Dzięki za odpowiedź ; ) właściwie to nie działało mi dlatego, że w ifie wpisałem jeszcze warunki na x, w Twoim przypadku kandydatX, a nie są one potrzebne, bo z założenia są w obszarze. No zależy względem której osi obszar jest normalny. Jeśli chodzi o dalszą część to mam i ze swojej strony dodam, że
s = dx*dy * (s/N) jest najdokładniejszym sposobem mi znanym i sprawdza się ( potwierdziłem dużą ilością przykładów ; ) i porównań do innych rodzajów tego algorytmu )

pozdro i dzięki