最近我发现了一个有趣的问题,问题如上。问题简化一下,大概就是一圆内有四个随机生成的点,求四个点在同一半圆的概率。
问题分析:
- 怎么判断四点是否在同一半圆内?
- 怎么判断一个点是否在一个凸多边形内?
- 怎么通过四个点生成一个凸多边形?
- 怎么定量判断一点是否在一个三角形内?
先解决第一个问题,那我们都知道,任意一天圆的直径都能将一个圆平分,而任意一条直径都过圆心,所以只要圆心不在四点所围成的凸多边形,那为什么是凸多边形呢?后面会有解答。
解决第一个问题后,随之衍生出第二个问题,先明确凸多边形是啥意思?凸多边形(Convex Polygon)指如果把一个多边形的所有边中,任意一条边向两方无限延长成为一直线时,其他各边都在此直线的同旁,那么这个多边形就叫做凸多边形。简单的说,就是内角都得小于180才叫凸多边形。
再看一什么是非凸多边形
可以看到,图中圆内四点所形成的就是一个非凸多边形,此时圆心确实也在该多边形的外面,但是却不能找到一条直径让这四点分布在同一侧,所以我们强调要凸多边形。
解决了上述问题,那就要思考一下给四个点怎么能生成一个凸多边形?我就想到,任一多边形都能由几个三角形拼接而成。所以,任给四个点,我们可以先生成C43也就是4个三角形,然后把这四个三角形围成的区域并起来就可以形成一个凸多边形,当然后面也证实这确实是我们所需要的凸多边形。类似下面这样。
上图中,圆内四点ABCD,圆心O,然后是有ABC,ABD,ACD,BCD四个三角形,可以看到这四个三角形区域并起来围成的区域是三角形ABD,圆心位于该多边形内(三角形也是凸多边形),这时就得到了正确的结论。
到现在问题已经转化到如何判断一点是否在一个三角形内了。因为,圆内生成的四点可以三三生成四个三角形,那么我们只需要保证圆心不在任意一个生成的三角形内就可以保证该四点可被一条直径划分到同一侧,即全部分布在一个半圆内。当然相反地,圆心只要分布在至少一个三角形内就说明没有一条直径能够将这四点划分在同一侧,即非全部分布在一个半圆内。那么到底该怎么判断勒?看下图
上图呢,是一个任意的三角形ABC和内外两点P1,P2,将P1,P2都与三角形的三个顶点相连又可以得到三个三角形,观察得到,内点P1形成的三个三角形面积之和就等于ABC的面积,而外点P2形成的三个三角形面积之和比ABC的面积大一个P2BC的面积。所以得出结论就是只需要计算点与三角形顶点所形成的三个三角形面积之和再与原三角形面积比较即可得到该点是否在三角形内。
写到这里,如果不用概率论与数理统计的知识去解决,而单纯的用计算机模拟的话,问题已经变得很简单了。继问题四往后推,也就是需要实现一个求三角形面积的函数,实现如下
double Striangle(double *A,double *B,double *C)
{
//输入三点ABC坐标(二维数组)
//输出所围三角形面积
double deta_x1,deta_x2,deta_y1,deta_y2;
double Ax,Ay,Bx,By,Cx,Cy;
double len_a,len_b,sineC;
double S=-1;
Ax=A[0];
Ay=A[1];
Bx=B[0];
By=B[1];
Cx=C[0];
Cy=C[1];
deta_x1=(Ax-Bx);
deta_x2=(Ax-Cx);
deta_y1=(Ay-By);
deta_y2=(Ay-Cy);
if(deta_x1*deta_x2!=0 && deta_y1/deta_x1==deta_y2/deta_x2) {
return 0;
} else if(deta_x1==0 && deta_y1==0) {
return 0;
} else if(deta_x2==0 && deta_y2==0) {
return 0;
} else {
len_a=length(B,C);
len_b=length(A,C);
sineC=sine_angle(A,B,C);
S=0.5*len_a*len_b*sineC;
return S;
}
}
可以看到在运行这个函数之前,还需要再实现两个函数求两点距离和三角形某个角正弦值,实现如下
double length(double *A,double *B)
{
//输入两点AB坐标(二维数组)
//输出两点距离(模长)
double Ax,Ay,Bx,By,R;
Ax=A[0];
Ay=A[1];
Bx=B[0];
By=B[1];
R=sqrt(pow(Ax-Bx,2) + pow(Ay-By,2));
return R;
}
double sine_angle(double *A,double *B,double *C)
{
//输入三点ABC坐标(二维数组)
//输出角C的正弦值
double len_a,len_b,len_c;
double cosine2ab,cosine,sine;
len_a=length(B,C);
len_b=length(A,C);
len_c=length(A,B);
cosine2ab=pow(len_a,2)+pow(len_b,2)-pow(len_c,2);
if(len_a*len_b!=0) {
cosine=cosine2ab / (2*len_a*len_b);
} else {
cosine=1;
}
sine=sqrt(1-pow(cosine,2));
return sine;
}
到这里就没啥难度了,下面再说一下怎么随机生成四点坐标,实现如下
void rand_circle_coord(*CC)
{
//输入一个空二维数组
//将随机数赋值给传入的指针
//srand(time(0));
CC[0]=rand()%(2*eps+1) - eps;
CC[1]=rand()%(2*eps+1) - eps;
}
void rand_circle_coords(double *A,double *B,double *C,double *D)
{
//输入四个二维数组
//将随机坐标复制给传入的四个指针
int intA[2],intB[2],intC[2],intD[2];
double O[2]= {0,0};
//srand(time(0));
rand_circle_coord(intA);
A[0]=intA[0] / (double)eps;
A[1]=intA[1] / (double)eps;
while(length(A,O)>=1.00) {
rand_circle_coord(intA);
A[0]=intA[0] / (double)eps;
A[1]=intA[1] / (double)eps;
}
//printf("A是(%lf,%lf)\n",A[0],A[1]);
rand_circle_coord(intB);
B[0]=intB[0] / (double)eps;
B[1]=intB[1] / (double)eps;
while(length(A,O)>=1.00) {
rand_circle_coord(intB);
B[0]=intB[0] / (double)eps;
B[1]=intB[1] / (double)eps;
}
//printf("B是(%lf,%lf)\n",B[0],B[1]);
rand_circle_coord(intC);
C[0]=intC[0] / (double)eps;
C[1]=intC[1] / (double)eps;
while(length(C,O)>=1.00) {
rand_circle_coord(intC);
C[0]=intC[0] / (double)eps;
C[1]=intC[1] / (double)eps;
}
//printf("C是(%lf,%lf)\n",C[0],C[1]);
rand_circle_coord(intD);
D[0]=intD[0] / (double)eps;
D[1]=intD[1] / (double)eps;
while(length(D,O)>=1.00) {
rand_circle_coord(intD);
D[0]=intD[0] / (double)eps;
D[1]=intD[1] / (double)eps;
}
//printf("D是(%lf,%lf)\n",D[0],D[1]);
}
到这里基本的一些函数就全部实现好了,剩下的就是将这些函数联合起来用,达到模拟的效果,整个代码如下
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
#define eps 1000000
double length(double *A,double *B);
double sine_angle(double *A,double *B,double *C);
double Striangle(double *A,double *B,double *C);
int judgeP_in_triangle(double *A,double *B,double *C,double *P);
int judgeO_in_triangles(double *A,double *B,double *C,double *D);
void rand_circle_coord(int *CC);
void rand_circle_coords(double *A,double *B,double *C,double *D);
int main()
{
double A[2];
double B[2];
double C[2];
double D[2];
int i,j,P=0,R[100],epochs=1000;
srand(time(0));
for(j=0; j<100; j++) {
R[j]=0;
for(i=0; i<epochs; i++) {
rand_circle_coords(A,B,C,D);
if(judgeO_in_triangles(A,B,C,D)) {
R[j]+=1;
}
}
}
for(j=0;j<100;j++){
P+=R[j];
}
printf("\n结果是%lf",(double)P / (epochs*100));
return 0;
}
double length(double *A,double *B)
{
//输入两点AB坐标(二维数组)
//输出两点距离(模长)
double Ax,Ay,Bx,By,R;
Ax=A[0];
Ay=A[1];
Bx=B[0];
By=B[1];
R=sqrt(pow(Ax-Bx,2) + pow(Ay-By,2));
return R;
}
double sine_angle(double *A,double *B,double *C)
{
//输入三点ABC坐标(二维数组)
//输出角C的正弦值
double len_a,len_b,len_c;
double cosine2ab,cosine,sine;
len_a=length(B,C);
len_b=length(A,C);
len_c=length(A,B);
cosine2ab=pow(len_a,2)+pow(len_b,2)-pow(len_c,2);
if(len_a*len_b!=0) {
cosine=cosine2ab / (2*len_a*len_b);
} else {
cosine=1;
}
sine=sqrt(1-pow(cosine,2));
return sine;
}
double Striangle(double *A,double *B,double *C)
{
//输入三点ABC坐标(二维数组)
//输出所围三角形面积
double deta_x1,deta_x2,deta_y1,deta_y2;
double Ax,Ay,Bx,By,Cx,Cy;
double len_a,len_b,sineC;
double S=-1;
Ax=A[0];
Ay=A[1];
Bx=B[0];
By=B[1];
Cx=C[0];
Cy=C[1];
deta_x1=(Ax-Bx);
deta_x2=(Ax-Cx);
deta_y1=(Ay-By);
deta_y2=(Ay-Cy);
if(deta_x1*deta_x2!=0 && deta_y1/deta_x1==deta_y2/deta_x2) {
return 0;
} else if(deta_x1==0 && deta_y1==0) {
return 0;
} else if(deta_x2==0 && deta_y2==0) {
return 0;
} else {
len_a=length(B,C);
len_b=length(A,C);
sineC=sine_angle(A,B,C);
S=0.5*len_a*len_b*sineC;
return S;
}
}
int judgeP_in_triangle(double *A,double *B,double *C,double *P)
{
int flag;
double Spab,Spac,Spbc,Sabc,deta_S;
Spab=Striangle(A,B,P);
Spac=Striangle(A,C,P);
Spbc=Striangle(B,C,P);
Sabc=Striangle(A,B,C);
deta_S=Spab+Spac+Spbc-Sabc;
if(deta_S<=1e-06) {
flag=1;
} else {
flag=0;
}
return flag;
}
int judgeO_in_triangles(double *A,double *B,double *C,double *D)
{
double O[2]= {0,0};
int flag1,flag2,flag3,flag4;
flag1=!judgeP_in_triangle(A,B,C,O);
flag2=!judgeP_in_triangle(A,B,D,O);
flag3=!judgeP_in_triangle(A,C,D,O);
flag4=!judgeP_in_triangle(B,C,D,O);
if(flag1 && flag2 && flag3 && flag4) {
return 1;
} else {
return 0;
}
}
void rand_circle_coord(int *CC)
{
//输入一个空二维数组
//将随机数赋值给传入的指针
//srand(time(0));
CC[0]=rand()%(2*eps+1) - eps;
CC[1]=rand()%(2*eps+1) - eps;
}
void rand_circle_coords(double *A,double *B,double *C,double *D)
{
//输入四个二维数组
//将随机坐标复制给传入的四个指针
int intA[2],intB[2],intC[2],intD[2];
double O[2]= {0,0};
//srand(time(0));
rand_circle_coord(intA);
A[0]=intA[0] / (double)eps;
A[1]=intA[1] / (double)eps;
while(length(A,O)>=1.00) {
rand_circle_coord(intA);
A[0]=intA[0] / (double)eps;
A[1]=intA[1] / (double)eps;
}
//printf("A是(%lf,%lf)\n",A[0],A[1]);
rand_circle_coord(intB);
B[0]=intB[0] / (double)eps;
B[1]=intB[1] / (double)eps;
while(length(A,O)>=1.00) {
rand_circle_coord(intB);
B[0]=intB[0] / (double)eps;
B[1]=intB[1] / (double)eps;
}
//printf("B是(%lf,%lf)\n",B[0],B[1]);
rand_circle_coord(intC);
C[0]=intC[0] / (double)eps;
C[1]=intC[1] / (double)eps;
while(length(C,O)>=1.00) {
rand_circle_coord(intC);
C[0]=intC[0] / (double)eps;
C[1]=intC[1] / (double)eps;
}
//printf("C是(%lf,%lf)\n",C[0],C[1]);
rand_circle_coord(intD);
D[0]=intD[0] / (double)eps;
D[1]=intD[1] / (double)eps;
while(length(D,O)>=1.00) {
rand_circle_coord(intD);
D[0]=intD[0] / (double)eps;
D[1]=intD[1] / (double)eps;
}
//printf("D是(%lf,%lf)\n",D[0],D[1]);
}
这里我们模拟1000次,重复实验100次然后取平均结果如下
大概是0.5的概率,一半一半吧。
好了。整个问题我的想法就是这样了。如有错误欢迎指正,一起进步,加油💪💪💪