668 lines · plain
1// (C) Copyright John Maddock 2006.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_TOOLS_REMEZ_HPP7#define BOOST_MATH_TOOLS_REMEZ_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include "solve.hpp"14#include <boost/math/tools/minima.hpp>15#include <boost/math/tools/roots.hpp>16#include <boost/math/tools/polynomial.hpp>17#include <boost/function/function1.hpp>18#include <boost/scoped_array.hpp>19#include <boost/math/constants/constants.hpp>20#include <boost/math/policies/policy.hpp>21 22namespace boost{ namespace math{ namespace tools{23 24namespace detail{25 26//27// The error function: the difference between F(x) and28// the current approximation. This is the function29// for which we must find the extrema.30//31template <class T>32struct remez_error_function33{34 typedef boost::function1<T, T const &> function_type;35public:36 remez_error_function(37 function_type f_, 38 const polynomial<T>& n, 39 const polynomial<T>& d, 40 bool rel_err)41 : f(f_), numerator(n), denominator(d), rel_error(rel_err) {}42 43 T operator()(const T& z)const44 {45 T y = f(z);46 T abs = y - (numerator.evaluate(z) / denominator.evaluate(z));47 T err;48 if(rel_error)49 {50 if(y != 0)51 err = abs / fabs(y);52 else if(0 == abs)53 {54 // we must be at a root, or it's not recoverable:55 BOOST_MATH_ASSERT(0 == abs);56 err = 0;57 }58 else59 {60 // We have a divide by zero!61 // Lets assume that f(x) is zero as a result of62 // internal cancellation error that occurs as a result63 // of shifting a root at point z to the origin so that64 // the approximation can be "pinned" to pass through65 // the origin: in that case it really66 // won't matter what our approximation calculates here67 // as long as it's a small number, return the absolute error:68 err = abs;69 }70 }71 else72 err = abs;73 return err;74 }75private:76 function_type f;77 polynomial<T> numerator;78 polynomial<T> denominator;79 bool rel_error;80};81//82// This function adapts the error function so that it's minima83// are the extrema of the error function. We can find the minima84// with standard techniques.85//86template <class T>87struct remez_max_error_function88{89 remez_max_error_function(const remez_error_function<T>& f)90 : func(f) {}91 92 T operator()(const T& x)93 {94 BOOST_MATH_STD_USING95 return -fabs(func(x));96 }97private:98 remez_error_function<T> func;99};100 101} // detail102 103template <class T>104class remez_minimax105{106public:107 typedef boost::function1<T, T const &> function_type;108 typedef boost::numeric::ublas::vector<T> vector_type;109 typedef boost::numeric::ublas::matrix<T> matrix_type;110 111 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);112 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);113 114 void reset(unsigned oN, unsigned oD, T a, T b, bool pin = true, bool rel_err = false, int sk = 0, int bits = 0);115 void reset(unsigned oN, unsigned oD, T a, T b, bool pin, bool rel_err, int sk, int bits, const vector_type& points);116 117 void set_brake(int b)118 {119 BOOST_MATH_ASSERT(b < 100);120 BOOST_MATH_ASSERT(b >= 0);121 m_brake = b;122 }123 124 T iterate();125 126 polynomial<T> denominator()const;127 polynomial<T> numerator()const;128 129 vector_type const& chebyshev_points()const130 {131 return control_points;132 }133 134 vector_type const& zero_points()const135 {136 return zeros;137 }138 139 T error_term()const140 {141 return solution[solution.size() - 1];142 }143 T max_error()const144 {145 return m_max_error;146 }147 T max_change()const148 {149 return m_max_change;150 }151 void rotate()152 {153 --orderN;154 ++orderD;155 }156 void rescale(T a, T b)157 {158 T scale = (b - a) / (max - min);159 for(unsigned i = 0; i < control_points.size(); ++i)160 {161 control_points[i] = (control_points[i] - min) * scale + a;162 }163 min = a;164 max = b;165 }166private:167 168 void init_chebyshev();169 170 function_type func; // The function to approximate.171 vector_type control_points; // Current control points to be used for the next iteration.172 vector_type solution; // Solution from the last iteration contains all unknowns including the error term.173 vector_type zeros; // Location of points of zero error from last iteration, plus the two end points.174 vector_type maxima; // Location of maxima of the error function, actually contains the control points used for the last iteration.175 T m_max_error; // Maximum error found in last approximation.176 T m_max_change; // Maximum change in location of control points after last iteration.177 unsigned orderN; // Order of the numerator polynomial.178 unsigned orderD; // Order of the denominator polynomial.179 T min, max; // End points of the range to optimise over.180 bool rel_error; // If true optimise for relative not absolute error.181 bool pinned; // If true the approximation is "pinned" to go through the origin.182 unsigned unknowns; // Total number of unknowns.183 int m_precision; // Number of bits precision to which the zeros and maxima are found.184 T m_max_change_history[2]; // Past history of changes to control points.185 int m_brake; // amount to break by in percentage points.186 int m_skew; // amount to skew starting points by in percentage points: -100-100187};188 189#ifndef BRAKE190#define BRAKE 0191#endif192#ifndef SKEW193#define SKEW 0194#endif195 196template <class T>197void remez_minimax<T>::init_chebyshev()198{199 BOOST_MATH_STD_USING200 //201 // Fill in the zeros:202 //203 unsigned terms = pinned ? orderD + orderN : orderD + orderN + 1;204 205 for(unsigned i = 0; i < terms; ++i)206 {207 T cheb = cos((2 * terms - 1 - 2 * i) * constants::pi<T>() / (2 * terms));208 cheb += 1;209 cheb /= 2;210 if(m_skew != 0)211 {212 T p = static_cast<T>(200 + m_skew) / 200;213 cheb = pow(cheb, p);214 }215 cheb *= (max - min);216 cheb += min;217 zeros[i+1] = cheb;218 }219 zeros[0] = min;220 zeros[unknowns] = max;221 // perform a regular interpolation fit:222 matrix_type A(terms, terms);223 vector_type b(terms);224 // fill in the y values:225 for(unsigned i = 0; i < b.size(); ++i)226 {227 b[i] = func(zeros[i+1]);228 }229 // fill in powers of x evaluated at each of the control points:230 unsigned offsetN = pinned ? 0 : 1;231 unsigned offsetD = offsetN + orderN;232 unsigned maxorder = (std::max)(orderN, orderD);233 for(unsigned i = 0; i < b.size(); ++i)234 {235 T x0 = zeros[i+1];236 T x = x0;237 if(!pinned)238 A(i, 0) = 1;239 for(unsigned j = 0; j < maxorder; ++j)240 {241 if(j < orderN)242 A(i, j + offsetN) = x;243 if(j < orderD)244 {245 A(i, j + offsetD) = -x * b[i];246 }247 x *= x0;248 }249 }250 //251 // Now go ahead and solve the expression to get our solution:252 //253 vector_type l_solution = boost::math::tools::solve(A, b);254 // need to add a "fake" error term:255 l_solution.resize(unknowns);256 l_solution[unknowns-1] = 0;257 solution = l_solution;258 //259 // Now find all the extrema of the error function:260 //261 detail::remez_error_function<T> Err(func, this->numerator(), this->denominator(), rel_error);262 detail::remez_max_error_function<T> Ex(Err);263 m_max_error = 0;264 //int max_err_location = 0;265 for(unsigned i = 0; i < unknowns; ++i)266 {267 std::pair<T, T> r = brent_find_minima(Ex, zeros[i], zeros[i+1], m_precision);268 maxima[i] = r.first;269 T rel_err = fabs(r.second);270 if(rel_err > m_max_error)271 {272 m_max_error = fabs(r.second);273 //max_err_location = i;274 }275 }276 control_points = maxima;277}278 279template <class T>280void remez_minimax<T>::reset(281 unsigned oN, 282 unsigned oD, 283 T a, 284 T b, 285 bool pin, 286 bool rel_err, 287 int sk,288 int bits)289{290 control_points = vector_type(oN + oD + (pin ? 1 : 2));291 solution = control_points;292 zeros = vector_type(oN + oD + (pin ? 2 : 3));293 maxima = control_points;294 orderN = oN;295 orderD = oD;296 rel_error = rel_err;297 pinned = pin;298 m_skew = sk;299 min = a;300 max = b;301 m_max_error = 0;302 unknowns = orderN + orderD + (pinned ? 1 : 2);303 // guess our initial control points:304 control_points[0] = min;305 control_points[unknowns - 1] = max;306 T interval = (max - min) / (unknowns - 1);307 T spot = min + interval;308 for(unsigned i = 1; i < control_points.size(); ++i)309 {310 control_points[i] = spot;311 spot += interval;312 }313 solution[unknowns - 1] = 0;314 m_max_error = 0;315 if(bits == 0)316 {317 // don't bother about more than float precision:318 m_precision = (std::min)(24, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);319 }320 else321 {322 // can't be more accurate than half the bits of T:323 m_precision = (std::min)(bits, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);324 }325 m_max_change_history[0] = m_max_change_history[1] = 1;326 init_chebyshev();327 // do one iteration whatever:328 //iterate();329}330 331template <class T>332inline remez_minimax<T>::remez_minimax(333 typename remez_minimax<T>::function_type f, 334 unsigned oN, 335 unsigned oD, 336 T a, 337 T b, 338 bool pin, 339 bool rel_err, 340 int sk,341 int bits)342 : func(f) 343{344 m_brake = 0;345 reset(oN, oD, a, b, pin, rel_err, sk, bits);346}347 348template <class T>349void remez_minimax<T>::reset(350 unsigned oN, 351 unsigned oD, 352 T a, 353 T b, 354 bool pin, 355 bool rel_err, 356 int sk,357 int bits,358 const vector_type& points)359{360 control_points = vector_type(oN + oD + (pin ? 1 : 2));361 solution = control_points;362 zeros = vector_type(oN + oD + (pin ? 2 : 3));363 maxima = control_points;364 orderN = oN;365 orderD = oD;366 rel_error = rel_err;367 pinned = pin;368 m_skew = sk;369 min = a;370 max = b;371 m_max_error = 0;372 unknowns = orderN + orderD + (pinned ? 1 : 2);373 control_points = points;374 solution[unknowns - 1] = 0;375 m_max_error = 0;376 if(bits == 0)377 {378 // don't bother about more than float precision:379 m_precision = (std::min)(24, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);380 }381 else382 {383 // can't be more accurate than half the bits of T:384 m_precision = (std::min)(bits, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);385 }386 m_max_change_history[0] = m_max_change_history[1] = 1;387 // do one iteration whatever:388 //iterate();389}390 391template <class T>392inline remez_minimax<T>::remez_minimax(393 typename remez_minimax<T>::function_type f, 394 unsigned oN, 395 unsigned oD, 396 T a, 397 T b, 398 bool pin, 399 bool rel_err, 400 int sk,401 int bits,402 const vector_type& points)403 : func(f)404{405 m_brake = 0;406 reset(oN, oD, a, b, pin, rel_err, sk, bits, points);407}408 409template <class T>410T remez_minimax<T>::iterate()411{412 BOOST_MATH_STD_USING413 matrix_type A(unknowns, unknowns);414 vector_type b(unknowns);415 416 // fill in evaluation of f(x) at each of the control points:417 for(unsigned i = 0; i < b.size(); ++i)418 {419 // take care that none of our control points are at the origin:420 if(pinned && (control_points[i] == 0))421 {422 if(i)423 control_points[i] = control_points[i-1] / 3;424 else425 control_points[i] = control_points[i+1] / 3;426 }427 b[i] = func(control_points[i]);428 }429 430 T err_err;431 unsigned convergence_count = 0;432 do{433 // fill in powers of x evaluated at each of the control points:434 int sign = 1;435 unsigned offsetN = pinned ? 0 : 1;436 unsigned offsetD = offsetN + orderN;437 unsigned maxorder = (std::max)(orderN, orderD);438 T Elast = solution[unknowns - 1];439 440 for(unsigned i = 0; i < b.size(); ++i)441 {442 T x0 = control_points[i];443 T x = x0;444 if(!pinned)445 A(i, 0) = 1;446 for(unsigned j = 0; j < maxorder; ++j)447 {448 if(j < orderN)449 A(i, j + offsetN) = x;450 if(j < orderD)451 {452 T mult = rel_error ? T(b[i] - sign * fabs(b[i]) * Elast): T(b[i] - sign * Elast);453 A(i, j + offsetD) = -x * mult;454 }455 x *= x0;456 }457 // The last variable to be solved for is the error term, 458 // sign changes with each control point:459 T E = rel_error ? T(sign * fabs(b[i])) : T(sign);460 A(i, unknowns - 1) = E;461 sign = -sign;462 }463 464 #ifdef BOOST_MATH_INSTRUMENT465 for(unsigned i = 0; i < b.size(); ++i)466 std::cout << b[i] << " ";467 std::cout << "\n\n";468 for(unsigned i = 0; i < b.size(); ++i)469 {470 for(unsigned j = 0; j < b.size(); ++ j)471 std::cout << A(i, j) << " ";472 std::cout << "\n";473 }474 std::cout << std::endl;475 #endif476 //477 // Now go ahead and solve the expression to get our solution:478 //479 solution = boost::math::tools::solve(A, b);480 481 err_err = (Elast != 0) ? T(fabs((fabs(solution[unknowns-1]) - fabs(Elast)) / fabs(Elast))) : T(1);482 }while(orderD && (convergence_count++ < 80) && (err_err > 0.001));483 484 //485 // Perform a sanity check to verify that the solution to the equations486 // is not so much in error as to be useless. The matrix inversion can487 // be very close to singular, so this can be a real problem.488 //489 vector_type sanity = prod(A, solution);490 for(unsigned i = 0; i < b.size(); ++i)491 {492 T err = fabs((b[i] - sanity[i]) / fabs(b[i]));493 if(err > sqrt(epsilon<T>()))494 {495 std::cerr << "Sanity check failed: more than half the digits in the found solution are in error." << std::endl;496 }497 }498 499 //500 // Next comes another sanity check, we want to verify that all the control501 // points do actually alternate in sign, in practice we may have 502 // additional roots in the error function that cause this to fail.503 // Failure here is always fatal: even though this code attempts to correct504 // the problem it usually only postpones the inevitable.505 //506 polynomial<T> num, denom;507 num = this->numerator();508 denom = this->denominator();509 T e1 = b[0] - num.evaluate(control_points[0]) / denom.evaluate(control_points[0]);510#ifdef BOOST_MATH_INSTRUMENT511 std::cout << e1;512#endif513 for(unsigned i = 1; i < b.size(); ++i)514 {515 T e2 = b[i] - num.evaluate(control_points[i]) / denom.evaluate(control_points[i]);516#ifdef BOOST_MATH_INSTRUMENT517 std::cout << " " << e2;518#endif519 if(e2 * e1 > 0)520 {521 std::cerr << std::flush << "Basic sanity check failed: Error term does not alternate in sign, non-recoverable error may follow..." << std::endl;522 T perturbation = 0.05;523 do{524 T point = control_points[i] * (1 - perturbation) + control_points[i-1] * perturbation;525 e2 = func(point) - num.evaluate(point) / denom.evaluate(point);526 if(e2 * e1 < 0)527 {528 control_points[i] = point;529 break;530 }531 perturbation += 0.05;532 }while(perturbation < 0.8);533 534 if((e2 * e1 > 0) && (i + 1 < b.size()))535 {536 perturbation = 0.05;537 do{538 T point = control_points[i] * (1 - perturbation) + control_points[i+1] * perturbation;539 e2 = func(point) - num.evaluate(point) / denom.evaluate(point);540 if(e2 * e1 < 0)541 {542 control_points[i] = point;543 break;544 }545 perturbation += 0.05;546 }while(perturbation < 0.8);547 }548 549 }550 e1 = e2;551 }552 553#ifdef BOOST_MATH_INSTRUMENT554 for(unsigned i = 0; i < solution.size(); ++i)555 std::cout << solution[i] << " ";556 std::cout << std::endl << this->numerator() << std::endl;557 std::cout << this->denominator() << std::endl;558 std::cout << std::endl;559#endif560 561 //562 // The next step is to find all the intervals in which our maxima563 // lie:564 //565 detail::remez_error_function<T> Err(func, this->numerator(), this->denominator(), rel_error);566 zeros[0] = min;567 zeros[unknowns] = max;568 for(unsigned i = 1; i < control_points.size(); ++i)569 {570 eps_tolerance<T> tol(m_precision);571 std::uintmax_t max_iter = 1000;572 std::pair<T, T> p = toms748_solve(573 Err, 574 control_points[i-1], 575 control_points[i], 576 tol, 577 max_iter);578 zeros[i] = (p.first + p.second) / 2;579 //zeros[i] = bisect(Err, control_points[i-1], control_points[i], m_precision);580 }581 //582 // Now find all the extrema of the error function:583 //584 detail::remez_max_error_function<T> Ex(Err);585 m_max_error = 0;586 //int max_err_location = 0;587 for(unsigned i = 0; i < unknowns; ++i)588 {589 std::pair<T, T> r = brent_find_minima(Ex, zeros[i], zeros[i+1], m_precision);590 maxima[i] = r.first;591 T rel_err = fabs(r.second);592 if(rel_err > m_max_error)593 {594 m_max_error = fabs(r.second);595 //max_err_location = i;596 }597 }598 //599 // Almost done now! we just need to set our control points600 // to the extrema, and calculate how much each point has changed601 // (this will be our termination condition):602 //603 swap(control_points, maxima);604 m_max_change = 0;605 //int max_change_location = 0;606 for(unsigned i = 0; i < unknowns; ++i)607 {608 control_points[i] = (control_points[i] * (100 - m_brake) + maxima[i] * m_brake) / 100;609 T change = fabs((control_points[i] - maxima[i]) / control_points[i]);610#if 0611 if(change > m_max_change_history[1])612 {613 // divergence!!! try capping the change:614 std::cerr << "Possible divergent step, change will be capped!!" << std::endl;615 change = m_max_change_history[1];616 if(control_points[i] < maxima[i])617 control_points[i] = maxima[i] - change * maxima[i];618 else619 control_points[i] = maxima[i] + change * maxima[i];620 }621#endif622 if(change > m_max_change)623 {624 m_max_change = change;625 //max_change_location = i;626 }627 }628 //629 // store max change information:630 //631 m_max_change_history[0] = m_max_change_history[1];632 m_max_change_history[1] = fabs(m_max_change);633 634 return m_max_change;635}636 637template <class T>638polynomial<T> remez_minimax<T>::numerator()const639{640 boost::scoped_array<T> a(new T[orderN + 1]);641 if(pinned)642 a[0] = 0;643 unsigned terms = pinned ? orderN : orderN + 1;644 for(unsigned i = 0; i < terms; ++i)645 a[pinned ? i+1 : i] = solution[i];646 return boost::math::tools::polynomial<T>(&a[0], orderN);647}648 649template <class T>650polynomial<T> remez_minimax<T>::denominator()const651{652 unsigned terms = orderD + 1;653 unsigned offsetD = pinned ? orderN : (orderN + 1);654 boost::scoped_array<T> a(new T[terms]);655 a[0] = 1;656 for(unsigned i = 0; i < orderD; ++i)657 a[i+1] = solution[i + offsetD];658 return boost::math::tools::polynomial<T>(&a[0], orderD);659}660 661 662}}} // namespaces663 664#endif // BOOST_MATH_TOOLS_REMEZ_HPP665 666 667 668