From 0b7a5c1a9ed29706f487f01a448396b1db7e8c62 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Sun, 19 Mar 2017 08:41:47 -0600
Subject: [PATCH 01/16] Add 1D bisection root-finder

---
 src/IntervalRootFinding.jl | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl
index 31c171f9..73bfb763 100644
--- a/src/IntervalRootFinding.jl
+++ b/src/IntervalRootFinding.jl
@@ -14,6 +14,7 @@ export
     Root, is_unique,
     find_roots,
     find_roots_midpoint,
+    bisection,
     bisect
 
 import Base: ⊆, show
@@ -21,6 +22,8 @@ import Base: ⊆, show
 const derivative = ForwardDiff.derivative
 const D = derivative
 
+
+# Root object:
 immutable Root{T<:Real}
     interval::Interval{T}
     status::Symbol
@@ -34,7 +37,26 @@ is_unique{T}(root::Root{T}) = root.status == :unique
 ⊆(a::Root, b::Root) = a.interval ⊆ b.interval
 
 
+# Common functionality:
+
+doc"""Returns the midpoint of the interval x, slightly shifted in case
+the midpoint is an exact root"""
+function guarded_mid{T}(f::Function, x::Interval{T})
+    m = mid(x)
+
+    if f(m) == zero(T)                      # midpoint exactly a root
+        α = convert(T, 0.46875)      # close to 0.5, but exactly representable as a floating point
+        m = α*x.lo + (one(T)-α)*x.hi   # displace to another point in the interval
+    end
+
+    @assert m ∈ x
+
+    return m
+end
+
+
 include("bisect.jl")
+include("bisection.jl")
 include("newton.jl")
 include("krawczyk.jl")
 

From 844d135d88d67bb7724afdf0d14bda4d48d33704 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Sun, 19 Mar 2017 08:48:33 -0600
Subject: [PATCH 02/16] Add bisection.jl file

---
 src/root_finding/bisection.jl | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)
 create mode 100644 src/root_finding/bisection.jl

diff --git a/src/root_finding/bisection.jl b/src/root_finding/bisection.jl
new file mode 100644
index 00000000..21b8298c
--- /dev/null
+++ b/src/root_finding/bisection.jl
@@ -0,0 +1,19 @@
+
+function bisection{T}(f, X::Interval{T}; tolerance=1e-3)
+
+    image = f(X)
+
+    if 0 ∉ image
+        return Root{T}[]  # guaranteed that no zero in the interval
+    end
+
+    if diam(X) < tolerance
+        return [Root{T}(X, :unknown)]
+    end
+
+    X1, X2 = bisect(X)
+
+    return [bisection(f, X1, tolerance=tolerance);
+            bisection(f, X2, tolerance=tolerance)]
+
+end

From c0801364842594f8f16d365ddcd3b0cde43609eb Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Sun, 19 Mar 2017 11:55:16 -0600
Subject: [PATCH 03/16] Rewrite bisection so that applies directly to
 IntervalBox

---
 src/root_finding/bisection.jl | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/root_finding/bisection.jl b/src/root_finding/bisection.jl
index 21b8298c..23502ecf 100644
--- a/src/root_finding/bisection.jl
+++ b/src/root_finding/bisection.jl
@@ -1,10 +1,15 @@
 
-function bisection{T}(f, X::Interval{T}; tolerance=1e-3)
+doc"""
+    bisection(f, X; tolerance=1e-3)
+
+Find possible roots of the function `f` inside the `Interval` or `IntervalBox` `X`.
+"""
+function bisection(f, X::Union{Interval, IntervalBox}; tolerance=1e-3)
 
     image = f(X)
 
-    if 0 ∉ image
-        return Root{T}[]  # guaranteed that no zero in the interval
+    if !(zero(X) ⊆ image)
+        return Root{T}[]
     end
 
     if diam(X) < tolerance

From eb287e7a762d53927f9475ba34e8621f41a06612 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Sun, 19 Mar 2017 12:56:54 -0600
Subject: [PATCH 04/16] Modifications to make bisection work with
 higher-dimensional functions

---
 src/IntervalRootFinding.jl    |  5 ++---
 src/root_finding/bisection.jl | 10 ++++++++--
 2 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl
index 73bfb763..32a486f2 100644
--- a/src/IntervalRootFinding.jl
+++ b/src/IntervalRootFinding.jl
@@ -22,10 +22,9 @@ import Base: ⊆, show
 const derivative = ForwardDiff.derivative
 const D = derivative
 
-
 # Root object:
-immutable Root{T<:Real}
-    interval::Interval{T}
+immutable Root{T}
+    interval::T
     status::Symbol
 end
 
diff --git a/src/root_finding/bisection.jl b/src/root_finding/bisection.jl
index 23502ecf..6d588a22 100644
--- a/src/root_finding/bisection.jl
+++ b/src/root_finding/bisection.jl
@@ -4,20 +4,26 @@ doc"""
 
 Find possible roots of the function `f` inside the `Interval` or `IntervalBox` `X`.
 """
-function bisection(f, X::Union{Interval, IntervalBox}; tolerance=1e-3)
+function bisection{T}(f, X::T; tolerance=1e-3, debug=false)
 
     image = f(X)
 
+    debug && @show X, image
+
     if !(zero(X) ⊆ image)
         return Root{T}[]
     end
 
     if diam(X) < tolerance
-        return [Root{T}(X, :unknown)]
+        return [Root(X, :unknown)]
     end
 
     X1, X2 = bisect(X)
 
+    if debug
+        @show X1, X2
+    end
+
     return [bisection(f, X1, tolerance=tolerance);
             bisection(f, X2, tolerance=tolerance)]
 

From 11c36fdfa0058b7aee083089fda8fa9878b8e0da Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Mon, 20 Mar 2017 12:48:17 -0600
Subject: [PATCH 05/16] Restrict signature of bisection

---
 src/root_finding/bisection.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/root_finding/bisection.jl b/src/root_finding/bisection.jl
index 6d588a22..97404eb6 100644
--- a/src/root_finding/bisection.jl
+++ b/src/root_finding/bisection.jl
@@ -4,7 +4,7 @@ doc"""
 
 Find possible roots of the function `f` inside the `Interval` or `IntervalBox` `X`.
 """
-function bisection{T}(f, X::T; tolerance=1e-3, debug=false)
+function bisection{T<:Union{Interval,IntervalBox}}(f, X::T; tolerance=1e-3, debug=false)
 
     image = f(X)
 

From c903baa57f242d8ddb9e179e781695cbd9b44892 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Mon, 20 Mar 2017 12:51:25 -0600
Subject: [PATCH 06/16] Restrict type for Root object

---
 src/IntervalRootFinding.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl
index 32a486f2..da3fd0f2 100644
--- a/src/IntervalRootFinding.jl
+++ b/src/IntervalRootFinding.jl
@@ -23,7 +23,7 @@ const derivative = ForwardDiff.derivative
 const D = derivative
 
 # Root object:
-immutable Root{T}
+immutable Root{T<:Union{Interval,IntervalBox}}
     interval::T
     status::Symbol
 end

From e9d6e37b19b8f47ecca825dcb6ea8876fbd031a4 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Mon, 20 Mar 2017 13:11:01 -0600
Subject: [PATCH 07/16] Add tests for bisection on infinite 1D and 2D intervals

---
 test/bisect.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/bisect.jl b/test/bisect.jl
index 72553d8a..09599ff5 100644
--- a/test/bisect.jl
+++ b/test/bisect.jl
@@ -2,7 +2,7 @@ using IntervalArithmetic, IntervalRootFinding
 using Base.Test
 
 
-@testset "Bisection tests" begin
+@testset "`bisect` function" begin
     X = 0..1
     @test bisect(X) == (0..0.5, 0.5..1)
     @test bisect(X, 0.25) == (0..0.25, 0.25..1)

From d5f58ae438664d43912bb4eb057eb8e8c0d1b809 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Mon, 20 Mar 2017 13:15:43 -0600
Subject: [PATCH 08/16] Fix failing test

---
 test/findroots.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/findroots.jl b/test/findroots.jl
index d8f15339..ccdc5198 100644
--- a/test/findroots.jl
+++ b/test/findroots.jl
@@ -86,7 +86,7 @@ function_list = [
                                             for i in 1:length(roots)
                                                 root = roots[i]
 
-                                                @test isa(root, Root{prec[1]})
+                                                @test isa(root, Root)
                                                 @test is_unique(root)
                                                 @test true_roots[i] ⊆ root.interval
                                             end

From 46c1c2d57ce35a9344457efa69207aea8922cb89 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Mon, 20 Mar 2017 18:36:41 -0600
Subject: [PATCH 09/16] Add bisection test file

---
 test/root_finding_tests/bisection.jl | 34 ++++++++++++++++++++++++++++
 1 file changed, 34 insertions(+)
 create mode 100644 test/root_finding_tests/bisection.jl

diff --git a/test/root_finding_tests/bisection.jl b/test/root_finding_tests/bisection.jl
new file mode 100644
index 00000000..e602e89e
--- /dev/null
+++ b/test/root_finding_tests/bisection.jl
@@ -0,0 +1,34 @@
+using ValidatedNumerics, ValidatedNumerics.RootFinding
+using Base.Test
+
+@testset "Bisection tests" begin
+    @testset "Single variable" begin
+        f(x) = sin(x) + cos(x/2) + x/5
+
+        roots = bisection(f, -∞..∞, tolerance=1e-3)
+        roots2 = [root.interval for root in roots]
+
+        @test roots2 ==  [  Interval(-0.8408203124999999, -0.8398437499999999),
+                            Interval(3.6425781249999996, 3.6435546874999996),
+                            Interval(6.062499999999999, 6.063476562499999)
+                         ]
+
+    end
+
+    @testset "Two variables" begin
+        @intervalbox g(x, y) = (x^2 + y^2 - 1, y - x)
+
+        X = (-∞..∞) × (-∞..∞)
+        roots = bisection(g, X, tolerance=1e-3)
+        roots2 = [root.interval for root in roots]
+        
+        @test roots2 == [IntervalBox(Interval(-0.7080078124999999, -0.7070312499999999), Interval(-0.7080078124999999, -0.7070312499999999)),
+ IntervalBox(Interval(-0.7080078124999999, -0.7070312499999999), Interval(-0.7070312499999999, -0.7060546874999999)),
+ IntervalBox(Interval(-0.7070312499999999, -0.7060546874999999), Interval(-0.7080078124999999, -0.7070312499999999)),
+ IntervalBox(Interval(0.7060546874999999, 0.7070312499999999), Interval(0.7070312499999999, 0.7080078124999999)),
+ IntervalBox(Interval(0.7070312499999999, 0.7080078124999999), Interval(0.7060546874999999, 0.7070312499999999)),
+ IntervalBox(Interval(0.7070312499999999, 0.7080078124999999), Interval(0.7070312499999999, 0.7080078124999999))]
+
+    end
+
+end

From a2aae73372105c393ded1e45a97c266b3f30d1c7 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Mon, 20 Mar 2017 18:41:14 -0600
Subject: [PATCH 10/16] Modify type of empty vector of Roots

---
 src/krawczyk.jl               | 4 ++--
 src/newton.jl                 | 4 ++--
 src/root_finding/bisection.jl | 2 +-
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/krawczyk.jl b/src/krawczyk.jl
index 99fbf0e4..974565dc 100644
--- a/src/krawczyk.jl
+++ b/src/krawczyk.jl
@@ -66,10 +66,10 @@ function krawczyk{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=
     # Maximum level of bisection
     level >= maxlevel && return [Root(x, :unknown)]
 
-    isempty(x) && return Root{T}[]  # [(x, :none)]
+    isempty(x) && return Root{typeof(x)}[]  # [(x, :none)]
     Kx = K(f, f_prime, x) ∩ x
 
-    isempty(Kx ∩ x) && return Root{T}[]  # [(x, :none)]
+    isempty(Kx ∩ x) && return Root{typeof(x)}[]  # [(x, :none)]
 
     if Kx ⪽ x  # isinterior
         debug && (print("Refining "); @show(x))
diff --git a/src/newton.jl b/src/newton.jl
index cdce49bd..1f33f398 100644
--- a/src/newton.jl
+++ b/src/newton.jl
@@ -64,7 +64,7 @@ function newton{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
 
     debug && (print("Entering newton:"); @show(level); @show(x))
 
-    isempty(x) && return Root{T}[]  #[(x, :none)]
+    isempty(x) && return Root{typeof(x)}[]  #[(x, :none)]
 
     z = zero(x.lo)
     tolerance = abs(tolerance)
@@ -79,7 +79,7 @@ function newton{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
         Nx = N(f, x, deriv)
         debug && @show(Nx, Nx ⊆ x, Nx ∩ x)
 
-        isempty(Nx ∩ x) && return Root{T}[]
+        isempty(Nx ∩ x) && return Root{typeof(x)}[]
 
         if Nx ⊆ x
             debug && (print("Refining "); @show(x))
diff --git a/src/root_finding/bisection.jl b/src/root_finding/bisection.jl
index 97404eb6..8eb3fd3b 100644
--- a/src/root_finding/bisection.jl
+++ b/src/root_finding/bisection.jl
@@ -11,7 +11,7 @@ function bisection{T<:Union{Interval,IntervalBox}}(f, X::T; tolerance=1e-3, debu
     debug && @show X, image
 
     if !(zero(X) ⊆ image)
-        return Root{T}[]
+        return Root{typeof(X)}[]
     end
 
     if diam(X) < tolerance

From eec1811dfbb1ca665b51b7237009f609d6bf2464 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Sun, 23 Apr 2017 19:46:45 -0500
Subject: [PATCH 11/16] Move bisection files to correct place

---
 src/{root_finding => }/bisection.jl        | 0
 test/{root_finding_tests => }/bisection.jl | 0
 2 files changed, 0 insertions(+), 0 deletions(-)
 rename src/{root_finding => }/bisection.jl (100%)
 rename test/{root_finding_tests => }/bisection.jl (100%)

diff --git a/src/root_finding/bisection.jl b/src/bisection.jl
similarity index 100%
rename from src/root_finding/bisection.jl
rename to src/bisection.jl
diff --git a/test/root_finding_tests/bisection.jl b/test/bisection.jl
similarity index 100%
rename from test/root_finding_tests/bisection.jl
rename to test/bisection.jl

From 884c6ca867642636d33739dc4cda66ee6d4b041e Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Sun, 23 Apr 2017 20:01:38 -0500
Subject: [PATCH 12/16] Add complex roots example

---
 examples/complex_roots.jl | 13 +++++++++++++
 1 file changed, 13 insertions(+)
 create mode 100644 examples/complex_roots.jl

diff --git a/examples/complex_roots.jl b/examples/complex_roots.jl
new file mode 100644
index 00000000..100c61e4
--- /dev/null
+++ b/examples/complex_roots.jl
@@ -0,0 +1,13 @@
+# Example of finding roots of a complex function
+# f(z) = z^4 + z - 2  (book of Hansen-Walster, chap. 10)
+
+using IntervalRootFinding, IntervalArithmetic
+
+@intervalbox f(x, y) = (z = x + im*y;
+                        z2 = z^4 + z - 2;
+                        reim(z2)
+                        )
+
+roots = bisection(f, (-10..10) × (-10..10) )
+
+

From 49dbd294a22ff2e2e57acc5f32b25273b76d7b05 Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Sun, 23 Apr 2017 20:13:19 -0500
Subject: [PATCH 13/16] Complex bisection returns complex intervals

---
 examples/complex_roots.jl | 22 +++++++++++++++++-----
 1 file changed, 17 insertions(+), 5 deletions(-)

diff --git a/examples/complex_roots.jl b/examples/complex_roots.jl
index 100c61e4..084434f1 100644
--- a/examples/complex_roots.jl
+++ b/examples/complex_roots.jl
@@ -3,11 +3,23 @@
 
 using IntervalRootFinding, IntervalArithmetic
 
-@intervalbox f(x, y) = (z = x + im*y;
-                        z2 = z^4 + z - 2;
-                        reim(z2)
-                        )
+# Find complex roots by bisection
+function complex_bisection(f, X::IntervalBox)
 
-roots = bisection(f, (-10..10) × (-10..10) )
+    function g(X::IntervalBox)
+        x, y = X
+        z = x + im*y;
+        zz = f(z)
+        x, y = reim(zz)
+        return IntervalBox(x, y)
+    end
 
+    roots = bisection(g, X)
+
+    return [root.interval[1] + im*root.interval[2] for root in roots]
+end
+
+f(z) = z^4 + z - 2
+
+complex_bisection(f, (-100..100) × (-100..100) )
 

From 888825fc27ce09c99f93c85df4d39288e87916ae Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Mon, 24 Apr 2017 01:06:41 -0500
Subject: [PATCH 14/16] Add method for complex_bisection with two intervals
 instead of IntervalBox

---
 examples/complex_roots.jl | 31 ++++++++++++++++++++++++++++---
 1 file changed, 28 insertions(+), 3 deletions(-)

diff --git a/examples/complex_roots.jl b/examples/complex_roots.jl
index 084434f1..eedcbf27 100644
--- a/examples/complex_roots.jl
+++ b/examples/complex_roots.jl
@@ -3,9 +3,26 @@
 
 using IntervalRootFinding, IntervalArithmetic
 
-# Find complex roots by bisection
+"""
+    complex_bisection(f, X)
+
+Find complex roots of $f: \mathbb{C} \to \mathbb{C}$.
+
+Inputs:
+
+- `f`: function that takes $z \in \mathbb{C}$ and returns another
+complex number.
+
+- `X`: An `IntervalBox` specifying the bounds on the real and imaginary parts
+of `z`.
+
+"""
 function complex_bisection(f, X::IntervalBox)
 
+    """
+    Make a 2D real version of the complex function `f` suitable for `bisection`,
+    i.e. that accepts an `IntervalBox` and returns an `IntervalBox`
+    """
     function g(X::IntervalBox)
         x, y = X
         z = x + im*y;
@@ -16,10 +33,18 @@ function complex_bisection(f, X::IntervalBox)
 
     roots = bisection(g, X)
 
-    return [root.interval[1] + im*root.interval[2] for root in roots]
+    return g, [root.interval[1] + im*root.interval[2] for root in roots]
 end
 
+"""
+    complex_bisection(f, x, y)
+
+Version in which the bounds are specified as two separate `Interval`s.
+"""
+complex_bisection(f, x::Interval, y::Interval) = complex_bisection(f, x × y)
+
 f(z) = z^4 + z - 2
 
-complex_bisection(f, (-100..100) × (-100..100) )
+L = 10
+g, roots = complex_bisection(f, -L..L, -L..L)
 

From b3ef09fe591bf171e0bc821cc6833e645074935e Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Mon, 8 May 2017 14:08:26 -0500
Subject: [PATCH 15/16] Improve complex_bisection using Complex constructor

---
 examples/complex_roots.jl | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/examples/complex_roots.jl b/examples/complex_roots.jl
index eedcbf27..8c4e14de 100644
--- a/examples/complex_roots.jl
+++ b/examples/complex_roots.jl
@@ -3,15 +3,14 @@
 
 using IntervalRootFinding, IntervalArithmetic
 
-"""
+doc"""
     complex_bisection(f, X)
 
 Find complex roots of $f: \mathbb{C} \to \mathbb{C}$.
 
 Inputs:
 
-- `f`: function that takes $z \in \mathbb{C}$ and returns another
-complex number.
+- `f`: function $f: \mathbb{C} \to \mathbb{C}$, i.e. a function that accepts a complex number and returns another complex number.
 
 - `X`: An `IntervalBox` specifying the bounds on the real and imaginary parts
 of `z`.
@@ -24,16 +23,15 @@ function complex_bisection(f, X::IntervalBox)
     i.e. that accepts an `IntervalBox` and returns an `IntervalBox`
     """
     function g(X::IntervalBox)
-        x, y = X
-        z = x + im*y;
+        z = Complex(X...)
         zz = f(z)
-        x, y = reim(zz)
-        return IntervalBox(x, y)
+
+        return IntervalBox(reim(zz))
     end
 
     roots = bisection(g, X)
 
-    return g, [root.interval[1] + im*root.interval[2] for root in roots]
+    return g, [Complex(root.interval...) for root in roots]
 end
 
 """
@@ -48,3 +46,4 @@ f(z) = z^4 + z - 2
 L = 10
 g, roots = complex_bisection(f, -L..L, -L..L)
 
+println(roots)

From 8734c278cec53cd082c1267ca448be587ec8d2ce Mon Sep 17 00:00:00 2001
From: David Sanders <dpsanders@gmail.com>
Date: Mon, 8 May 2017 14:50:56 -0500
Subject: [PATCH 16/16] Add bisection for vector of IntervalBox

---
 src/bisection.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/bisection.jl b/src/bisection.jl
index 8eb3fd3b..433d3476 100644
--- a/src/bisection.jl
+++ b/src/bisection.jl
@@ -28,3 +28,10 @@ function bisection{T<:Union{Interval,IntervalBox}}(f, X::T; tolerance=1e-3, debu
             bisection(f, X2, tolerance=tolerance)]
 
 end
+
+
+function bisection{T<:Union{Interval,IntervalBox}}(f, V::Vector{T}; tolerance=1e-3, debug=false)
+
+    return vcat([bisection(f, X, tolerance=tolerance, debug=debug) for X in V])
+
+end