main.cpp

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

            << std::setw(15) << std::left << boost::math::tools::real_cast<T>(err) << 
            boost::math::tools::real_cast<T>(cheb_err) << std::endl;
      }
      std::string msg = "Max Error found at ";
      msg += name;
      msg += " precision = ";
      msg.append(62 - 17 - msg.size(), ' ');
      std::cout << msg << std::setprecision(6) << "Poly: " << std::setw(20) << std::left
         << boost::math::tools::real_cast<T>(max_error) << "Cheb: " << boost::math::tools::real_cast<T>(cheb_max_error) << std::endl;
   }
   else
   {
      std::cout << "Nothing to test: try converging an approximation first!!!" << std::endl;
   }
}

void test_float(const char*, const char*)
{
   do_test(float(0), "float");
}

void test_double(const char*, const char*)
{
   do_test(double(0), "double");
}

void test_long(const char*, const char*)
{
   do_test((long double)(0), "long double");
}

void test_all(const char*, const char*)
{
   do_test(float(0), "float");
   do_test(double(0), "double");
   do_test((long double)(0), "long double");
}

template <class T>
void do_test_n(T, const char* name, unsigned count)
{
   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]));
      }

      boost::math::ntl::RR max_error(0), max_cheb_error(0);
      boost::math::ntl::RR step = (b - a) / count;

      //
      // Do the tests at the zeros:
      //
      std::cout << "Starting tests at " << name << " precision...\n";
      std::cout << "Absissa        Error (poly)   Error (Cheb)\n";
      for(boost::math::ntl::RR x = a; x <= b; x += step)
      {
         boost::math::ntl::RR true_result = the_function(x);
         T absissa = boost::math::tools::real_cast<T>(x);
         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 > max_cheb_error)
            max_cheb_error = cheb_err;
         std::cout << std::setprecision(6) << std::setw(15) << std::left << boost::math::tools::real_cast<double>(absissa)
            << (test_result < true_result ? "-" : "") << std::setw(20) << std::left 
            << boost::math::tools::real_cast<double>(err) 
            << boost::math::tools::real_cast<double>(cheb_err) << std::endl;
      }
      std::string msg = "Max Error found at ";
      msg += name;
      msg += " precision = ";
      msg.append(62 - 17 - msg.size(), ' ');
      std::cout << msg << "Poly: " << std::setprecision(6) 
         << std::setw(15) << std::left 
         << boost::math::tools::real_cast<T>(max_error) 
         << "Cheb: " << boost::math::tools::real_cast<T>(max_cheb_error) << std::endl;
   }
   else
   {
      std::cout << "Nothing to test: try converging an approximation first!!!" << std::endl;
   }
}

void test_n(unsigned n)
{
   do_test_n(boost::math::ntl::RR(), "boost::math::ntl::RR", n);
}

void test_float_n(unsigned n)
{
   do_test_n(float(0), "float", n);
}

void test_double_n(unsigned n)
{
   do_test_n(double(0), "double", n);
}

void test_long_n(unsigned n)
{
   do_test_n((long double)(0), "long double", n);
}

void rotate(const char*, const char*)
{
   if(p_remez)
   {
      p_remez->rotate();
   }
   else
   {
      std::cerr << "Nothing to rotate" << std::endl;
   }
}

void rescale(const char*, const char*)
{
   if(p_remez)
   {
      p_remez->rescale(a, b);
   }
   else
   {
      std::cerr << "Nothing to rescale" << std::endl;
   }
}

void graph_poly(const char*, const char*)
{
   int i = 50;
   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<boost::math::ntl::RR> n, d;
      n = p_remez->numerator();
      d = p_remez->denominator();

      boost::math::ntl::RR max_error(0);
      boost::math::ntl::RR step = (b - a) / i;

      std::cout << "Evaluating Numerator...\n";
      boost::math::ntl::RR val;
      for(val = a; val <= b; val += step)
         std::cout << n.evaluate(val) << std::endl;
      std::cout << "Evaluating Denominator...\n";
      for(val = a; val <= b; val += step)
         std::cout << d.evaluate(val) << std::endl;
   }
   else
   {
      std::cout << "Nothing to test: try converging an approximation first!!!" << std::endl;
   }
}

int test_main(int, char* [])
{
   std::string line;
   real_parser<long double/*boost::math::ntl::RR*/ > const rr_p;
   while(std::getline(std::cin, line))
   {
      if(parse(line.c_str(), str_p("quit"), space_p).full)
         return 0;
      if(false == parse(line.c_str(), 
         (

            str_p("range")[assign_a(started, false)] && real_p[assign_a(a)] && real_p[assign_a(b)]
      ||
            str_p("relative")[assign_a(started, false)][assign_a(rel_error, true)]
      ||
            str_p("absolute")[assign_a(started, false)][assign_a(rel_error, false)]
      ||
            str_p("pin")[assign_a(started, false)] && str_p("true")[assign_a(pin, true)]
      ||
            str_p("pin")[assign_a(started, false)] && str_p("false")[assign_a(pin, false)]
      ||
            str_p("pin")[assign_a(started, false)] && str_p("1")[assign_a(pin, true)]
      ||
            str_p("pin")[assign_a(started, false)] && str_p("0")[assign_a(pin, false)]
      ||
            str_p("pin")[assign_a(started, false)][assign_a(pin, true)]
      ||
            str_p("order")[assign_a(started, false)] && uint_p[assign_a(orderN)] && uint_p[assign_a(orderD)]
      ||
            str_p("order")[assign_a(started, false)] && uint_p[assign_a(orderN)]
      ||
            str_p("target-precision") && uint_p[assign_a(target_precision)]
      ||
            str_p("working-precision")[assign_a(started, false)] && uint_p[assign_a(working_precision)]
      ||
            str_p("variant")[assign_a(started, false)] && int_p[assign_a(variant)]
      ||
            str_p("skew")[assign_a(started, false)] && int_p[assign_a(skew)]
      ||
            str_p("brake") && int_p[assign_a(brake)]
      ||
            str_p("step") && int_p[&step_some]
      ||
            str_p("step")[&step]
      ||
            str_p("poly")[&graph_poly]
      ||
            str_p("info")[&show]
      ||
            str_p("graph") && uint_p[&do_graph]
      ||
            str_p("graph")[&graph]
      ||
            str_p("x-offset") && real_p[assign_a(x_offset)]
      ||
            str_p("x-scale") && real_p[assign_a(x_scale)]
      ||
            str_p("y-offset") && str_p("auto")[assign_a(auto_offset_y, true)]
      ||
            str_p("y-offset") && real_p[assign_a(y_offset)][assign_a(auto_offset_y, false)]
      ||
            str_p("test") && str_p("float") && uint_p[&test_float_n]
      ||
            str_p("test") && str_p("float")[&test_float]
      ||
            str_p("test") && str_p("double") && uint_p[&test_double_n]
      ||
            str_p("test") && str_p("double")[&test_double]
      ||
            str_p("test") && str_p("long") && uint_p[&test_long_n]
      ||
            str_p("test") && str_p("long")[&test_long]
      ||
            str_p("test") && str_p("all")[&test_all]
      ||
            str_p("test") && uint_p[&test_n]
      ||
            str_p("rotate")[&rotate]
      ||
            str_p("rescale") && real_p[assign_a(a)] && real_p[assign_a(b)] && epsilon_p[&rescale]

         ), space_p).full)
      {
         std::cout << "Unable to parse directive: \"" << line << "\"" << std::endl;
      }
      else
      {
         std::cout << "Variant              = " << variant << std::endl;
         std::cout << "range                = [" << a << "," << b << "]" << std::endl;
         std::cout << "Relative Error       = " << rel_error << std::endl;
         std::cout << "Pin to Origin        = " << pin << std::endl;
         std::cout << "Order (Num/Denom)    = " << orderN << "/" << orderD << std::endl;
         std::cout << "Target Precision     = " << target_precision << std::endl;
         std::cout << "Working Precision    = " << working_precision << std::endl;
         std::cout << "Skew                 = " << skew << std::endl;
         std::cout << "Brake                = " << brake << std::endl;
         std::cout << "X Offset             = " << x_offset << std::endl;
         std::cout << "X scale              = " << x_scale << std::endl;
         std::cout << "Y Offset             = ";
         if(auto_offset_y)
            std::cout << "Auto (";
         std::cout << y_offset;
         if(auto_offset_y)
            std::cout << ")";
         std::cout << std::endl;
     }
   }
   return 0;
}