0
\$\begingroup\$

I have been working on a code to distribute randomly generated points uniformly in a circle of radius 1 centered at origin. I have tried two variations of the code, yet they both yield one major problem.

Even for as many as 500 points, there is a particular region (quadrant 1, (0,1) on X Axis, and about (0,0.5) on Y axis) which is literally empty and it seems like that is a gaping hole in an otherwise perfectly well distributed points in the circle (the points seems to be distributed uniformly in the region outside this). The following are the two codes I tried out:

#include<iostream>
#include<random>
#include<time.h>
#include<fstream>
#define pi atan(1)*4
using namespace std;
int main()
{
    int N;  cin>>N;
    double X,Y,x,y,a,b;
    srand(time(NULL));
    for(int i=0; i<N; i++)
    {
        do
            {   X=(double)rand()/(double)RAND_MAX;    
                Y=(double)rand()/(double)RAND_MAX;
                x=sqrt(X);    
                y = 2*(pi)*Y;       
                a=x*cos(y);   
                b=x*sin(y);
            }while(((a*a)+(b*b))<=1);
    }
    return 0;
}

The second code is:

#include<iostream>
#include<random>
#include<time.h>
#include<fstream>
#define pi atan(1)*4
using namespace std;
int main()
{
    int N;  cin>>N;
    double X,Y,x,y,a,b;
    random_device R;    
    mt19937 G(R());     
    uniform_real_distribution<double> D(0,1);
    for(int i=0; i<N; i++)
    {
        do
            {   
                X=D(G);     Y=D(G);     //uniform real distribution in (0,1)
                x=sqrt(X);    
                y = 2*(pi)*Y;        
                a=x*cos(y);   
                b=x*sin(y);
            }while(((a*a)+(b*b))<=1);

    }
    return 0;
}

What is the possible reason for the problem I'm facing?

Update: For some reason, all points started getting generated at the boundary if I use ((aa)+(bb)<1) in the do-while loops. The problem I mentioned earlier has ((x * x)+(y * y)<1) which is clearly wrong, my bad. I replaced the do-while of ((a * a)+(b * b)<1) with a goto as follows

A:
X=D(G);     Y=D(G);     //uniform real distribution in (0,1)
x=sqrt(X);    
y = 2*(pi)*Y;        
a=x*cos(y);   
b=x*sin(y);
if((a*a)+(b*b)>1)
{  gota A:   }

and now it seems to be working fine. What could be the reason that goto makes it work but do-while does not?

\$\endgroup\$

closed as off-topic by JS1, Graipher, alecxe, Edward, pacmaninbw Jun 9 '17 at 13:39

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "Questions containing broken code or asking for advice about code not yet written are off-topic, as the code is not ready for review. After the question has been edited to contain working code, we will consider reopening it." – JS1, Graipher, alecxe, Edward, pacmaninbw
If this question can be reworded to fit the rules in the help center, please edit the question.

  • \$\begingroup\$ Neither program produces any output, so it's not clear what you expect to see (or whether that is actually what's produced). Are you intending to produce points on a unit circle (i.e. r==1) or in a unit circle (r < 1)? \$\endgroup\$ – Toby Speight Jun 12 '17 at 8:18
3
\$\begingroup\$

You can paramterize a circle via a radius and the angle. Therefore, for the cycle with radius r around (h,w) you can calculate the x and y coordinates with the following equations

x = h + r * sin(theta)

y = w + r *cos(theta)

Therefore your code should look like that

random_device R;    
mt19937 G(R());   
uniform_real_distribution<double> D(0,360);

for(int i=0; i<N; i++) {
    angle = D(G);
    x = sin(angle);
    y = cos(angle);
}

UPDATE

Toby is right, that this only adds points to the unit circle, rather than the full area of the circle. The solution is rather simple, add a random variable for the radius

unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
mt19937 rngAngle(seed);  
seed = std::chrono::system_clock::now().time_since_epoch().count(); 
mt19937 rngRadius(seed);   
uniform_real_distribution<double> angleDistribution(0,360);
uniform_real_distribution<double> radiusDistribution(0,1);

for(int i=0; i<N; i++) {
    angle  = angleDistribution(rngAngle);
    radius = radiusDistribution(rngRadius);
    x = radius * sin(angle);
    y = radius * cos(angle);
}

As you have seen, I used some better names too and seeded both random number generators independently with the current time.

\$\endgroup\$
  • \$\begingroup\$ Whilst that's true, it appears that the code in the question is generating points in the unit disc: (r, theta) = (sqrt(rnd(1,1)), rnd(0, 2π)); That's not going to be uniform in the area, though - you need something like (1/rnd(1, ∞), rnd(0, 2π)) for that. \$\endgroup\$ – Toby Speight Jun 9 '17 at 13:01
  • \$\begingroup\$ Thats is exactly what he meant? However, the solution is rather easy, as you only need to add an second random variable that goes from 0 to radius and multiply the coordinates with it. \$\endgroup\$ – miscco Jun 9 '17 at 14:23
  • \$\begingroup\$ The original code does have two random variables, but it's not clear what was intended (especially with that "while hypot <= 1" loop). But the "r" variable there is near-uniformly distributed (which would be wrong if you want uniformly-distributed points). \$\endgroup\$ – Toby Speight Jun 9 '17 at 14:33

Not the answer you're looking for? Browse other questions tagged or ask your own question.