For Reddit Math Problem

integrate.py

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