ceph/src/boost/libs/math/include_private/boost/math/tools/remez.hpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_TOOLS_REMEZ_HPP
   7 #define BOOST_MATH_TOOLS_REMEZ_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/tools/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
  22 namespace boost{ namespace math{ namespace tools{
  23
  24 namespace detail{
  25
  26 //
  27 // The error function: the difference between F(x) and
  28 // the current approximation.  This is the function
  29 // for which we must find the extema.
  30 //
  31 template <class T>
  32 struct remez_error_function
  33 {
  34    typedef boost::function1<T, T const &> function_type;
  35 public:
  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)const
  44    {
  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_ASSERT(0 == abs);
  56             err = 0;
  57          }
  58          else
  59          {
  60             // We have a divide by zero!
  61             // Lets assume that f(x) is zero as a result of
  62             // internal cancellation error that occurs as a result
  63             // of shifting a root at point z to the origin so that
  64             // the approximation can be "pinned" to pass through
  65             // the origin: in that case it really
  66             // won't matter what our approximation calculates here
  67             // as long as it's a small number, return the absolute error:
  68             err = abs;
  69          }
  70       }
  71       else
  72          err = abs;
  73       return err;
  74    }
  75 private:
  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 minima
  83 // are the extema of the error function.  We can find the minima
  84 // with standard techniques.
  85 //
  86 template <class T>
  87 struct remez_max_error_function
  88 {
  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_USING
  95       return -fabs(func(x));
  96    }
  97 private:
  98    remez_error_function<T> func;
  99 };
 100
 101 } // detail
 102
 103 template <class T>
 104 class remez_minimax
 105 {
 106 public:
 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_ASSERT(b < 100);
 120       BOOST_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()const
 130    {
 131       return control_points;
 132    }
 133
 134    vector_type const& zero_points()const
 135    {
 136       return zeros;
 137    }
 138
 139    T error_term()const
 140    {
 141       return solution[solution.size() - 1];
 142    }
 143    T max_error()const
 144    {
 145       return m_max_error;
 146    }
 147    T max_change()const
 148    {
 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    }
 166 private:
 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-100
 187 };
 188
 189 #ifndef BRAKE
 190 #define BRAKE 0
 191 #endif
 192 #ifndef SKEW
 193 #define SKEW 0
 194 #endif
 195
 196 template <class T>
 197 void remez_minimax<T>::init_chebyshev()
 198 {
 199    BOOST_MATH_STD_USING
 200    //
 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
 279 template <class T>
 280 void 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    else
 321    {
 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
 331 template <class T>
 332 inline 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
 348 template <class T>
 349 void 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    else
 382    {
 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
 391 template <class T>
 392 inline 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
 409 template <class T>
 410 T remez_minimax<T>::iterate()
 411 {
 412    BOOST_MATH_STD_USING
 413    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          else
 425             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_INSTRUMENT
 465       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    #endif
 476       //
 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 equations
 486    // is not so much in error as to be useless.  The matrix inversion can
 487    // 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 control
 501    // 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 correct
 504    // 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_INSTRUMENT
 511    std::cout << e1;
 512 #endif
 513    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_INSTRUMENT
 517       std::cout << " " << e2;
 518 #endif
 519       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_INSTRUMENT
 554    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 #endif
 560
 561    //
 562    // The next step is to find all the intervals in which our maxima
 563    // 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       boost::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 points
 600    // to the extrema, and calculate how much each point has changed
 601    // (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 0
 611       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          else
 619             control_points[i] = maxima[i] + change * maxima[i];
 620       }
 621 #endif
 622       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
 637 template <class T>
 638 polynomial<T> remez_minimax<T>::numerator()const
 639 {
 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
 649 template <class T>
 650 polynomial<T> remez_minimax<T>::denominator()const
 651 {
 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 }}} // namespaces
 663
 664 #endif // BOOST_MATH_TOOLS_REMEZ_HPP
 665
 666
 667