//  (C) Copyright John Maddock 2007.
//  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)

#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <sstream>

int max_order = 20;
const char* path_prefix = "..\\..\\..\\boost\\math\\tools\\detail\\polynomial_";
const char* path_prefix2 = "..\\..\\..\\boost\\math\\tools\\detail\\rational_";

const char* copyright_string = 
"//  (C) Copyright John Maddock 2007.\n"
"//  Use, modification and distribution are subject to the\n"
"//  Boost Software License, Version 1.0. (See accompanying file\n"
"//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)\n"
"//\n"
"//  This file is machine generated, do not edit by hand\n\n";


void print_polynomials(int max_order)
{
   for(int i = 2; i <= max_order; ++i)
   {
      std::stringstream filename;
      filename << path_prefix << "horner1_" << i << ".hpp";
      std::ofstream ofs(filename.str().c_str());
      if(ofs.bad())
         break;
      //
      // Output the boilerplate at the top of the header:
      //
      ofs << copyright_string <<
         "// Polynomial evaluation using Horners rule\n"
         "#ifndef BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n"
         "#define BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n\n"
         "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<0>*)\n"
         "{\n"
         "   return static_cast<V>(0);\n"
         "}\n"
         "\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<1>*)\n"
         "{\n"
         "   return static_cast<V>(a[0]);\n"
         "}\n\n";

      for(int order = 2; order <= i; ++order)
      {
         ofs << 
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<" << order << ">*)\n"
         "{\n"
         "   return static_cast<V>(";
         
         for(int bracket = 2; bracket < order; ++bracket)
            ofs << "(";
         ofs << "a[" << order - 1 << "] * x + a[" << order - 2 << "]" ;
         for(int item = order - 3; item >= 0; --item)
         {
            ofs << ") * x + a[" << item << "]";
         }
         
         ofs << ");\n"
         "}\n\n";
      }
      //
      // And finally the boilerplate at the end of the header:
      //
      ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";

      filename.str("");
      filename << path_prefix << "horner2_" << i << ".hpp";
      ofs.close();
      ofs.open(filename.str().c_str());
      if(ofs.bad())
         break;
      //
      // Output the boilerplate at the top of the header:
      //
      ofs << copyright_string <<
         "// Polynomial evaluation using second order Horners rule\n"
         "#ifndef BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n"
         "#define BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n\n"
         "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<0>*)\n"
         "{\n"
         "   return static_cast<V>(0);\n"
         "}\n"
         "\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<1>*)\n"
         "{\n"
         "   return static_cast<V>(a[0]);\n"
         "}\n\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<2>*)\n"
         "{\n"
         "   return static_cast<V>(a[1] * x + a[0]);\n"
         "}\n\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<3>*)\n"
         "{\n"
         "   return static_cast<V>((a[2] * x + a[1]) * x + a[0]);\n"
         "}\n\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<4>*)\n"
         "{\n"
         "   return static_cast<V>(((a[3] * x + a[2]) * x + a[1]) * x + a[0]);\n"
         "}\n\n";

      for(int order = 5; order <= i; ++order)
      {
         ofs << 
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<" << order << ">*)\n"
         "{\n"
         "   V x2 = x * x;\n"
         "   return static_cast<V>(";

         if(order & 1)
         {
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "a[" << order - 1 << "] * x2 + a[" << order - 3 << "]" ;
            for(int item = order - 5; item >= 0; item -= 2)
            {
               ofs << ") * x2 + a[" << item << "]";
            }
            ofs << " + ";
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "a[" << order - 2 << "] * x2 + a[" << order - 4 << "]" ;
            for(int item = order - 6; item >= 0; item -= 2)
            {
               ofs << ") * x2 + a[" << item << "]";
            }
            ofs << ") * x";
         }
         else
         {
            for(int bracket = 0; bracket < (order - 1) / 2; ++bracket)
               ofs << "(";
            ofs << "a[" << order - 1 << "] * x2 + a[" << order - 3 << "]" ;
            for(int item = order - 5; item >= 0; item -= 2)
            {
               ofs << ") * x2 + a[" << item << "]";
            }
            ofs << ") * x + ";
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "a[" << order - 2 << "] * x2 + a[" << order - 4 << "]" ;
            for(int item = order - 6; item >= 0; item -= 2)
            {
               ofs << ") * x2 + a[" << item << "]";
            }
         }
         ofs << ");\n"
            "}\n\n";
      }
      //
      // And finally the boilerplate at the end of the header:
      //
      ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";


      filename.str("");
      filename << path_prefix << "horner3_" << i << ".hpp";
      ofs.close();
      ofs.open(filename.str().c_str());
      if(ofs.bad())
         break;
      //
      // Output the boilerplate at the top of the header:
      //
      ofs << copyright_string <<
         "// Unrolled polynomial evaluation using second order Horners rule\n"
         "#ifndef BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n"
         "#define BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n\n"
         "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<0>*)\n"
         "{\n"
         "   return static_cast<V>(0);\n"
         "}\n"
         "\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<1>*)\n"
         "{\n"
         "   return static_cast<V>(a[0]);\n"
         "}\n\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<2>*)\n"
         "{\n"
         "   return static_cast<V>(a[1] * x + a[0]);\n"
         "}\n\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<3>*)\n"
         "{\n"
         "   return static_cast<V>((a[2] * x + a[1]) * x + a[0]);\n"
         "}\n\n"
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<4>*)\n"
         "{\n"
         "   return static_cast<V>(((a[3] * x + a[2]) * x + a[1]) * x + a[0]);\n"
         "}\n\n";

      for(int order = 5; order <= i; ++order)
      {
         ofs << 
         "template <class T, class V>\n"
         "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<" << order << ">*)\n"
         "{\n"
         "   V x2 = x * x;\n"
         "   V t[2];\n";

         if(order & 1)
         {
            ofs << "   t[0] = static_cast<V>(a[" << order - 1 << "] * x2 + a[" << order - 3 << "]);\n" ;
            ofs << "   t[1] = static_cast<V>(a[" << order - 2 << "] * x2 + a[" << order - 4 << "]);\n" ;
            for(int item = order - 5; item >= 0; item -= 2)
            {
               ofs << "   t[0] *= x2;\n";
               if(item - 1 >= 0)
                  ofs << "   t[1] *= x2;\n";
               ofs << "   t[0] += static_cast<V>(a[" << item << "]);\n";
               if(item - 1 >= 0)
                  ofs << "   t[1] += static_cast<V>(a[" << item - 1 << "]);\n";
            }
            ofs << 
               "   t[1] *= x;\n"
               "   return t[0] + t[1];\n";
         }
         else
         {
            ofs << "   t[0] = a[" << order - 1 << "] * x2 + a[" << order - 3 << "];\n" ;
            ofs << "   t[1] = a[" << order - 2 << "] * x2 + a[" << order - 4 << "];\n" ;
            for(int item = order - 5; item >= 0; item -= 2)
            {
               ofs << "   t[0] *= x2;\n";
               if(item - 1 >= 0)
                  ofs << "   t[1] *= x2;\n";
               ofs << "   t[0] += static_cast<V>(a[" << item << "]);\n";
               if(item - 1 >= 0)
                  ofs <<  "   t[1] += static_cast<V>(a[" << item - 1 << "]);\n";
            }
            ofs << "   t[0] *= x;\n";
            ofs << "   return t[0] + t[1];\n";
         }
         ofs << "}\n\n";
      }
      //
      // And finally the boilerplate at the end of the header:
      //
      ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";
   }
}

void print_rationals(int max_order)
{
   for(int i = 2; i <= max_order; ++i)
   {
      std::stringstream filename;
      filename << path_prefix2 << "horner1_" << i << ".hpp";
      std::ofstream ofs(filename.str().c_str());
      if(ofs.bad())
         break;
      //
      // Output the boilerplate at the top of the header:
      //
      ofs << copyright_string <<
         "// Polynomial evaluation using Horners rule\n"
         "#ifndef BOOST_MATH_TOOLS_POLY_RAT_" << i << "_HPP\n"
         "#define BOOST_MATH_TOOLS_POLY_RAT_" << i << "_HPP\n\n"
         "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T*, const U*, const V&, const mpl::int_<0>*)\n"
         "{\n"
         "   return static_cast<V>(0);\n"
         "}\n"
         "\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const mpl::int_<1>*)\n"
         "{\n"
         "   return static_cast<V>(a[0]) / static_cast<V>(b[0]);\n"
         "}\n\n";

      for(int order = 2; order <= i; ++order)
      {
         ofs << 
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<" << order << ">*)\n"
         "{\n"
         "   if(x <= 1)\n"
         "     return static_cast<V>((";
         
         for(int bracket = 2; bracket < order; ++bracket)
            ofs << "(";
         ofs << "a[" << order - 1 << "] * x + a[" << order - 2 << "]" ;
         for(int item = order - 3; item >= 0; --item)
         {
            ofs << ") * x + a[" << item << "]";
         }

         ofs << ") / (";
         for(int bracket = 2; bracket < order; ++bracket)
            ofs << "(";
         ofs << "b[" << order - 1 << "] * x + b[" << order - 2 << "]" ;
         for(int item = order - 3; item >= 0; --item)
         {
            ofs << ") * x + b[" << item << "]";
         }
         
         ofs << "));\n   else\n   {\n      V z = 1 / x;\n      return static_cast<V>((";

         for(int bracket = order - 1; bracket > 1; --bracket)
            ofs << "(";
         ofs << "a[" << 0 << "] * z + a[" << 1 << "]" ;
         for(int item = 2; item <= order - 1; ++item)
         {
            ofs << ") * z + a[" << item << "]";
         }

         ofs << ") / (";
         for(int bracket = 2; bracket < order; ++bracket)
            ofs << "(";
         ofs << "b[" << 0 << "] * z + b[" << 1 << "]" ;
         for(int item = 2; item <= order - 1; ++item)
         {
            ofs << ") * z + b[" << item << "]";
         }
         
         ofs << "));\n   }\n";
         
         ofs << "}\n\n";
      }
      //
      // And finally the boilerplate at the end of the header:
      //
      ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";

      filename.str("");
      filename << path_prefix2 << "horner2_" << i << ".hpp";
      ofs.close();
      ofs.open(filename.str().c_str());
      if(ofs.bad())
         break;
      //
      // Output the boilerplate at the top of the header:
      //
      ofs << copyright_string <<
         "// Polynomial evaluation using second order Horners rule\n"
         "#ifndef BOOST_MATH_TOOLS_RAT_EVAL_" << i << "_HPP\n"
         "#define BOOST_MATH_TOOLS_RAT_EVAL_" << i << "_HPP\n\n"
         "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const mpl::int_<0>*)\n"
         "{\n"
         "   return static_cast<V>(0);\n"
         "}\n"
         "\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const mpl::int_<1>*)\n"
         "{\n"
         "   return static_cast<V>(a[0]) / static_cast<V>(b[0]);\n"
         "}\n\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<2>*)\n"
         "{\n"
         "   return static_cast<V>((a[1] * x + a[0]) / (b[1] * x + b[0]));\n"
         "}\n\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<3>*)\n"
         "{\n"
         "   return static_cast<V>(((a[2] * x + a[1]) * x + a[0]) / ((b[2] * x + b[1]) * x + b[0]));\n"
         "}\n\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<4>*)\n"
         "{\n"
         "   return static_cast<V>((((a[3] * x + a[2]) * x + a[1]) * x + a[0]) / (((b[3] * x + b[2]) * x + b[1]) * x + b[0]));\n"
         "}\n\n";

      for(int order = 5; order <= i; ++order)
      {
         ofs << 
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<" << order << ">*)\n"
         "{\n"
         "   if(x <= 1)\n   {\n"
         "      V x2 = x * x;\n"
         "      return static_cast<V>((";

         if(order & 1)
         {
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "a[" << order - 1 << "] * x2 + a[" << order - 3 << "]" ;
            for(int item = order - 5; item >= 0; item -= 2)
            {
               ofs << ") * x2 + a[" << item << "]";
            }
            ofs << " + ";
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "a[" << order - 2 << "] * x2 + a[" << order - 4 << "]" ;
            for(int item = order - 6; item >= 0; item -= 2)
            {
               ofs << ") * x2 + a[" << item << "]";
            }
            ofs << ") * x";
         }
         else
         {
            for(int bracket = 0; bracket < (order - 1) / 2; ++bracket)
               ofs << "(";
            ofs << "a[" << order - 1 << "] * x2 + a[" << order - 3 << "]" ;
            for(int item = order - 5; item >= 0; item -= 2)
            {
               ofs << ") * x2 + a[" << item << "]";
            }
            ofs << ") * x + ";
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "a[" << order - 2 << "] * x2 + a[" << order - 4 << "]" ;
            for(int item = order - 6; item >= 0; item -= 2)
            {
               ofs << ") * x2 + a[" << item << "]";
            }
         }
         ofs << ") / (";
         if(order & 1)
         {
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "b[" << order - 1 << "] * x2 + b[" << order - 3 << "]" ;
            for(int item = order - 5; item >= 0; item -= 2)
            {
               ofs << ") * x2 + b[" << item << "]";
            }
            ofs << " + ";
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "b[" << order - 2 << "] * x2 + b[" << order - 4 << "]" ;
            for(int item = order - 6; item >= 0; item -= 2)
            {
               ofs << ") * x2 + b[" << item << "]";
            }
            ofs << ") * x";
         }
         else
         {
            for(int bracket = 0; bracket < (order - 1) / 2; ++bracket)
               ofs << "(";
            ofs << "b[" << order - 1 << "] * x2 + b[" << order - 3 << "]" ;
            for(int item = order - 5; item >= 0; item -= 2)
            {
               ofs << ") * x2 + b[" << item << "]";
            }
            ofs << ") * x + ";
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "b[" << order - 2 << "] * x2 + b[" << order - 4 << "]" ;
            for(int item = order - 6; item >= 0; item -= 2)
            {
               ofs << ") * x2 + b[" << item << "]";
            }
         }

         ofs << "));\n   }\n   else\n   {\n      V z = 1 / x;\n      V z2 = 1 / (x * x);\n      return static_cast<V>((";

         if(order & 1)
         {
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "a[" << 0 << "] * z2 + a[" << 2 << "]" ;
            for(int item = 4; item < order; item += 2)
            {
               ofs << ") * z2 + a[" << item << "]";
            }
            ofs << " + ";
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "a[" << 1 << "] * z2 + a[" << 3 << "]" ;
            for(int item = 5; item < order; item += 2)
            {
               ofs << ") * z2 + a[" << item << "]";
            }
            ofs << ") * z";
         }
         else
         {
            for(int bracket = 0; bracket < (order - 1) / 2; ++bracket)
               ofs << "(";
            ofs << "a[" << 0 << "] * z2 + a[" << 2 << "]" ;
            for(int item = 4; item < order; item += 2)
            {
               ofs << ") * z2 + a[" << item << "]";
            }
            ofs << ") * z + ";
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "a[" << 1 << "] * z2 + a[" << 3 << "]" ;
            for(int item = 5; item < order; item += 2)
            {
               ofs << ") * z2 + a[" << item << "]";
            }
         }

         ofs << ") / (";

         if(order & 1)
         {
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "b[" << 0 << "] * z2 + b[" << 2 << "]" ;
            for(int item = 4; item < order; item += 2)
            {
               ofs << ") * z2 + b[" << item << "]";
            }
            ofs << " + ";
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "b[" << 1 << "] * z2 + b[" << 3 << "]" ;
            for(int item = 5; item < order; item += 2)
            {
               ofs << ") * z2 + b[" << item << "]";
            }
            ofs << ") * z";
         }
         else
         {
            for(int bracket = 0; bracket < (order - 1) / 2; ++bracket)
               ofs << "(";
            ofs << "b[" << 0 << "] * z2 + b[" << 2 << "]" ;
            for(int item = 4; item < order; item += 2)
            {
               ofs << ") * z2 + b[" << item << "]";
            }
            ofs << ") * z + ";
            for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
               ofs << "(";
            ofs << "b[" << 1 << "] * z2 + b[" << 3 << "]" ;
            for(int item = 5; item < order; item += 2)
            {
               ofs << ") * z2 + b[" << item << "]";
            }
         }
         ofs << "));\n   }\n";

         ofs << "}\n\n";
      }
      //
      // And finally the boilerplate at the end of the header:
      //
      ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";


      filename.str("");
      filename << path_prefix2 << "horner3_" << i << ".hpp";
      ofs.close();
      ofs.open(filename.str().c_str());
      if(ofs.bad())
         break;
      //
      // Output the boilerplate at the top of the header:
      //
      ofs << copyright_string <<
         "// Polynomial evaluation using second order Horners rule\n"
         "#ifndef BOOST_MATH_TOOLS_RAT_EVAL_" << i << "_HPP\n"
         "#define BOOST_MATH_TOOLS_RAT_EVAL_" << i << "_HPP\n\n"
         "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const mpl::int_<0>*)\n"
         "{\n"
         "   return static_cast<V>(0);\n"
         "}\n"
         "\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const mpl::int_<1>*)\n"
         "{\n"
         "   return static_cast<V>(a[0]) / static_cast<V>(b[0]);\n"
         "}\n\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<2>*)\n"
         "{\n"
         "   return static_cast<V>((a[1] * x + a[0]) / (b[1] * x + b[0]));\n"
         "}\n\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<3>*)\n"
         "{\n"
         "   return static_cast<V>(((a[2] * x + a[1]) * x + a[0]) / ((b[2] * x + b[1]) * x + b[0]));\n"
         "}\n\n"
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<4>*)\n"
         "{\n"
         "   return static_cast<V>((((a[3] * x + a[2]) * x + a[1]) * x + a[0]) / (((b[3] * x + b[2]) * x + b[1]) * x + b[0]));\n"
         "}\n\n";

      for(int order = 5; order <= i; ++order)
      {
         ofs << 
         "template <class T, class U, class V>\n"
         "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<" << order << ">*)\n"
         "{\n"
         "   if(x <= 1)\n   {\n"
         "      V x2 = x * x;\n"
         "      V t[4];\n";

         if(order & 1)
         {
            ofs << "      t[0] = a[" << order - 1 << "] * x2 + a[" << order - 3 << "];\n" ;
            ofs << "      t[1] = a[" << order - 2 << "] * x2 + a[" << order - 4 << "];\n" ;
            ofs << "      t[2] = b[" << order - 1 << "] * x2 + b[" << order - 3 << "];\n" ;
            ofs << "      t[3] = b[" << order - 2 << "] * x2 + b[" << order - 4 << "];\n" ;
            for(int item = order - 5; item >= 0; item -= 2)
            {
               ofs << "      t[0] *= x2;\n";
               if(item - 1 >= 0)
                  ofs << "      t[1] *= x2;\n";
               ofs << "      t[2] *= x2;\n";
               if(item - 1 >= 0)
                  ofs << "      t[3] *= x2;\n";
               ofs << "      t[0] += static_cast<V>(a[" << item << "]);\n";
               if(item - 1 >= 0)
                  ofs << "      t[1] += static_cast<V>(a[" << item - 1 << "]);\n";
               ofs << "      t[2] += static_cast<V>(b[" << item << "]);\n";
               if(item - 1 >= 0)
                  ofs << "      t[3] += static_cast<V>(b[" << item - 1 << "]);\n";
            }
            ofs << "      t[1] *= x;\n";
            ofs << "      t[3] *= x;\n";
         }
         else
         {
            ofs << "      t[0] = a[" << order - 1 << "] * x2 + a[" << order - 3 << "];\n" ;
            ofs << "      t[1] = a[" << order - 2 << "] * x2 + a[" << order - 4 << "];\n" ;
            ofs << "      t[2] = b[" << order - 1 << "] * x2 + b[" << order - 3 << "];\n" ;
            ofs << "      t[3] = b[" << order - 2 << "] * x2 + b[" << order - 4 << "];\n" ;
            for(int item = order - 5; item >= 0; item -= 2)
            {
               ofs << "      t[0] *= x2;\n";
               if(item - 1 >= 0)
                  ofs << "      t[1] *= x2;\n";
               ofs << "      t[2] *= x2;\n";
               if(item - 1 >= 0)
                  ofs << "      t[3] *= x2;\n";
               ofs << "      t[0] += static_cast<V>(a[" << item << "]);\n";
               if(item - 1 >= 0)
                  ofs << "      t[1] += static_cast<V>(a[" << item - 1 << "]);\n";
               ofs << "      t[2] += static_cast<V>(b[" << item << "]);\n";
               if(item - 1 >= 0)
                  ofs << "      t[3] += static_cast<V>(b[" << item - 1 << "]);\n";
            }
            ofs << "      t[0] *= x;\n";
            ofs << "      t[2] *= x;\n";
         }
         ofs << "      return (t[0] + t[1]) / (t[2] + t[3]);\n";

         ofs << "   }\n   else\n   {\n      V z = 1 / x;\n      V z2 = 1 / (x * x);\n      V t[4];\n";

         if(order & 1)
         {
            ofs << "      t[0] = a[" << 0 << "] * z2 + a[" << 2 << "];\n" ;
            ofs << "      t[1] = a[" << 1 << "] * z2 + a[" << 3 << "];\n" ;
            ofs << "      t[2] = b[" << 0 << "] * z2 + b[" << 2 << "];\n" ;
            ofs << "      t[3] = b[" << 1 << "] * z2 + b[" << 3 << "];\n" ;
            for(int item = 4; item < order; item += 2)
            {
               ofs << "      t[0] *= z2;\n";
               if(item + 1 < order)
                  ofs << "      t[1] *= z2;\n";
               ofs << "      t[2] *= z2;\n";
               if(item + 1 < order)
                  ofs << "      t[3] *= z2;\n";
               ofs << "      t[0] += static_cast<V>(a[" << item << "]);\n";
               if(item + 1 < order)
                  ofs << "      t[1] += static_cast<V>(a[" << item + 1 << "]);\n";
               ofs << "      t[2] += static_cast<V>(b[" << item << "]);\n";
               if(item + 1 < order)
                  ofs << "      t[3] += static_cast<V>(b[" << item + 1 << "]);\n";
            }
            ofs << "      t[1] *= z;\n";
            ofs << "      t[3] *= z;\n";
         }
         else
         {
            ofs << "      t[0] = a[" << 0 << "] * z2 + a[" << 2 << "];\n" ;
            ofs << "      t[1] = a[" << 1 << "] * z2 + a[" << 3 << "];\n" ;
            ofs << "      t[2] = b[" << 0 << "] * z2 + b[" << 2 << "];\n" ;
            ofs << "      t[3] = b[" << 1 << "] * z2 + b[" << 3 << "];\n" ;
            for(int item = 4; item < order; item += 2)
            {
               ofs << "      t[0] *= z2;\n";
               if(item + 1 < order)
                  ofs << "      t[1] *= z2;\n";
               ofs << "      t[2] *= z2;\n";
               if(item + 1 < order)
                  ofs << "      t[3] *= z2;\n";
               ofs << "      t[0] += static_cast<V>(a[" << item << "]);\n";
               if(item + 1 < order)
                  ofs << "      t[1] += static_cast<V>(a[" << item + 1 << "]);\n";
               ofs << "      t[2] += static_cast<V>(b[" << item << "]);\n";
               if(item + 1 < order)
                  ofs << "      t[3] += static_cast<V>(b[" << item + 1 << "]);\n";
            }
            ofs << "      t[0] *= z;\n";
            ofs << "      t[2] *= z;\n";
         }
         ofs << "      return (t[0] + t[1]) / (t[2] + t[3]);\n   }\n";

         ofs << "}\n\n";
      }
      //
      // And finally the boilerplate at the end of the header:
      //
      ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";
   }
}

int main()
{
   for(int i = 2; i <= max_order; ++i)
   {
      print_polynomials(i);
      print_rationals(i);
   }
   return 0;
}