def midpoint_double1(f, a, b, c, d, nx, ny):
	hx = (b-a)/float(nx)
	hy = (d-c)/float(ny)
	I = 0
	for i in range(nx):
		for j in range(ny):
			xi = a + hx/2 + i*hx
			yj = c + hy/2 + j*hy
			I += hx*hy*f(xi,yj)
	return I

def midpoint(f, a, b, n):
	h = float(b-a)/n
	result = 0
	for i in range(n):
		result += f((a + h/2.0) + i*h)
	result *= h
	return result

def midpoint_double2(f, a, b, c, d, nx, ny):
	def g(x):
		return midpoint(lambda y: f(x,y), c, d, ny)

	return midopoint(g, a, b, nx)

def test_midpoint_double():
	"""Test that a linear function is integrated exactly."""
	def f(x,y):
		return 2*x + y

	a = 0; b = 2; c = 2; d = 3
	import sympy
	x, y = sympy.symbols('x y')
	I_expected = sympy.integrate(f(x,y), (x,a,b), (y,c,d))
	# Test three cases: nx < ny, nx = ny, nx > ny
	for nx, ny in (3,5), (4,4), (5,3):
		I_computed1 = midpoint_double1(f, a, b ,c ,d, nx, ny)
		I_computed2 = midpoint_double2(f, a, b, c, d, nx, ny)
		tol = 1E-14
		#print I_expected, I_computed1, I_computed2
		assert abs(I_computed1 - I_expected) < tol
		assert abs(I_computed2 - I_expected) < tol
