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