Self adjusting numerical integration precision

Question

I've written a program to numerically calculate the area under a curve for a given interval and want to be able to give a precision goal.

Let $f: [a, b] \rightarrow \mathbb{R}$ with $a,b \in \mathbb{R}$ and $a<b$ be an integrable function, quadrature(f, a, b) the Gauss-Legendre-Quadrature, mid(a,b)=a + 1/2 (b-a) and rel_error(v, w) = abs(1 - abs(v/w)) where abs(x) = |x|

I then calculate the area as follows:

  function integrate(f, a, b)
    total = quadrature(f, a, b)
    left  = quadrature(f, a, mid(a,b))
    right = quadrature(f, mid(a,b), b)

    if: rel_error(total, left+right) > epsilon
      return integrate(f, left, a, mid(a,b)) + integrate(f, right, mid(a,b), b)
    endif

    return total

  function integrate(f, total, a, b)
    left  = quadrature(f, a, mid(a, b))
    right = quadrature(f, mid(a, b), b)

    if: rel_error(total, left+right) > epsilon
      return integrate(f, left, a, mid(a,b)) + integrate(f, right, mid(a,b), b)
    endif

    return total

Note: I excluded some validations and a max_recursion limit for the sake of clarity.

To calculate an area I start with the top function, which may recursively call the bottom function, i.e. integrate(x^2, 0, 10) to approximate $\int \limits_0^{10} x^2 \mathrm{d}x$

My questions are:

Is this a solid way that will converge? How does $\epsilon$ relate to the overall precision? What value should I choose for $\epsilon$ so that integrate(f, a, b) differs by at most $10^{-6}$ from the exact value? You may assume $n=10$ for the quadrature.