module integ
implicit none
integer nrec

contains

function adapt(f,a,b,acc,eps) result(q)
interface
function f(x)
real*8 f,x
end function f
end interface
real*8 q,a,b,acc,eps,f2,f3
f2=f(a+2*(b-a)/6)
f3=f(a+4*(b-a)/6)
nrec=0
q=adapt23(f,a,b,acc,eps,f2,f3)
end function adapt

recursive function adapt23(f,a,b,acc,eps,f2,f3) result(q)
interface
function f(x)
real*8 f,x
end function f
end interface
real*8 a,b,acc,eps,f1,f2,f3,f4,q,r,err
nrec=nrec+1
if(nrec>1000000)then
	q=(a-a)/(a-a) ! NaN
	return
else
	f1=f(a+(b-a)/6)
	f4=f(a+5*(b-a)/6)
	q=(2*f1+f2+f3+2*f4)/6*(b-a)
	r=(f1+f4+f2+f3)/4*(b-a)
	if(abs(q-r)/2<acc+eps*abs(q))then
		return
	else
		q=adapt23(f,a,(a+b)/2,acc/sqrt(2.),eps,f1,f2) &
		+adapt23(f,(a+b)/2,b,acc/sqrt(2.),eps,f3,f4)
		return
	endif
endif
end function adapt23

end module integ