diff --git a/Makefile.am b/Makefile.am index 27d9dcda..a7c8b7fc 100644 --- a/Makefile.am +++ b/Makefile.am @@ -279,6 +279,7 @@ UTIL_FILES = \ region.cpp \ region.hpp \ region.ipp \ + roots/bisection.hpp \ si.cpp \ signal.cpp \ signal.hpp \ @@ -433,6 +434,7 @@ TEST_BIN = \ test/rational \ test/region \ test/ripemd \ + test/roots/bisection \ test/sha1 \ test/sha2 \ test/signal \ diff --git a/roots/bisection.hpp b/roots/bisection.hpp new file mode 100644 index 00000000..b6e8935f --- /dev/null +++ b/roots/bisection.hpp @@ -0,0 +1,69 @@ +/* + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * Copyright 2016 Danny Robson + */ + +#ifndef __UTIL_ROOTS_BISECTION_HPP +#define __UTIL_ROOTS_BISECTION_HPP + +#include "../maths.hpp" + +#include + +namespace util { namespace roots { + /// find a root of a function using the bisection method + /// + /// the user is responsible for ensuring there is in fact a root. there + /// is no iteration limit currently, so poorly considered arguments could + /// result in an infinite loop. + /// + /// a: lower bound + /// b: upper bound + /// f: real function + /// tolerance: minimum range termination condition + template + T + bisection (T a, T b, T (*f)(T), T tolerance) + { + // to simplify the exit conditions we assume that a..b is an + // increasing range. + if (a > b) + std::swap (a, b); + + while (b - a > tolerance) { + // calculate the midpoint + auto c = (a + b) / T{2}; + auto f_c = f (c); + + // return early if we happened across the root + if (exactly_zero (f_c)) + return c; + + // check which interval we fell into, and divide for the next + // iteration + auto f_a = f (a); + if (samesign (f_a, f_c)) + a = c; + else + b = c; + } + + // we didn't find it exactly, but it's somewhere in this range so take + // a punt and return the middle of the range. not sure if this is + // valid mathematically. + return (a + b) / T{2}; + } +} } + +#endif diff --git a/test/roots/bisection.cpp b/test/roots/bisection.cpp new file mode 100644 index 00000000..10a5e645 --- /dev/null +++ b/test/roots/bisection.cpp @@ -0,0 +1,40 @@ +#include "roots/bisection.hpp" + +#include "tap.hpp" + +#include "maths.hpp" + +using util::pow; + + +constexpr float order2 (float x) { return x * x + 3 * x - 7.f; } +constexpr float order4 (float x) { return 10 * pow (x, 4u) + -270 * pow (x, 2u) + -140 * pow (x, 1u) + +1200; } + +struct { + float (*func) (float); + float lo; + float hi; + float root; + const char *msg; +} TESTS[] = { + { order2, 0.f, 3.f, 0.5f * (std::sqrt (37.f) - 3.f), "order-2 bisection" }, + { order4, 0.f, 5.f, 2.f, "order-4 bisection" } +}; + + +int +main (void) +{ + util::TAP::logger tap; + + for (const auto &t: TESTS) { + constexpr float TOLERANCE = 0.00001f; + auto root = util::roots::bisection (t.lo, t.hi, t.func, TOLERANCE); + tap.expect_eq (root, t.root, t.msg); + } + + return tap.status (); +}