diff options
Diffstat (limited to 'boost/math/tools/remez.hpp')
| -rw-r--r-- | boost/math/tools/remez.hpp | 667 |
1 files changed, 0 insertions, 667 deletions
diff --git a/boost/math/tools/remez.hpp b/boost/math/tools/remez.hpp deleted file mode 100644 index afcbe25afa..0000000000 --- a/boost/math/tools/remez.hpp +++ /dev/null @@ -1,667 +0,0 @@ -// (C) Copyright John Maddock 2006. -// Use, modification and distribution are subject to the -// Boost Software License, Version 1.0. (See accompanying file -// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) - -#ifndef BOOST_MATH_TOOLS_REMEZ_HPP -#define BOOST_MATH_TOOLS_REMEZ_HPP - -#ifdef _MSC_VER -#pragma once -#endif - -#include <boost/math/tools/solve.hpp> -#include <boost/math/tools/minima.hpp> -#include <boost/math/tools/roots.hpp> -#include <boost/math/tools/polynomial.hpp> -#include <boost/function/function1.hpp> -#include <boost/scoped_array.hpp> -#include <boost/math/constants/constants.hpp> -#include <boost/math/policies/policy.hpp> - -namespace boost{ namespace math{ namespace tools{ - -namespace detail{ - -// -// The error function: the difference between F(x) and -// the current approximation. This is the function -// for which we must find the extema. -// -template <class T> -struct remez_error_function -{ - typedef boost::function1<T, T const &> function_type; -public: - remez_error_function( - function_type f_, - const polynomial<T>& n, - const polynomial<T>& d, - bool rel_err) - : f(f_), numerator(n), denominator(d), rel_error(rel_err) {} - - T operator()(const T& z)const - { - T y = f(z); - T abs = y - (numerator.evaluate(z) / denominator.evaluate(z)); - T err; - if(rel_error) - { - if(y != 0) - err = abs / fabs(y); - else if(0 == abs) - { - // we must be at a root, or it's not recoverable: - BOOST_ASSERT(0 == abs); - err = 0; - } - else - { - // We have a divide by zero! - // Lets assume that f(x) is zero as a result of - // internal cancellation error that occurs as a result - // of shifting a root at point z to the origin so that - // the approximation can be "pinned" to pass through - // the origin: in that case it really - // won't matter what our approximation calculates here - // as long as it's a small number, return the absolute error: - err = abs; - } - } - else - err = abs; - return err; - } -private: - function_type f; - polynomial<T> numerator; - polynomial<T> denominator; - bool rel_error; -}; -// -// This function adapts the error function so that it's minima -// are the extema of the error function. We can find the minima -// with standard techniques. -// -template <class T> -struct remez_max_error_function -{ - remez_max_error_function(const remez_error_function<T>& f) - : func(f) {} - - T operator()(const T& x) - { - BOOST_MATH_STD_USING - return -fabs(func(x)); - } -private: - remez_error_function<T> func; -}; - -} // detail - -template <class T> -class remez_minimax -{ -public: - typedef boost::function1<T, T const &> function_type; - typedef boost::numeric::ublas::vector<T> vector_type; - typedef boost::numeric::ublas::matrix<T> matrix_type; - - remez_minimax(function_type f, unsigned oN, unsigned oD, T a, T b, bool pin = true, bool rel_err = false, int sk = 0, int bits = 0); - remez_minimax(function_type f, unsigned oN, unsigned oD, T a, T b, bool pin, bool rel_err, int sk, int bits, const vector_type& points); - - void reset(unsigned oN, unsigned oD, T a, T b, bool pin = true, bool rel_err = false, int sk = 0, int bits = 0); - void reset(unsigned oN, unsigned oD, T a, T b, bool pin, bool rel_err, int sk, int bits, const vector_type& points); - - void set_brake(int b) - { - BOOST_ASSERT(b < 100); - BOOST_ASSERT(b >= 0); - m_brake = b; - } - - T iterate(); - - polynomial<T> denominator()const; - polynomial<T> numerator()const; - - vector_type const& chebyshev_points()const - { - return control_points; - } - - vector_type const& zero_points()const - { - return zeros; - } - - T error_term()const - { - return solution[solution.size() - 1]; - } - T max_error()const - { - return m_max_error; - } - T max_change()const - { - return m_max_change; - } - void rotate() - { - --orderN; - ++orderD; - } - void rescale(T a, T b) - { - T scale = (b - a) / (max - min); - for(unsigned i = 0; i < control_points.size(); ++i) - { - control_points[i] = (control_points[i] - min) * scale + a; - } - min = a; - max = b; - } -private: - - void init_chebyshev(); - - function_type func; // The function to approximate. - vector_type control_points; // Current control points to be used for the next iteration. - vector_type solution; // Solution from the last iteration contains all unknowns including the error term. - vector_type zeros; // Location of points of zero error from last iteration, plus the two end points. - vector_type maxima; // Location of maxima of the error function, actually contains the control points used for the last iteration. - T m_max_error; // Maximum error found in last approximation. - T m_max_change; // Maximum change in location of control points after last iteration. - unsigned orderN; // Order of the numerator polynomial. - unsigned orderD; // Order of the denominator polynomial. - T min, max; // End points of the range to optimise over. - bool rel_error; // If true optimise for relative not absolute error. - bool pinned; // If true the approximation is "pinned" to go through the origin. - unsigned unknowns; // Total number of unknowns. - int m_precision; // Number of bits precision to which the zeros and maxima are found. - T m_max_change_history[2]; // Past history of changes to control points. - int m_brake; // amount to break by in percentage points. - int m_skew; // amount to skew starting points by in percentage points: -100-100 -}; - -#ifndef BRAKE -#define BRAKE 0 -#endif -#ifndef SKEW -#define SKEW 0 -#endif - -template <class T> -void remez_minimax<T>::init_chebyshev() -{ - BOOST_MATH_STD_USING - // - // Fill in the zeros: - // - unsigned terms = pinned ? orderD + orderN : orderD + orderN + 1; - - for(unsigned i = 0; i < terms; ++i) - { - T cheb = cos((2 * terms - 1 - 2 * i) * constants::pi<T>() / (2 * terms)); - cheb += 1; - cheb /= 2; - if(m_skew != 0) - { - T p = static_cast<T>(200 + m_skew) / 200; - cheb = pow(cheb, p); - } - cheb *= (max - min); - cheb += min; - zeros[i+1] = cheb; - } - zeros[0] = min; - zeros[unknowns] = max; - // perform a regular interpolation fit: - matrix_type A(terms, terms); - vector_type b(terms); - // fill in the y values: - for(unsigned i = 0; i < b.size(); ++i) - { - b[i] = func(zeros[i+1]); - } - // fill in powers of x evaluated at each of the control points: - unsigned offsetN = pinned ? 0 : 1; - unsigned offsetD = offsetN + orderN; - unsigned maxorder = (std::max)(orderN, orderD); - for(unsigned i = 0; i < b.size(); ++i) - { - T x0 = zeros[i+1]; - T x = x0; - if(!pinned) - A(i, 0) = 1; - for(unsigned j = 0; j < maxorder; ++j) - { - if(j < orderN) - A(i, j + offsetN) = x; - if(j < orderD) - { - A(i, j + offsetD) = -x * b[i]; - } - x *= x0; - } - } - // - // Now go ahead and solve the expression to get our solution: - // - vector_type l_solution = boost::math::tools::solve(A, b); - // need to add a "fake" error term: - l_solution.resize(unknowns); - l_solution[unknowns-1] = 0; - solution = l_solution; - // - // Now find all the extrema of the error function: - // - detail::remez_error_function<T> Err(func, this->numerator(), this->denominator(), rel_error); - detail::remez_max_error_function<T> Ex(Err); - m_max_error = 0; - int max_err_location = 0; - for(unsigned i = 0; i < unknowns; ++i) - { - std::pair<T, T> r = brent_find_minima(Ex, zeros[i], zeros[i+1], m_precision); - maxima[i] = r.first; - T rel_err = fabs(r.second); - if(rel_err > m_max_error) - { - m_max_error = fabs(r.second); - max_err_location = i; - } - } - control_points = maxima; -} - -template <class T> -void remez_minimax<T>::reset( - unsigned oN, - unsigned oD, - T a, - T b, - bool pin, - bool rel_err, - int sk, - int bits) -{ - control_points = vector_type(oN + oD + (pin ? 1 : 2)); - solution = control_points; - zeros = vector_type(oN + oD + (pin ? 2 : 3)); - maxima = control_points; - orderN = oN; - orderD = oD; - rel_error = rel_err; - pinned = pin; - m_skew = sk; - min = a; - max = b; - m_max_error = 0; - unknowns = orderN + orderD + (pinned ? 1 : 2); - // guess our initial control points: - control_points[0] = min; - control_points[unknowns - 1] = max; - T interval = (max - min) / (unknowns - 1); - T spot = min + interval; - for(unsigned i = 1; i < control_points.size(); ++i) - { - control_points[i] = spot; - spot += interval; - } - solution[unknowns - 1] = 0; - m_max_error = 0; - if(bits == 0) - { - // don't bother about more than float precision: - m_precision = (std::min)(24, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2); - } - else - { - // can't be more accurate than half the bits of T: - m_precision = (std::min)(bits, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2); - } - m_max_change_history[0] = m_max_change_history[1] = 1; - init_chebyshev(); - // do one iteration whatever: - //iterate(); -} - -template <class T> -inline remez_minimax<T>::remez_minimax( - typename remez_minimax<T>::function_type f, - unsigned oN, - unsigned oD, - T a, - T b, - bool pin, - bool rel_err, - int sk, - int bits) - : func(f) -{ - m_brake = 0; - reset(oN, oD, a, b, pin, rel_err, sk, bits); -} - -template <class T> -void remez_minimax<T>::reset( - unsigned oN, - unsigned oD, - T a, - T b, - bool pin, - bool rel_err, - int sk, - int bits, - const vector_type& points) -{ - control_points = vector_type(oN + oD + (pin ? 1 : 2)); - solution = control_points; - zeros = vector_type(oN + oD + (pin ? 2 : 3)); - maxima = control_points; - orderN = oN; - orderD = oD; - rel_error = rel_err; - pinned = pin; - m_skew = sk; - min = a; - max = b; - m_max_error = 0; - unknowns = orderN + orderD + (pinned ? 1 : 2); - control_points = points; - solution[unknowns - 1] = 0; - m_max_error = 0; - if(bits == 0) - { - // don't bother about more than float precision: - m_precision = (std::min)(24, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2); - } - else - { - // can't be more accurate than half the bits of T: - m_precision = (std::min)(bits, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2); - } - m_max_change_history[0] = m_max_change_history[1] = 1; - // do one iteration whatever: - //iterate(); -} - -template <class T> -inline remez_minimax<T>::remez_minimax( - typename remez_minimax<T>::function_type f, - unsigned oN, - unsigned oD, - T a, - T b, - bool pin, - bool rel_err, - int sk, - int bits, - const vector_type& points) - : func(f) -{ - m_brake = 0; - reset(oN, oD, a, b, pin, rel_err, sk, bits, points); -} - -template <class T> -T remez_minimax<T>::iterate() -{ - BOOST_MATH_STD_USING - matrix_type A(unknowns, unknowns); - vector_type b(unknowns); - - // fill in evaluation of f(x) at each of the control points: - for(unsigned i = 0; i < b.size(); ++i) - { - // take care that none of our control points are at the origin: - if(pinned && (control_points[i] == 0)) - { - if(i) - control_points[i] = control_points[i-1] / 3; - else - control_points[i] = control_points[i+1] / 3; - } - b[i] = func(control_points[i]); - } - - T err_err; - unsigned convergence_count = 0; - do{ - // fill in powers of x evaluated at each of the control points: - int sign = 1; - unsigned offsetN = pinned ? 0 : 1; - unsigned offsetD = offsetN + orderN; - unsigned maxorder = (std::max)(orderN, orderD); - T Elast = solution[unknowns - 1]; - - for(unsigned i = 0; i < b.size(); ++i) - { - T x0 = control_points[i]; - T x = x0; - if(!pinned) - A(i, 0) = 1; - for(unsigned j = 0; j < maxorder; ++j) - { - if(j < orderN) - A(i, j + offsetN) = x; - if(j < orderD) - { - T mult = rel_error ? (b[i] - sign * fabs(b[i]) * Elast): (b[i] - sign * Elast); - A(i, j + offsetD) = -x * mult; - } - x *= x0; - } - // The last variable to be solved for is the error term, - // sign changes with each control point: - T E = rel_error ? sign * fabs(b[i]) : sign; - A(i, unknowns - 1) = E; - sign = -sign; - } - - #ifdef BOOST_MATH_INSTRUMENT - for(unsigned i = 0; i < b.size(); ++i) - std::cout << b[i] << " "; - std::cout << "\n\n"; - for(unsigned i = 0; i < b.size(); ++i) - { - for(unsigned j = 0; j < b.size(); ++ j) - std::cout << A(i, j) << " "; - std::cout << "\n"; - } - std::cout << std::endl; - #endif - // - // Now go ahead and solve the expression to get our solution: - // - solution = boost::math::tools::solve(A, b); - - err_err = (Elast != 0) ? fabs((fabs(solution[unknowns-1]) - fabs(Elast)) / fabs(Elast)) : 1; - }while(orderD && (convergence_count++ < 80) && (err_err > 0.001)); - - // - // Perform a sanity check to verify that the solution to the equations - // is not so much in error as to be useless. The matrix inversion can - // be very close to singular, so this can be a real problem. - // - vector_type sanity = prod(A, solution); - for(unsigned i = 0; i < b.size(); ++i) - { - T err = fabs((b[i] - sanity[i]) / fabs(b[i])); - if(err > sqrt(epsilon<T>())) - { - std::cerr << "Sanity check failed: more than half the digits in the found solution are in error." << std::endl; - } - } - - // - // Next comes another sanity check, we want to verify that all the control - // points do actually alternate in sign, in practice we may have - // additional roots in the error function that cause this to fail. - // Failure here is always fatal: even though this code attempts to correct - // the problem it usually only postpones the inevitable. - // - polynomial<T> num, denom; - num = this->numerator(); - denom = this->denominator(); - T e1 = b[0] - num.evaluate(control_points[0]) / denom.evaluate(control_points[0]); -#ifdef BOOST_MATH_INSTRUMENT - std::cout << e1; -#endif - for(unsigned i = 1; i < b.size(); ++i) - { - T e2 = b[i] - num.evaluate(control_points[i]) / denom.evaluate(control_points[i]); -#ifdef BOOST_MATH_INSTRUMENT - std::cout << " " << e2; -#endif - if(e2 * e1 > 0) - { - std::cerr << std::flush << "Basic sanity check failed: Error term does not alternate in sign, non-recoverable error may follow..." << std::endl; - T perturbation = 0.05; - do{ - T point = control_points[i] * (1 - perturbation) + control_points[i-1] * perturbation; - e2 = func(point) - num.evaluate(point) / denom.evaluate(point); - if(e2 * e1 < 0) - { - control_points[i] = point; - break; - } - perturbation += 0.05; - }while(perturbation < 0.8); - - if((e2 * e1 > 0) && (i + 1 < b.size())) - { - perturbation = 0.05; - do{ - T point = control_points[i] * (1 - perturbation) + control_points[i+1] * perturbation; - e2 = func(point) - num.evaluate(point) / denom.evaluate(point); - if(e2 * e1 < 0) - { - control_points[i] = point; - break; - } - perturbation += 0.05; - }while(perturbation < 0.8); - } - - } - e1 = e2; - } - -#ifdef BOOST_MATH_INSTRUMENT - for(unsigned i = 0; i < solution.size(); ++i) - std::cout << solution[i] << " "; - std::cout << std::endl << this->numerator() << std::endl; - std::cout << this->denominator() << std::endl; - std::cout << std::endl; -#endif - - // - // The next step is to find all the intervals in which our maxima - // lie: - // - detail::remez_error_function<T> Err(func, this->numerator(), this->denominator(), rel_error); - zeros[0] = min; - zeros[unknowns] = max; - for(unsigned i = 1; i < control_points.size(); ++i) - { - eps_tolerance<T> tol(m_precision); - boost::uintmax_t max_iter = 1000; - std::pair<T, T> p = toms748_solve( - Err, - control_points[i-1], - control_points[i], - tol, - max_iter); - zeros[i] = (p.first + p.second) / 2; - //zeros[i] = bisect(Err, control_points[i-1], control_points[i], m_precision); - } - // - // Now find all the extrema of the error function: - // - detail::remez_max_error_function<T> Ex(Err); - m_max_error = 0; - int max_err_location = 0; - for(unsigned i = 0; i < unknowns; ++i) - { - std::pair<T, T> r = brent_find_minima(Ex, zeros[i], zeros[i+1], m_precision); - maxima[i] = r.first; - T rel_err = fabs(r.second); - if(rel_err > m_max_error) - { - m_max_error = fabs(r.second); - max_err_location = i; - } - } - // - // Almost done now! we just need to set our control points - // to the extrema, and calculate how much each point has changed - // (this will be our termination condition): - // - swap(control_points, maxima); - m_max_change = 0; - int max_change_location = 0; - for(unsigned i = 0; i < unknowns; ++i) - { - control_points[i] = (control_points[i] * (100 - m_brake) + maxima[i] * m_brake) / 100; - T change = fabs((control_points[i] - maxima[i]) / control_points[i]); -#if 0 - if(change > m_max_change_history[1]) - { - // divergence!!! try capping the change: - std::cerr << "Possible divergent step, change will be capped!!" << std::endl; - change = m_max_change_history[1]; - if(control_points[i] < maxima[i]) - control_points[i] = maxima[i] - change * maxima[i]; - else - control_points[i] = maxima[i] + change * maxima[i]; - } -#endif - if(change > m_max_change) - { - m_max_change = change; - max_change_location = i; - } - } - // - // store max change information: - // - m_max_change_history[0] = m_max_change_history[1]; - m_max_change_history[1] = fabs(m_max_change); - - return m_max_change; -} - -template <class T> -polynomial<T> remez_minimax<T>::numerator()const -{ - boost::scoped_array<T> a(new T[orderN + 1]); - if(pinned) - a[0] = 0; - unsigned terms = pinned ? orderN : orderN + 1; - for(unsigned i = 0; i < terms; ++i) - a[pinned ? i+1 : i] = solution[i]; - return boost::math::tools::polynomial<T>(&a[0], orderN); -} - -template <class T> -polynomial<T> remez_minimax<T>::denominator()const -{ - unsigned terms = orderD + 1; - unsigned offsetD = pinned ? orderN : (orderN + 1); - boost::scoped_array<T> a(new T[terms]); - a[0] = 1; - for(unsigned i = 0; i < orderD; ++i) - a[i+1] = solution[i + offsetD]; - return boost::math::tools::polynomial<T>(&a[0], orderD); -} - - -}}} // namespaces - -#endif // BOOST_MATH_TOOLS_REMEZ_HPP - - - |