main.cpp

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

//  (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)

#define BOOST_UBLAS_TYPE_CHECK_EPSILON (type_traits<real_type>::type_sqrt (boost::math::tools::epsilon <real_type>()))
#define BOOST_UBLAS_TYPE_CHECK_MIN (type_traits<real_type>::type_sqrt ( boost::math::tools::min_value<real_type>()))
#define BOOST_UBLAS_NDEBUG

#include <boost/math/bindings/rr.hpp>
namespace std{
using boost::math::ntl::pow;
} // workaround for spirit parser.
#include <boost/math/tools/remez.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/math/special_functions/binomial.hpp>
#include <boost/spirit/core.hpp>
#include <boost/spirit/actor.hpp>
#include <boost/lexical_cast.hpp>
#include <iostream>
#include <iomanip>
#include <string>
#include <boost/test/included/test_exec_monitor.hpp> // for test_main

extern boost::math::ntl::RR f(const boost::math::ntl::RR& x, int variant);
extern void show_extra(
   const boost::math::tools::polynomial<boost::math::ntl::RR>& n, 
   const boost::math::tools::polynomial<boost::math::ntl::RR>& d, 
   const boost::math::ntl::RR& x_offset, 
   const boost::math::ntl::RR& y_offset, 
   int variant);

using namespace boost::spirit;

boost::math::ntl::RR a(0), b(1);   // range to optimise over
bool rel_error(true);
bool pin(false);
int orderN(3);
int orderD(1);
int target_precision = boost::math::tools::digits<long double>();
int working_precision = target_precision * 2;
bool started(false);
int variant(0);
int skew(0);
int brake(50);
boost::math::ntl::RR x_offset(0), y_offset(0), x_scale(1);
bool auto_offset_y;

boost::shared_ptr<boost::math::tools::remez_minimax<boost::math::ntl::RR> > p_remez;

boost::math::ntl::RR the_function(const boost::math::ntl::RR& val)
{
   return f(x_scale * (val + x_offset), variant) + y_offset;
}

void step_some(unsigned count)
{
   try{
      NTL::RR::SetPrecision(working_precision);
      if(!started)
      {
         //
         // If we have an automatic y-offset calculate it now:
         //
         if(auto_offset_y)
         {
            boost::math::ntl::RR fa, fb, fm;
            fa = f(x_scale * (a + x_offset), variant);
            fb = f(x_scale * (b + x_offset), variant);
            fm = f(x_scale * ((a+b)/2 + x_offset), variant);
            y_offset = -(fa + fb + fm) / 3;
            NTL::RR::SetOutputPrecision(5);
            std::cout << "Setting auto-y-offset to " << y_offset << std::endl;
         }
         //
         // Truncate offsets to float precision:
         //
         x_offset = NTL::RoundToPrecision(x_offset.value(), 20);
         y_offset = NTL::RoundToPrecision(y_offset.value(), 20);
         //
         // Construct new Remez state machine:
         //
         p_remez.reset(new boost::math::tools::remez_minimax<boost::math::ntl::RR>(
            &the_function, 
            orderN, orderD, 
            a, b, 
            pin, 
            rel_error, 
            skew, 
            working_precision));
         std::cout << "Max error in interpolated form: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(p_remez->max_error()) << std::endl;
         //
         // Signal that we've started:
         //
         started = true;
      }
      unsigned i;
      for(i = 0; i < count; ++i)
      {
         std::cout << "Stepping..." << std::endl;
         p_remez->set_brake(brake);
         boost::math::ntl::RR r = p_remez->iterate();
         NTL::RR::SetOutputPrecision(3);
         std::cout 
            << "Maximum Deviation Found:                     " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(p_remez->max_error()) << std::endl
            << "Expected Error Term:                         " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(p_remez->error_term()) << std::endl
            << "Maximum Relative Change in Control Points:   " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(r) << std::endl;
      }
   }
   catch(const std::exception& e)
   {
      std::cout << "Step failed with exception: " << e.what() << std::endl;
   }
}

void step(const char*, const char*)
{
   step_some(1);
}

void show(const char*, const char*)
{
   NTL::RR::SetPrecision(working_precision);
   if(started)
   {
      boost::math::tools::polynomial<boost::math::ntl::RR> n = p_remez->numerator();
      boost::math::tools::polynomial<boost::math::ntl::RR> d = p_remez->denominator();
      std::vector<boost::math::ntl::RR> cn = n.chebyshev();
      std::vector<boost::math::ntl::RR> cd = d.chebyshev();
      int prec = 2 + (target_precision * 3010LL)/10000;
      std::cout << std::scientific << std::setprecision(prec);
      NTL::RR::SetOutputPrecision(prec);
      boost::numeric::ublas::vector<boost::math::ntl::RR> v = p_remez->zero_points();
      
      std::cout << "  Zeros = {\n";
      unsigned i;
      for(i = 0; i < v.size(); ++i)
      {
         std::cout << "    " << v[i] << std::endl;
      }
      std::cout << "  }\n";

      v = p_remez->chebyshev_points();
      std::cout << "  Chebeshev Control Points = {\n";
      for(i = 0; i < v.size(); ++i)
      {
         std::cout << "    " << v[i] << std::endl;
      }
      std::cout << "  }\n";

      std::cout << "X offset: " << x_offset << std::endl;
      std::cout << "X scale:  " << x_scale << std::endl;
      std::cout << "Y offset: " << y_offset << std::endl;

      std::cout << "P = {";
      for(i = 0; i < n.size(); ++i)
      {
         std::cout << "    " << n[i] << "L," << std::endl;
      }
      std::cout << "  }\n";

      std::cout << "Q = {";
      for(i = 0; i < d.size(); ++i)
      {
         std::cout << "    " << d[i] << "L," << std::endl;
      }
      std::cout << "  }\n";

      std::cout << "CP = {";
      for(i = 0; i < cn.size(); ++i)
      {
         std::cout << "    " << cn[i] << "L," << std::endl;
      }
      std::cout << "  }\n";

      std::cout << "CQ = {";
      for(i = 0; i < cd.size(); ++i)
      {
         std::cout << "    " << cd[i] << "L," << std::endl;
      }
      std::cout << "  }\n";

      show_extra(n, d, x_offset, y_offset, variant);
   }
   else
   {
      std::cerr << "Nothing to display" << std::endl;
   }
}

void do_graph(unsigned points)
{
   NTL::RR::SetPrecision(working_precision);
   boost::math::ntl::RR step = (b - a) / (points - 1);
   boost::math::ntl::RR x = a;
   while(points > 1)
   {
      NTL::RR::SetOutputPrecision(10);
      std::cout << std::setprecision(10) << std::setw(30) << std::left 
         << boost::lexical_cast<std::string>(x) << the_function(x) << std::endl;
      --points;
      x += step;
   }
   std::cout << std::setprecision(10) << std::setw(30) << std::left 
      << boost::lexical_cast<std::string>(b) << the_function(b) << std::endl;
}

void graph(const char*, const char*)
{
   do_graph(3);
}

template <class T>
void do_test(T, const char* name)
{
   boost::math::ntl::RR::SetPrecision(working_precision);
   if(started)
   {
      //
      // We want to test the approximation at fixed precision:
      // either float, double or long double.  Begin by getting the
      // polynomials:
      //
      boost::math::tools::polynomial<T> n, d;
      boost::math::tools::polynomial<boost::math::ntl::RR> nr, dr;
      nr = p_remez->numerator();
      dr = p_remez->denominator();
      n = nr;
      d = dr;

      std::vector<boost::math::ntl::RR> cn1, cd1;
      cn1 = nr.chebyshev();
      cd1 = dr.chebyshev();
      std::vector<T> cn, cd;
      for(unsigned i = 0; i < cn1.size(); ++i)
      {
         cn.push_back(boost::math::tools::real_cast<T>(cn1[i]));
      }
      for(unsigned i = 0; i < cd1.size(); ++i)
      {
         cd.push_back(boost::math::tools::real_cast<T>(cd1[i]));
      }
      //
      // We'll test at the Chebeshev control points which is where
      // (in theory) the largest deviation should occur.  For good
      // measure we'll test at the zeros as well:
      //
      boost::numeric::ublas::vector<boost::math::ntl::RR> 
         zeros(p_remez->zero_points()),
         cheb(p_remez->chebyshev_points());

      boost::math::ntl::RR max_error(0), cheb_max_error(0);

      //
      // Do the tests at the zeros:
      //
      std::cout << "Starting tests at " << name << " precision...\n";
      std::cout << "Absissa        Error (Poly)   Error (Cheb)\n";
      for(unsigned i = 0; i < zeros.size(); ++i)
      {
         boost::math::ntl::RR true_result = the_function(zeros[i]);
         T absissa = boost::math::tools::real_cast<T>(zeros[i]);
         boost::math::ntl::RR test_result = n.evaluate(absissa) / d.evaluate(absissa);
         boost::math::ntl::RR cheb_result = boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa);
         boost::math::ntl::RR err, cheb_err;
         if(rel_error)
         {
            err = boost::math::tools::relative_error(test_result, true_result);
            cheb_err = boost::math::tools::relative_error(cheb_result, true_result);
         }
         else
         {
            err = fabs(test_result - true_result);
            cheb_err = fabs(cheb_result - true_result);
         }
         if(err > max_error)
            max_error = err;
         if(cheb_err > cheb_max_error)
            cheb_max_error = cheb_err;
         std::cout << std::setprecision(6) << std::setw(15) << std::left << absissa
            << std::setw(15) << std::left << boost::math::tools::real_cast<T>(err) << boost::math::tools::real_cast<T>(cheb_err) << std::endl;
      }
      //
      // Do the tests at the Chebeshev control points:
      //
      for(unsigned i = 0; i < cheb.size(); ++i)
      {
         boost::math::ntl::RR true_result = the_function(cheb[i]);
         T absissa = boost::math::tools::real_cast<T>(cheb[i]);
         boost::math::ntl::RR test_result = n.evaluate(absissa) / d.evaluate(absissa);
         boost::math::ntl::RR cheb_result = boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa);
         boost::math::ntl::RR err, cheb_err;
         if(rel_error)
         {
            err = boost::math::tools::relative_error(test_result, true_result);
            cheb_err = boost::math::tools::relative_error(cheb_result, true_result);
         }
         else
         {
            err = fabs(test_result - true_result);
            cheb_err = fabs(cheb_result - true_result);
         }
         if(err > max_error)
            max_error = err;
         std::cout << std::setprecision(6) << std::setw(15) << std::left << absissa