public
anonymous / integrate.py
Created

For Reddit Math Problem

  • Download Gist
integrate.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
import random;
 
def magsq(v):
return sum(x**2 for x in v)
 
def sub(a,b):
return (x[0]-x[1] for x in zip(a,b))
 
def in_sph(r0,r,pt):
return magsq(sub(pt,r0)) <= r**2
 
def rand_pt(domain):
return (random.uniform(domain[0], domain[1]),
random.uniform(domain[2], domain[3]),
random.uniform(domain[4], domain[5]))
 
def dom_vol(domain):
return (domain[1]-domain[0])*(domain[3]-domain[2])*(domain[5]-domain[4])
 
def inside(pt):
return in_sph((0,0,0),1.0,pt) or \
in_sph((1,0,0),1.0,pt) or \
in_sph((0,1,0),1.0,pt) or \
in_sph((1,1,0),1.0,pt) or \
in_sph((0,0,1),1.0,pt) or \
in_sph((1,0,1),1.0,pt) or \
in_sph((0,1,1),1.0,pt) or \
in_sph((1,1,1),1.0,pt)
 
Npts = 0
Ninside = 0
 
D = (-1, 2, -1, 2, -1, 2)
 
for N in xrange(1000000):
if inside(rand_pt(D)): Ninside += 1
Npts += 1
 
print 1.0*dom_vol(D)*Ninside/Npts

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.