def factorial(n,r):
	if n == 0:
		return r
	else:
		return n*factorial(n-1,r)

def taylor(x,A,F,y0,k):
	Y=y0
	for i in range(k):
		Y = Y + x**(i+1)
	Y = Y + x**(k+2+1)*A
	return Y

def step(F,k,y0, eps, p):
	alternate = 1
	if (eps/(k+1.) < 1.): #10-38 / 10-42
		h = (eps/(k+1.))**(1./(k+1.))
	else:
		h = (eps/(k+1.))**(1./(k)) #oops
	A = [y0[0]+0,y0[1]+h] #10-39
	Y = [taylor(0,A[0],F,y0[0],k), taylor(h,A[1],F,y0[1],k)] #10-40
	while (Y[0]<A[0] or Y[1]>A[1]):
		print "bad", "A", A,"Y", Y, "h", h
		#bad luck, Y is not contained in A
		if alternate == 1:
			print "modding A"
			A = [taylor(0,A[0],F,y0[0],k), taylor(h,A[1],F,y0[1],k)] #10-40
			Y = [taylor(0,A[0],F,y0[0],k), taylor(h,A[1],F,y0[1],k)] #10-40
		else:
			print "halving h"
			h = h/2. 
			Y = [taylor(0,A[0],F,y0[0],k), taylor(h,A[1],F,y0[1],k)] #10-40
		alternate = 3 - alternate #swap alternative
	for point in p:
		if point >= 0 and point <= h:
			print p
			print Y(p)
	print "good", "A", A,"Y", Y, "h", h
	y0 = [taylor(h,A[0],F,y0[0],k), taylor(h,A[1],F,y0[1],k)] #10-40
	return (h,y0)

y0 = [1.,1.] #example Moore1965interval p. 104
F = 'EXPLICITLY_DEFINED'
k = 9
eps = 0.000000001
p = [1.0,1.1,1.2,1.3,1.4]
while (y0[0] < 3.):
	(progress,y0) = step(F,k,y0, eps, p)
	print y0