Merge pull request #75 from jverzani/allocate

Allocate

README.md

@@ -121,4 +121,4 @@ Some additional documentation can be read [here](http://nbviewer.ipython.org/url
 ## Polynomials
 
 Special methods for finding roots of polynomials have been moved to
-the `PolynomialZeros` package and its `polyroots(f, domain)` function.
+the `PolynomialZeros` package and its `poly_roots(f, domain)` function.

src/Roots.jl

@@ -166,7 +166,7 @@ fzeros(p) = Base.depwarn("""
 Calling fzeros with just a polynomial is deprecated.
 Either:
    * Specify an interval to search over: fzeros(p, a, b).
-   * Use the `realroots` function from `PolynomialZeros`                         
+   * Use the `poly_roots(p, Over.R)` call from `PolynomialZeros`                         
    * Use `Polynomials` or `PolynomialRoots` and filter. For example, 
 
 ```

src/bracketing.jl

@@ -52,6 +52,64 @@ function _middle(x::Real, y::Real)
 end
 
 
+"""
+
+    `Roots.bisection64(f, a, b)` (unexported)
+
+* `f`: a callable object, like a function
+
+* `a`, `b`: Real values specifying a *bracketing* interval (one with
+`f(a) * f(b) < 0`). These will be converted to `Float64` values.
+
+Runs the bisection method using midpoints determined by a trick
+leveraging 64-bit floating point numbers. After ensuring the
+intermediate bracketing interval does not straddle 0, the "midpoint"
+is half way between the two values onces converted to unsigned 64-bit
+integers. This means no more than 64 steps will be taken, fewer if `a`
+and `b` already share some bits.
+
+The process is guaranteed to return a value `c` with `f(c)` one of
+`0`, `Inf`, or `NaN`; *or* one of `f(prevfloat(c))*f(c) < 0` or
+`f(c)*f(nextfloat(c)) > 0` holding. 
+
+This function is a bit faster than the slightly more general 
+`find_zero(f, [a,b], Bisection())` call.
+
+Due to Jason Merrill.
+
+"""
+function bisection64(f, a::Number, b::Number)
+
+    a,b = Float64(a), Float64(b)
+
+    if a > b
+        b,a = a, b
+    end
+
+
+    m = _middle(a,b)
+    fa, fb = sign(f(a)), sign(f(b))
+
+    fa * fb > 0 && throw(ArgumentError(bracketing_error)) 
+    iszero(fa) || isnan(fa) || isinf(fa) && return a
+    iszero(fb) || isnan(fb) || isinf(fb) && return b
+
+    while a < m < b
+        fm = sign(f(m))
+
+        if iszero(fm) || isnan(fm) || isinf(fm)
+            return m
+        elseif fa * fm < 0
+            b,fb=m,fm
+        else
+            a,fa=m,fm
+        end
+        m = _middle(a,b)
+    end
+    return m
+end
+
+
 ####
 ## find_zero interface. We need to specialize for T<:Float64, and BigSomething
 const BigSomething = Union{BigFloat, BigInt}
@@ -64,9 +122,21 @@ type A42 <: AbstractBisection end
 ## As we don't have the `A42` algorithm implemented through `find_zero`, we adjust here.
 function find_zero{M<:AbstractBisection, T<:Real}(f, x0::Vector{T}, method::M; maxevals::Int=50, verbose::Bool=false, kwargs...)
     x = sort(float(x0))
+
+    if isinf(x[1])
+        x[1] = nextfloat(x[1])
+    end
+    if isinf(x[2])
+        x[2] = prevfloat(x[2])
+    end
+
     if eltype(x) <: Float64
-        prob, options = derivative_free_setup(method, DerivativeFree(f), x; verbose=verbose, maxevals=maxevals, kwargs...)
-        find_zero(prob, method, options)
+        if verbose
+            prob, options = derivative_free_setup(method, DerivativeFree(f, f(x0[1])), x; verbose=verbose, maxevals=maxevals, kwargs...)
+            find_zero(prob, method, options)
+        else # avoid overhead of generic calling method
+            bisection64(f, x0[1], x0[2])::eltype(x)  
+        end
     else
         a42(f, x[1], x[2]; maxeval=maxevals, verbose=verbose)
     end
@@ -91,11 +161,12 @@ end
 ## The tolerances can be set to 0, in which case, the termination occurs when `nextfloat(x0) = x2`.
 ## The bracket `[a,b]` must be bounded.
 
-function init_state{T <: Float64}(method::Bisection, fs, x::Vector{T}, bracket)
+function init_state{T <: Float64,R}(method::Bisection, fs::DerivativeFree{R}, x::Vector{T}, bracket)
     x0, x2 = sort(x[1:2])
     isinf(x0) && (x0 = nextfloat(x0))
     isinf(x2) && (x2 = prevfloat(x2))
-    @compat y0, y2 = fs.f.([x0, x2])
+    y0::R = fs.f(x0)
+    y2::R = fs.f(x2)
 
     sign(y0) * sign(y2) > 0 && throw(ArgumentError(bracketing_error))
 
@@ -154,6 +225,8 @@ end
 
 """
 
+    `Roots.a42(f, a, b; kwargs...)` (not exported)
+
 Finds the root of a continuous function within a provided
 interval [a, b], without requiring derivatives. It is based on algorithm 4.2
 described in: 1. G. E. Alefeld, F. A. Potra, and Y. Shi, "Algorithm 748:
@@ -415,6 +488,7 @@ end
 
 
 """
+
 Searches for zeros  of `f` in an interval [a, b].
 
 Basic algorithm used:
@@ -427,10 +501,12 @@ If there are many zeros relative to the number of points, the process
 is repeated with more points, in hopes of finding more zeros for
 oscillating functions.
 
+Called by `fzeros` or `Roots.find_zeros`.
+
 """
 function find_zeros(f, a::Real, b::Real, args...;
                     no_pts::Int=100,
-                    ftol::Real=10*eps(), reltol::Real=10*eps(),
+                    abstol::Real=10*eps(), reltol::Real=10*eps(), ## should be abstol, reltol as used. 
                     kwargs...)
 
     a, b = a < b ? (a,b) : (b,a)
@@ -441,13 +517,13 @@ function find_zeros(f, a::Real, b::Real, args...;
     ## Look in [ai, bi)
     for i in 1:(no_pts+1)
         ai,bi=xs[i:i+1]
-        if isapprox(f(ai), 0.0, rtol=reltol, atol=ftol)
+        if isapprox(f(ai), 0.0, rtol=reltol, atol=abstol)
             push!(rts, ai)
         elseif sign(f(ai)) * sign(f(bi)) < 0
             push!(rts, find_zero(f, [ai, bi], Bisection()))
         else
             try
-                x = find_zero(f, ai + (0.5)* (bi-ai), Order8(); maxevals=10, reltol=ftol, xreltol=reltol)
+                x = find_zero(f, ai + (0.5)* (bi-ai), Order8(); maxevals=10, abstol=abstol, reltol=reltol)
                 if ai < x < bi
                     push!(rts, x)
                 end
@@ -456,11 +532,11 @@ function find_zeros(f, a::Real, b::Real, args...;
         end
     end
     ## finally, b?
-    isapprox(f(b), 0.0, rtol=reltol, atol=ftol) && push!(rts, b)
+    isapprox(f(b), 0.0, rtol=reltol, atol=abstol) && push!(rts, b)
 
     ## redo if it appears function oscillates alot in this interval...
     if length(rts) > (1/4) * no_pts
-        return(find_zeros(f, a, b, args...; no_pts = 10*no_pts, ftol=ftol, reltol=reltol, kwargs...))
+        return(find_zeros(f, a, b, args...; no_pts = 10*no_pts, abstol=abstol, reltol=reltol, kwargs...))
     else
         return(sort(rts))
     end

src/find_zero.jl

@@ -13,29 +13,32 @@
 @compat abstract type UnivariateZeroMethod end
 
 # container for callable objects; not really necessary, but has some value.
-@compat abstract type CallableFunction end
-immutable DerivativeFree <: CallableFunction
+@compat abstract type CallableFunction{R} end
+immutable DerivativeFree{R} <: CallableFunction{R}
     f
+    fx0::R
 end
 
-immutable FirstDerivative <: CallableFunction
+immutable FirstDerivative{R} <: CallableFunction{R}
     f
     fp
+    fx0::R
 end
 
 
-immutable SecondDerivative <: CallableFunction
+immutable SecondDerivative{R} <: CallableFunction{R}
     f
     fp
     fpp
+    fx0::R
 end
 
 ## allows override for automatic derivatives, see Newton
-function callable_function(m::UnivariateZeroMethod, f)
-    !isa(f, Tuple) && return DerivativeFree(f)
-    length(f) == 1 && return DerivativeFree(f[1])
-    length(f) == 2 && return FirstDerivative(f[1], f[2])
-    SecondDerivative(f[1], f[2], f[3])
+function callable_function(m::UnivariateZeroMethod, f, x0)
+    !isa(f, Tuple) && return DerivativeFree(f, f(x0))
+    length(f) == 1 && return DerivativeFree(f[1], f(x0))
+    length(f) == 2 && return FirstDerivative(f[1], f[2], f(x0))
+    SecondDerivative(f[1], f[2], f[3], f(x0))
 end
 
 
@@ -146,14 +149,15 @@ function derivative_free_setup{T<:AbstractFloat}(method::Any, fs::CallableFuncti
     )
     x = float(x0)
     prob = UnivariateZeroProblem(fs, x, isa(bracket, Nullable) ? bracket : Nullable(convert(Vector{T}, bracket)))
-    options = univariate_zero_options(;xabstol=xabstol,
-                                      xreltol=xreltol,
-                                      abstol=abstol,
-                                      reltol=reltol,
-                                      maxevals=maxevals,
-                                      maxfnevals=maxfnevals,
-                                      verbose=verbose,
-                                      kwargs...)
+    # options = univariate_zero_options(;xabstol=xabstol,
+    #                                   xreltol=xreltol,
+    #                                   abstol=abstol,
+    #                                   reltol=reltol,
+    #                                   maxevals=maxevals,
+    #                                   maxfnevals=maxfnevals,
+    #                                   verbose=verbose,
+    #                                   kwargs...)
+    options = UnivariateZeroOptions(xabstol, xreltol, abstol, reltol, maxevals, maxfnevals, verbose)
     prob, options
 end
 
@@ -399,7 +403,7 @@ find_zero(f, 1.0, Steffensen()) # also Order2()
 ```
 """
 find_zero{T<:Number}(f, x0::Union{T,Vector{T}}, method::UnivariateZeroMethod; kwargs...) =
-    find_zero(method, callable_function(method, f), x0; kwargs...)
+    find_zero(method, callable_function(method, f, float(x0)), x0; kwargs...)
 
 ## some defaults for methods
 find_zero{T <: Number}(f, x0::T; kwargs...) = find_zero(f, x0, Order0(); kwargs...)
@@ -418,7 +422,7 @@ end
 
 ## old interface for fzero
 ## old keyword arguments (see ?fzero) handled in univariate_zero_options
-function derivative_free{T <: AbstractFloat}(f, x0::T; order::Int=0,
+@noinline function derivative_free{T <: AbstractFloat}(f, x0::T; order::Int=0,
                                              kwargs...)
 
     if order == 0

src/newton.jl

@@ -7,11 +7,11 @@
 type Newton <: UnivariateZeroMethod
 end
 
-function callable_function(method::Newton, f::Tuple)
+function callable_function(method::Newton, f::Tuple, x0)
     length(f) == 1 && return FirstDerivative(f[1], D(f[1]))
-    FirstDerivative(f[1], f[2])
+    FirstDerivative(f[1], f[2], f[1](x0))
 end
-callable_function(method::Newton, f::Any) = FirstDerivative(f, D(f))
+callable_function(method::Newton, f::Any, x0) = FirstDerivative(f, D(f), x0)
 
 function init_state{T}(method::Newton, fs, x0::T)
     state = UnivariateZeroState(x0,
@@ -143,12 +143,12 @@ newton(f, fp, x0; kwargs...) = find_zero((f, fp), x0, Newton(); kwargs...)
 type Halley <: UnivariateZeroMethod
 end
 
-function callable_function(method::Halley, f::Tuple)
-    length(f) == 1 && return SecondDerivative(f[1], D(f[1]), D(f[1],2))
-    length(f) == 2 && return SecondDerivative(f[1], f[2], D(f[2],1))
-    SecondDerivative(f[1], f[2], f[3])
+function callable_function(method::Halley, f::Tuple, x0)
+    length(f) == 1 && return SecondDerivative(f[1], D(f[1]), D(f[1],2), f[1](x0))
+    length(f) == 2 && return SecondDerivative(f[1], f[2], D(f[2],1), f[1](x0))
+    SecondDerivative(f[1], f[2], f[3], f[1](x0))
 end
-callable_function(method::Halley, f) = SecondDerivative(f, D(f), D(f, 2))
+callable_function(method::Halley, f, x0) = SecondDerivative(f, D(f), D(f, 2), f(x0))
 
 
 function update_state{T}(method::Halley, fs, o::UnivariateZeroState{T}, options::UnivariateZeroOptions)

test/test_fzero.jl

@@ -19,7 +19,7 @@ end
 fn, xstar, x0, br = x -> x^5 - x - 1, 1.1673039782614187, 1.0, [1.0, 2.0]
 @test fzero(fn, x0, order=1)  ≈ xstar
 
-@test_throws Roots.ConvergenceFailed fzero(fn, x0, order=1, maxeval=2)
+@test_throws Roots.ConvergenceFailed fzero(fn, x0, order=1, maxevals=2)
 #@test norm(fzero(fn, x0, order=1, ftol=1e-2) - xstar) > 1e-5
 #@test norm(fzero(fn, x0, order=1, xtol=1e-2) - xstar) > 1e-10