尤拉計劃695：隨機長方形

單位正方形內3個點，P1P2P3,作為對角的頂點，形成P1P2 P2P3 P1P3 3個長方形，求面積第二大的矩形面積的平均值。

沒找到數學公式,先用程式模擬一下

from random import random
def mid(m): 
 p=[random() for i in range(0,m*6)]
 s=0
 for i in range(0,m):
  x1,x2,x3,y1,y2,y3=p[i*6-6],p[i*6-5],p[i*6-4],p[i*6-3],p[i*6-2],p[i*6-1]
  s1=abs((x2-x1)*(y2-y1))
  s2=abs((x2-x3)*(y2-y3))
  s3=abs((x3-x1)*(y3-y1))
  s+=(s1+s2+s3-max(s1,s2,s3)-min(s1,s2,s3))
 return s/m  

def s695(n):
 m=10**6
 a=0
 for j in range(0,n//m):
  a+=mid(m)
 return a/(n//m)

結果不穩定

print(s695(10**7))
0.101732401264
print(s695(10**8))
0.101784964399
print(s695(10**9))
0.101781944286

再測試一下平均值,最大面積,最小面積,試圖看出規律

def avg(m): 
 p=[random() for i in range(0,m*6)]
 s=0
 for i in range(0,m):
  x1,x2,x3,y1,y2,y3=p[i*6-6],p[i*6-5],p[i*6-4],p[i*6-3],p[i*6-2],p[i*6-1]
  s1=abs((x2-x1)*(y2-y1))
  s2=abs((x2-x3)*(y2-y3))
  s3=abs((x3-x1)*(y3-y1))
  s+=(s1+s2+s3)/3.0
 return s/m  

def big(m): 
 p=[random() for i in range(0,m*6)]
 s=0
 for i in range(0,m):
  x1,x2,x3,y1,y2,y3=p[i*6-6],p[i*6-5],p[i*6-4],p[i*6-3],p[i*6-2],p[i*6-1]
  s1=abs((x2-x1)*(y2-y1))
  s2=abs((x2-x3)*(y2-y3))
  s3=abs((x3-x1)*(y3-y1))
  s+=max(s1,s2,s3)
 return s/m 

def lit(m): 
 p=[random() for i in range(0,m*6)]
 s=0
 for i in range(0,m):
  x1,x2,x3,y1,y2,y3=p[i*6-6],p[i*6-5],p[i*6-4],p[i*6-3],p[i*6-2],p[i*6-1]
  s1=abs((x2-x1)*(y2-y1))
  s2=abs((x2-x3)*(y2-y3))
  s3=abs((x3-x1)*(y3-y1))
  s+=min(s1,s2,s3)
 return s/m

print(avg(10**6))
0.11112487950547068
print(big(10**6))
0.20462836607020585
print(lit(10**6))
0.026774104377932796

把上述python程式改成c程式，時間可以縮短一半左右。（n=10^8，pypy 2.7 40s, gcc7.1 20s）

#include <cstdio>                     
#include <cstdlib>                    
#include <cmath>                      
#include <ctime>                      
#define max(x,y,z) ((x>=y?x:y)>=z?y:z)
#define min(x,y,z) ((x<=y?x:y)<=z?y:z)
double mid(int m,double p[])
{
 srand (time(NULL));
 for (int i=0; i<m*6; i++)
  p[i]=double(rand())/RAND_MAX;
 double s=0.0,s1,s2,s3;
 double x1,x2,x3,y1,y2,y3;
 for (int i=0; i<m; i++)
 {
  x1=p[i*6-6],x2=p[i*6-5],x3=p[i*6-4],y1=p[i*6-3],y2=p[i*6-2],y3=p[i*6-1];
  s1=fabs((x2-x1)*(y2-y1));
  s2=fabs((x2-x3)*(y2-y3));
  s3=fabs((x3-x1)*(y3-y1));
  s+=(s1+s2+s3-max(s1,s2,s3)-min(s1,s2,s3));
 }
 return s/m;
}
int main()
{
 int m=1000000,n=1000000*100;
 double a=0;
 static double p[1000000*6];
 for (int j=0; j<int(n/m); j++)
  a+=mid(m,p);
 printf("%lf",n/m));
}

另一種模擬方法，把正方形切割成更小的方塊，逐個點羅列計算，此法不執行隨機數操作，甚至更慢（m=30，7.29e8個組合用時6s），但這裡最後兩重迴圈是傻算，其實2點確定以後，第3點落在什麼範圍內，哪個長方形第二大是能夠確定的。

double mid_b(int m)//m<=100
{
 int s1,s2,s3;
 double s=0.0;
 for (int x1=0; x1<m; x1++)
  for (int x2=0; x2<m; x2++)
   for (int y1=0; y1<m; y1++)
    for (int y2=0; y2<m; y2++)
    {
     s1=abs((x2-x1)*(y2-y1));
     for (int x3=0; x3<m; x3++)
      for (int y3=0; y3<m; y3++)
      {
       s2=abs((x2-x3)*(y2-y3));
       s3=abs((x3-x1)*(y3-y1));
       s+=(s1+s2+s3-max(s1,s2,s3)-min(s1,s2,s3));
      }
    }
 return double(s)/pow(double(m),6.0+2);
}

如果長方形p1p2確定並設p2的座標大於p1，並且讓它成為面積第2的那個，那麼如果p3在長方形p2頂點的右下角,雙曲線x*y=s（p1p2面積）上方的點都滿足，同理左上方的雙曲線下方的也都滿足。這個面積的平均值就變成雙曲線與“座標軸”圍起面積的比例。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
def g(a,b,a1,b1):
 x=[]
 y=[]
 for i in range(0,100):
  for j in range(0,100):
   s=abs((i-a)*(j-b))
   s1=abs((i-a1)*(j-b1))
   if s1>=abs((a-a1)*(b-b1)) and s<=abs((a-a1)*(b-b1)) or s1<=abs((a-a1)*(b-b1)) and s1>=abs((a-a1)*(b-b1)):
    x.append(i)
    y.append(j)
 fig = plt.figure()
 #plt.axis('scaled')
 plt.gca().set_aspect('equal', adjustable='box')
 ax = plt.subplot(111)
 ax.add_patch(patches.Rectangle((a, b), a1-a, b1-b))
 plt.scatter(x,y,c = 'r',marker = '.')
 plt.xlabel('X')
 plt.ylabel('Y')
 plt.show()                                                                                                 

呼叫上述python畫圖函式，可以得到，在100x100的方格內，p1p2為(30,30,60,60)時，滿足條件p3點的範圍。

p1p2為(30,30,40,40)時，滿足條件p3點的範圍。

p1p2為(30,40,60,50)時，滿足條件p3點的範圍。