Смоделировать алгоритм непрерывной случайной величины распределенной с плотностью

//Ãåíåðàòîð ñëó÷àéíûõ ÷èñåë ñ çàäàííîé ïëîòíîñòüþ ðàñïðåäåëåíèÿ
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iomanip>
 
using namespace std;
 
double p(double x){//ïëîòíîñòü ðàñïðåäåëåíèÿ âåðîÿòíîñòè
    
    const double C = 2;
    const double MI_PI = 3.1415926535;
    
    return C*exp(-MI_PI*x*x);
}
 
double Max(double a, double b, int n){//íàèáîëüøåå çíà÷åíèå íà çàäàííîé ôóíêöèè
    
    double dx = (b - a)/(n - 1);
    double max = p(a);
    for (double x = a + dx; x < b + dx/2; x += dx){
        double m = p(x);
        if (m > max)
            max = m;
    }
    return max;
}
 
double MyDRand(double a, double b, double eps){//ðàâíîìåðíîå ðàñïðåäåëåíèå
    
    long long A = a/eps;
    long long B = b/eps;
    long long C = 1;
    
    while ((B - A + 1) > RAND_MAX){
        
        A /= 10;
        B /= 10;
        C *= 10;
    }
    return (rand()%(B - A + 1) + A)*(C*eps);
}
 
double Randome(double a, double b, double eps, int n){//ðàñïð. ñ çàäàííîé ïëîòí.
    
    double fmax = Max(a, b, n);
    
    double x, y;
    do{
         x = MyDRand(a, b, eps);
         y = MyDRand(0, 1, eps);
    }while (p(x)/fmax < y);
    
    return x;
}
 
int main(){
    
    srand(time(0));
    
    const double EPS = 1e-3;
    const int N = 1000;
    
    double a, b;
    cout << "a = "; cin >> a;
    do{
        cout << "b = "; 
        cin >> b; 
    }while (b <= a);
    
    const int M = 20;
    int out[M] = {};
    const int K = 1000;
    for (int i = 0; i < K; ++i){
        
        double r = Randome(a, b, EPS, N);
        int p = (r - a)*(M - 1)/(b - a);
        ++out[p];
    }
    
    cout << "K = " << K << endl;
    
    for (int i = 0; i < M; ++i)
        cout << setw(5) << fixed << setprecision(2) << a + i*(b - a)/M << " - " 
             << setw(5) << fixed << setprecision(2) << a + (i + 1)*(b - a)/M
             << " : " 
             << setw(5) << fixed << setprecision(3) << out[i]/(K + 0.) << endl;
}