generate_rational_code.cpp

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

         "{\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;
}