diff --git a/Project.toml b/Project.toml
index d5622af1..90080a13 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,12 +1,12 @@
 name = "Roots"
 uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-version = "0.8.4"
+version = "1.0.0"
 
 [deps]
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [compat]
-SpecialFunctions = "0.8, 0.9"
+SpecialFunctions = "0.8, 0.9, 0.10"
 julia = "1.0"
 
 [extras]
diff --git a/src/bracketing.jl b/src/bracketing.jl
index ae3892d8..3985e042 100644
--- a/src/bracketing.jl
+++ b/src/bracketing.jl
@@ -90,8 +90,9 @@ function _init_state(method::AbstractBisection, fs, xs, fxs)
 
     x0, x1 = xs
     fx0, fx1 = fxs
-    state = UnivariateZeroState(x1, x0, [x1],
-                                fx1, fx0, [fx1],
+
+    state = UnivariateZeroState(x1, x0, zero(x1)/zero(x1)*oneunit(x1), [x1],
+                                fx1, fx0, fx1, [fx1],
                                 0, 2,
                                 false, false, false, false,
                                 "")
@@ -124,8 +125,10 @@ function init_state!(state::UnivariateZeroState{T,S}, M::AbstractBisection, fs,
             m, fm = x1, fx1
         end
         state.f_converged = true
-        state.xn1 = m
-        state.fxn1 = fm
+        state.xstar  = m
+        state.fxstar =  fm
+        #state.xn1 = m
+        #state.fxn1 = fm
         return state
     end
 
@@ -137,8 +140,15 @@ function init_state!(state::UnivariateZeroState{T,S}, M::AbstractBisection, fs,
         incfn(state)
         if sign(fx0) * sign(fm) < 0
             x1, fx1 = m, fm
-        else
+        elseif  sign(fx0) * sign(fm) > 0
             x0, fx0 = m, fm
+        else
+            state.message = "Exact zero found"
+            state.xstar = m
+            state.fxstar = fm
+            state.f_converged = true
+            state.x_converged = true
+            return state
         end
     end
 
@@ -234,7 +244,7 @@ function assess_convergence(M::Bisection, state::UnivariateZeroState{T,S}, optio
 
     if x_converged
         state.message=""
-        state.xn1 = xm
+        state.xstar = xm
         state.x_converged = true
         return true
     end
@@ -246,8 +256,8 @@ function assess_convergence(M::Bisection, state::UnivariateZeroState{T,S}, optio
 
     if f_converged
         state.message = ""
-        state.xn1 = xm
-        state.fxn1 = fm
+        state.xstar = xm
+        state.fxstar = fm
         state.f_converged = f_converged
         return true
     end
@@ -266,14 +276,17 @@ function assess_convergence(M::BisectionExact, state::UnivariateZeroState{T,S},
     for (c,fc) in ((x0,y0), (xm,ym), (x1, y1))
         if iszero(fc) || isnan(fc) #|| isinf(fc)
             state.f_converged = true
-            state.xn1 = c
-            state.fxn1 = fc
+            state.x_converged = true
+            state.xstar = c
+            state.fxstar = fc
             return true
         end
     end
 
     x0 < xm < x1 && return false
 
+    state.xstar = x1
+    state.fxstar = ym
     state.x_converged = true
     return true
 end
@@ -288,7 +301,7 @@ function update_state(method::Union{Bisection, BisectionExact}, fs, o::Univariat
     m::T = o.m[1]
     ym::S = o.fm[1] #fs(m)
 
-    if y0 * ym < 0
+    if sign(y0) * sign(ym) < 0
         o.xn1, o.fxn1 = m, ym
     else
         o.xn0, o.fxn0 = m, ym
@@ -336,7 +349,7 @@ function find_zero(fs, x0, method::M;
 
     verbose && show_trace(method, nothing, state, l)
 
-    state.xn1
+    state.xstar
 
 end
 
@@ -456,8 +469,8 @@ function init_state(M::AbstractAlefeldPotraShi, f, xs)
     end
     fu, fv = promote(f(u), f(v))
     isbracket(fu, fv) || throw(ArgumentError(bracketing_error))
-    state = UnivariateZeroState(v, u, [v, v], ## x1, x0, d, [ee]
-                                fv, fu, [fv,fv], ## fx1, fx0, d, [fe]
+    state = UnivariateZeroState(v, u, zero(v)/zero(v)*oneunit(v), [v, v], ## x1, x0, d, [ee]
+                                fv, fu, fv, [fv,fv], ## fx1, fx0, d, [fe]
                                 0, 2,
                                 false, false, false, false,
                                 "")
@@ -502,14 +515,14 @@ end
     default_tolerances(::AbstractAlefeldPotraShi, T, S)
 
 The default tolerances for Alefeld, Potra, and Shi methods are
-`xatol=zero(T)`, `xrtol=2eps(T)`, `atol= zero(S), and rtol=zero(S)`, with
+`xatol=zero(T)`, `xrtol=eps(T)/2`, `atol= zero(S), and rtol=zero(S)`, with
 appropriate units; `maxevals=45`, `maxfnevals = Inf`; and `strict=true`.
 
 """
 default_tolerances(M::AbstractAlefeldPotraShi) = default_tolerances(M, Float64, Float64)
 function default_tolerances(::AbstractAlefeldPotraShi, ::Type{T}, ::Type{S}) where {T,S}
     xatol = zero(T)
-    xrtol = 2 * eps(one(T))
+    xrtol = eps(one(T))/2
     atol = zero(float(one(S))) * oneunit(S)
     rtol = zero(float(one(S))) * one(S)
     maxevals = 45
@@ -536,8 +549,10 @@ function check_zero(::AbstractBracketing, state, c, fc)
     elseif iszero(fc)
         state.f_converged=true
         state.message *= "Exact zero found. "
-        state.xn1 = c
-        state.fxn1 = fc
+        state.xstar = c
+        state.fxstar =  fc
+#        state.xn1 = c
+#        state.fxn1 = fc
         return true
     end
     return false
@@ -545,7 +560,12 @@ end
 
 function assess_convergence(method::AbstractAlefeldPotraShi, state::UnivariateZeroState{T,S}, options) where {T,S}
 
-    (state.stopped || state.x_converged || state.f_converged) && return true
+    if state.stopped || state.x_converged || state.f_converged
+        if isnan(state.xstar)
+            state.xstar, state.fxstar = state.xn1, state.fxn1
+        end
+        return true
+    end
 
     if state.steps > options.maxevals
         state.stopped = true
@@ -565,21 +585,24 @@ function assess_convergence(method::AbstractAlefeldPotraShi, state::UnivariateZe
 
     if abs(fu) <= maximum(promote(options.abstol, abs(u) * oneunit(fu) / oneunit(u) * options.reltol))
         state.f_converged = true
-        state.xn1=u
-        state.fxn1=fu
+        state.xstar=u
+        state.fxstar=fu
         if iszero(fu)
+
+            state.x_converged
             state.message *= "Exact zero found. "
         end
         return true
     end
 
     a,b = state.xn0, state.xn1
-    tol = maximum(promote(options.xabstol, max(abs(a),abs(b)) * options.xreltol))
+    mx = max(abs(a), abs(b))
+    tol = maximum(promote(options.xabstol, mx * options.xreltol, sign(options.xreltol) *   eps(mx)))
 
     if abs(b-a) <= 2tol
         # use smallest of a,b,m
-        state.xn1 = u
-        state.fxn1 = fu
+        state.xstar = u
+        state.fxstar = fu
         state.x_converged = true
         return true
     end
@@ -755,8 +778,8 @@ function init_state(M::Brent, f, xs)
 
 
 
-    state = UnivariateZeroState(b, a, [a, a], ## x1, x0, c, d
-                                fb, fa, [fa, one(fa)], ## fx1, fx0, fc, mflag
+    state = UnivariateZeroState(b, a, zero(b)/zero(b)*oneunit(b), [a, a], ## x1, x0, c, d
+                                fb, fa, fb, [fa, one(fa)], ## fx1, fx0, fc, mflag
                                 0, 2,
                                 false, false, false, false,
                                 "")
@@ -944,3 +967,30 @@ end
         $f(fa, fb, fx)
     end
 end
+
+
+
+
+"""
+    find_bracket(f, x0, method=A42(); kwargs...)
+
+For bracketing methods returns an approximate root, the last bracketing interval used, and a flag indicating if an exact zero was found as a named tuple.
+
+With the default tolerances, one  of these should be the case: `exact`  is `true` (indicating termination  of the algorithm due to an exact zero  being identified) or the length of `bracket` is less or equal than `2eps(maximum(abs.(bracket)))`. In the `BisectionExact` case, the 2 could  be replaced by 1, as the bracket, `(a,b)` will satisfy  `nextfloat(a) == b `;  the Alefeld,  Potra, and Shi algorithms don't quite have  that promise.
+
+"""
+function find_bracket(fs, x0, method::M=A42(); kwargs...) where {M <: Union{AbstractAlefeldPotraShi, BisectionExact}}
+    x = adjust_bracket(x0)
+    T = eltype(x[1])
+    F = callable_function(fs)
+    state = init_state(method, F, x)
+    options = init_options(method, state; kwargs...)
+
+    # check if tolerances are exactly 0
+    iszero_tol = iszero(options.xabstol) && iszero(options.xreltol) && iszero(options.abstol) && iszero(options.reltol)
+
+    find_zero(method, F, options, state, NullTracks())
+
+    (xstar=state.xstar, bracket=(state.xn0, state.xn1), exact=iszero(state.fxstar))
+
+end
diff --git a/src/derivative_free.jl b/src/derivative_free.jl
index f427cea6..f9932ab5 100644
--- a/src/derivative_free.jl
+++ b/src/derivative_free.jl
@@ -62,8 +62,8 @@ end
 function init_state(method::AbstractSecant, fs, x::Union{Tuple, Vector})
     x0, x1 = promote(float(x[1]), float(x[2]))
     fx0, fx1 = fs(x0), fs(x1)
-    state = UnivariateZeroState(x1, x0, eltype(x1)[],
-                                fx1, fx0, eltype(fx1)[],
+    state = UnivariateZeroState(x1, x0, zero(x1)/zero(x1)*oneunit(x1), eltype(x1)[],
+                                fx1, fx0, fx1, eltype(fx1)[],
                                 0, 2,
                                 false, false, false, false, "")
 
diff --git a/src/find_zero.jl b/src/find_zero.jl
index 3f193602..a34f7ae9 100644
--- a/src/find_zero.jl
+++ b/src/find_zero.jl
@@ -39,9 +39,11 @@ abstract type  AbstractUnivariateZeroState end
 mutable struct UnivariateZeroState{T,S} <: AbstractUnivariateZeroState where {T,S}
     xn1::T
     xn0::T
+    xstar::T
     m::Vector{T}
     fxn1::S
     fxn0::S
+    fxstar::S
     fm::Vector{S}
     steps::Int
     fnevals::Int
@@ -61,9 +63,9 @@ function init_state(method::Any, fs, x)
     x1 = float(x)
     fx1 = fs(x1); fnevals = 1
     T, S = eltype(x1), eltype(fx1)
-
-    state = UnivariateZeroState(x1, oneunit(x1) * (0*x1)/(0*x1), T[],
-                                fx1, oneunit(fx1) * (0*fx1)/(0*fx1), S[],
+    zT, zS =  oneunit(x1) * (0*x1)/(0*x1), oneunit(fx1) * (0*fx1)/(0*fx1)
+    state = UnivariateZeroState(x1, zT, zT/Zt*oneunit(x1), T[],
+                                fx1, zS, zS, S[],
                                 0, fnevals,
                                 false, false, false, false,
                                 "")
@@ -112,8 +114,8 @@ end
 
 
 function Base.copy(state::UnivariateZeroState{T,S}) where {T, S}
-    UnivariateZeroState(state.xn1, state.xn0, copy(state.m),
-                        state.fxn1, state.fxn0, copy(state.fm),
+    UnivariateZeroState(state.xn1, state.xn0, state.xstar, copy(state.m),
+                        state.fxn1, state.fxn0, state.fxstar, copy(state.fm),
                         state.steps, state.fnevals,
                         state.stopped, state.x_converged,
                         state.f_converged, state.convergence_failed,
@@ -348,6 +350,9 @@ function assess_convergence(method::Any, state::UnivariateZeroState{T,S}, option
     fxn1 = state.fxn1
 
     if (state.x_converged || state.f_converged || state.stopped)
+        if isnan(state.xstar)
+            state.xstar, state.fxstar =  xn1, fxn1
+        end
         return true
     end
 
@@ -365,6 +370,7 @@ function assess_convergence(method::Any, state::UnivariateZeroState{T,S}, option
 
     # f(xstar) ≈ xstar * f'(xstar)*eps(), so we pass in lambda
     if _is_f_approx_0(fxn1, xn1, options.abstol, options.reltol)
+        state.xstar, state.fxstar = xn1, fxn1
         state.f_converged = true
         return true
     end
@@ -372,6 +378,7 @@ function assess_convergence(method::Any, state::UnivariateZeroState{T,S}, option
     # stop when xn1 ~ xn.
     # in find_zeros there is a check that f could be a zero with a relaxed tolerance
     if abs(xn1 - xn0) < max(options.xabstol, max(abs(xn1), abs(xn0)) * options.xreltol)
+        state.xstar, state.fxstar = xn1, fxn1
         state.message *= "x_n ≈ x_{n-1}. "
         state.x_converged = true
         return true
@@ -580,10 +587,22 @@ function find_zero(M::AbstractUnivariateZeroMethod,
     return decide_convergence(M, F, state, options)
 end
 
+function find_zero(M::AbstractUnivariateZeroMethod,
+                   F,
+                   state::AbstractUnivariateZeroState,
+                   l::AbstractTracks=NullTracks()
+                   )  #where {T<:Number, S<:Number}
+
+    options = init_options(M, state)
+    find_zero(M, F, options, state, l)
+end
+
+
+
 # state has stopped, this identifies if it has converged
 function decide_convergence(M::AbstractUnivariateZeroMethod,  F, state::UnivariateZeroState{T,S}, options) where {T,S}
-    xn1 = state.xn1
-    fxn1 = state.fxn1
+    xn1 = state.xstar
+    fxn1 = state.fxstar
 
     if (state.stopped || state.x_converged) && !(state.f_converged)
         ## stopped is a heuristic, x_converged can mask issues
@@ -611,6 +630,7 @@ function decide_convergence(M::AbstractUnivariateZeroMethod,  F, state::Univaria
         else
             xstar, fxstar = state.xn1, state.fxn1
             if _is_f_approx_0(fxstar, xstar, options.abstol, options.reltol, :relaxed)
+                state.xstar, state.fxstar = xstar, fxstar
                 msg = "Algorithm stopped early, but |f(xn)| < ϵ^(1/3), where ϵ depends on xn, rtol, and atol. "
                 state.message = state.message == "" ? msg : state.message * "\n\t" * msg
                 state.f_converged = true
@@ -621,10 +641,10 @@ function decide_convergence(M::AbstractUnivariateZeroMethod,  F, state::Univaria
     end
 
     if state.f_converged
-        return state.xn1
+        return state.xstar
     end
 
-    nan = NaN * state.xn1
+    nan = NaN * xn1
     if state.convergence_failed
         return nan
     end
@@ -685,6 +705,7 @@ function find_zero(M::AbstractUnivariateZeroMethod,
         ## did we find a zero or a bracketing interval?
         if iszero(state0.fxn1)
             copy!(state, state0)
+            state.xstar, state.fxstar = state.xn1, state.fxn1
             state.f_converged = true
             break
         elseif sign(state0.fxn0) * sign(state0.fxn1) < 0
diff --git a/src/newton.jl b/src/newton.jl
index c64dcd08..8701d90a 100644
--- a/src/newton.jl
+++ b/src/newton.jl
@@ -52,8 +52,9 @@ function init_state(method::AbstractNewtonLikeMethod, fs, x)
     fnevals = 1
     S = eltype(fx1)
 
-    state = UnivariateZeroState(x1, oneunit(x1) * (0*x1)/(0*x1), [Δ],
-                                fx1, oneunit(fx1) * (0*fx1)/(0*fx1), S[],
+    zT, zS = oneunit(x1) * (0*x1)/(0*x1), oneunit(fx1) * (0*fx1)/(0*fx1)
+    state = UnivariateZeroState(x1, zT, zT/zT*oneunit(x1), [Δ],
+                                fx1, zS, zS, S[],
                                 0, fnevals,
                                 false, false, false, false,
                                 "")
@@ -163,8 +164,9 @@ function init_state(method::AbstractHalleyLikeMethod, fs, x)
     S = eltype(fx1)
     fnevals = 3
 
-    state = UnivariateZeroState(x1, oneunit(x1) * (0*x1)/(0*x1), [Δ,ΔΔ],
-                                fx1, oneunit(fx1) * (0*fx1)/(0*fx1), S[],
+    zT, zS =  oneunit(x1) * (0*x1)/(0*x1), oneunit(fx1) * (0*fx1)/(0*fx1)
+    state = UnivariateZeroState(x1, zT, zT/zT*oneunit(x1), [Δ,ΔΔ],
+                                fx1, zS, zS, S[],
                                 0, fnevals,
                                 false, false, false, false,
                                 "")
diff --git a/test/test_bracketing.jl b/test/test_bracketing.jl
index c3a7565d..da64c2ae 100644
--- a/test/test_bracketing.jl
+++ b/test/test_bracketing.jl
@@ -229,6 +229,7 @@ avg(x) = sum(x)/length(x)
     ## brent has some failures
     Ms = [Roots.Brent()]
     results = [run_tests((f,b) -> find_zero(f, b, M), name="$M") for M in Ms]
+
     maxfailures = maximum([length(result.failures) for result in results])
     maxresidual = maximum([result.maxresidual for result in results])
     cnts = [result.evalcount for result in results]
@@ -251,6 +252,13 @@ avg(x) = sum(x)/length(x)
 
 end
 
+mutable struct Cnt
+    cnt::Int
+    f
+    Cnt(f) = new(0, f)
+end
+(f::Cnt)(x) = (f.cnt += 1; f.f(x))
+
 ## Some tests for FalsePosition methods
 @testset "FalsePosition" begin
     galadino_probs = [(x -> x^3 - 1, [.5, 1.5]),
diff --git a/test/test_find_zero.jl b/test/test_find_zero.jl
index 841dccf8..eb018081 100644
--- a/test/test_find_zero.jl
+++ b/test/test_find_zero.jl
@@ -317,3 +317,55 @@ end
     end
 
 end
+
+@testset "find_bracket test" begin
+    fs_xs = ((x -> x^5 - x - 1, (0,2)),
+             (sin, (3,4)),
+             (x -> exp(x) - x^4, (5, 20))
+             )
+
+    for (f, x0) in fs_xs
+        out = Roots.find_bracket(f, x0)
+        @test prod(f.(out.bracket)) <= 0
+    end
+
+    ## test size of  bracket
+    for M in  (Roots.BisectionExact(),  Roots.A42(), Roots.AlefeldPotraShi())
+        for (fn, b) in ((x  -> x == 0. ? 0. : x/exp(1/(x*x)), (-1.0, 4.0)),
+                        (x -> exp(-15*x)*(x-1) + x^15, (0., 1.)),
+                        (x -> (-40.0)*x*exp(-x),  (-9, 31))
+                        )
+            l =  Roots.find_bracket(fn, b, M)
+            @test  l.exact ||  abs(l.bracket[2] -  l.bracket[1])  <=  eps(l.xstar)
+        end
+    end
+
+    ## subnormal
+    x0 = nextfloat(nextfloat(0.0))
+    fn = x  ->  (x-x0)
+    b =  (0.0, 1.0)
+    m =  Roots.__middle(b...)  # 1e-154
+    M = Roots.BisectionExact()
+    l = Roots.find_bracket(fn,  b,  M)
+    ## Here l.exact=true, but  this checks the bracket size
+    @test  !(abs(-(l.bracket...)) <=  eps(l.xstar))
+    @test  abs(-(l.bracket...)) <=  m
+
+    M =  Roots.A42()
+    l = Roots.find_bracket(fn,  b,  M)
+    @test  !(abs(-(l.bracket...)) <=  eps(l.xstar))
+    @test  abs(-(l.bracket...)) <= m
+
+    M =  Roots.AlefeldPotraShi()
+    l = Roots.find_bracket(fn,  b,  M)
+    @test  !(abs(-(l.bracket...)) <=  eps(l.xstar))
+    @test  abs(-(l.bracket...)) <= m
+
+    x0 = nextfloat(0.0)
+    fn  = x ->  x <=  0.0 ?  x-x0  : x  +  x0
+    b = (-1.0, 1.0)
+    M = Roots.BisectionExact()
+    l = Roots.find_bracket(fn,  b,  M)
+    @test !l.exact &&   abs(-(l.bracket...))  <= maximum(eps.(l.bracket))
+
+end