summaryrefslogtreecommitdiffstats
path: root/src/boost/libs/math/tools/generate_rational_code.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/boost/libs/math/tools/generate_rational_code.cpp')
-rw-r--r--src/boost/libs/math/tools/generate_rational_code.cpp738
1 files changed, 738 insertions, 0 deletions
diff --git a/src/boost/libs/math/tools/generate_rational_code.cpp b/src/boost/libs/math/tools/generate_rational_code.cpp
new file mode 100644
index 00000000..4f4d62a0
--- /dev/null
+++ b/src/boost/libs/math/tools/generate_rational_code.cpp
@@ -0,0 +1,738 @@
+// (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*, 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*, 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*, 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*, 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"
+ "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*, 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"
+ "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;
+}
+
+
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-12-30 14:35:19 +0000