126 lines
3.0 KiB
Python
Executable File
126 lines
3.0 KiB
Python
Executable File
import time
|
|
import random
|
|
import math
|
|
from multiprocessing import Process, Value, Lock, Pool
|
|
|
|
def main():
|
|
start = time.time()
|
|
# number process
|
|
k=4
|
|
# points per process
|
|
n=250000000
|
|
|
|
unlinked=0
|
|
linked=0
|
|
|
|
pool = Pool(k)
|
|
|
|
results=[pool.apply_async(generate_data,args=(n,)) for i in range(k)]
|
|
|
|
for result in results:
|
|
result.wait()
|
|
output = result.get()
|
|
#print(output)
|
|
unlinked+=output[0]
|
|
linked+=output[1]
|
|
|
|
end = time.time()
|
|
print('Linked triangles: ' + str(linked))
|
|
print('Unlinked triangles: ' + str(unlinked))
|
|
print('Total time: ' + str(end-start))
|
|
p=1.0*linked/(n*k)
|
|
print('99% confidence interval for q: ' + str(p*4/18) + ' +/- ' +
|
|
str(2.575*4/18*math.sqrt(p*(1-p)/n)))
|
|
|
|
def generate_data(n):
|
|
unlinked=0
|
|
linked=0
|
|
for i in range(n):
|
|
if random_triangles()==0:
|
|
unlinked+=1
|
|
else:
|
|
linked+=1
|
|
return [unlinked,linked]
|
|
|
|
|
|
def generate_coords(len=6):
|
|
pts=[]
|
|
for i in range(len):
|
|
pts.append(random.random())
|
|
return pts
|
|
|
|
def random_triangles():
|
|
x=generate_coords()
|
|
y=generate_coords()
|
|
z=generate_coords()
|
|
#mathematica code to plot triangles
|
|
# print('Graphics3D[{Red,Line[{{'+str(x[0])+','+str(y[0])+','+str(z[0])+'},{'
|
|
# +str(x[1])+','+str(y[1])+','+str(z[1])+'},{'
|
|
# +str(x[2])+','+str(y[2])+','+str(z[2])+'},{'
|
|
# +str(x[0])+','+str(y[0])+','+str(z[0])+'}}],Blue,Line[{{'
|
|
# +str(x[3])+','+str(y[3])+','+str(z[3])+'},{'
|
|
# +str(x[4])+','+str(y[4])+','+str(z[4])+'},{'
|
|
# +str(x[5])+','+str(y[5])+','+str(z[5])+'},{'
|
|
# +str(x[3])+','+str(y[3])+','+str(z[3])+'}}]}]')
|
|
|
|
linknum = 0
|
|
for i in range(3):
|
|
for k in range(3):
|
|
j= (i+1) % 3
|
|
m=k+3
|
|
n= 3+( (k+1) % 3 )
|
|
linknum+=check_crossing([x[i],y[i],z[i]],[x[j],y[j],z[j]],
|
|
[x[m],y[m],z[m]],[x[n],y[n],z[n]])
|
|
# print(str(crossings))
|
|
return linknum
|
|
|
|
|
|
def check_crossing(p1,p2,q1,q2):
|
|
# solve linear system for when two lines cross
|
|
dxp = p2[0]-p1[0]
|
|
dyp = p2[1]-p1[1]
|
|
dxq = q2[0]-q1[0]
|
|
dyq = q2[1]-q1[1]
|
|
dxc = q1[0]-p1[0]
|
|
dyc = q1[1]-p1[1]
|
|
det=dyp*dxq-dxp*dyq
|
|
# det=(p2[0]-p1[0])*(q1[1]-q2[1])-(q1[0]-q2[0])*(p2[1]-p1[1])
|
|
t=( dxq*dyc-dyq*dxc ) / det
|
|
# t=( (q1[1]-q2[1])*(q1[0]-p1[0])+(q2[0]-q1[0])*(q1[1]-p1[1]) ) / (
|
|
# det )
|
|
s=( dxp*dyc-dyp*dxc ) / det
|
|
# s=( (p1[1]-p2[1])*(q1[0]-p1[0])+(p2[0]-p1[0])*(q1[1]-p1[1]) ) / (
|
|
# det )
|
|
|
|
# check if the lines cross inside the segments from p1 to p2 and q1 to q2
|
|
if (0<t<1) and (0<s<1):
|
|
# segments cross, find which is the over crossing
|
|
if (p1[2] + t*(p2[2]-p1[2])) > (q1[2]+s*(q2[2]-q1[2])):
|
|
overunder= 1
|
|
else:
|
|
overunder= -1
|
|
else:
|
|
overunder= 0
|
|
|
|
# also check the possibility of numerical error
|
|
#if abs((p1[2] + t*(p2[2]-p1[2])) - (q1[2]+s*(q2[2]-q1[2])))<0.0000000001:
|
|
# print('possible rounding error')
|
|
|
|
# now determine sign of crossing
|
|
if (det*overunder >0):
|
|
return 1
|
|
elif (det*overunder < 0):
|
|
return -1
|
|
else:
|
|
return 0
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|
|
|
|
# Linked triangles: 152402780
|
|
# Unlinked triangles: 847597220
|
|
# Total time: 11496.435719013214
|
|
# 99% confidence interval for q: 0.033867284444444444 +/- 1.300726300256763e-05
|
|
|