From 483eb2f56657e8e7f419ab1a4fab8dce9ade8609 Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Sat, 27 Apr 2024 20:24:20 +0200 Subject: Adding upstream version 14.2.21. Signed-off-by: Daniel Baumann --- .../libs/math/tools/generate_rational_code.cpp | 738 +++++++++++++++++++++ 1 file changed, 738 insertions(+) create mode 100644 src/boost/libs/math/tools/generate_rational_code.cpp (limited to 'src/boost/libs/math/tools/generate_rational_code.cpp') 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 +#include +#include +#include +#include + +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 \n" + "inline V evaluate_polynomial_c_imp(const T*, const V&, const mpl::int_<0>*)\n" + "{\n" + " return static_cast(0);\n" + "}\n" + "\n" + "template \n" + "inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<1>*)\n" + "{\n" + " return static_cast(a[0]);\n" + "}\n\n"; + + for(int order = 2; order <= i; ++order) + { + ofs << + "template \n" + "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<" << order << ">*)\n" + "{\n" + " return static_cast("; + + 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 \n" + "inline V evaluate_polynomial_c_imp(const T*, const V&, const mpl::int_<0>*)\n" + "{\n" + " return static_cast(0);\n" + "}\n" + "\n" + "template \n" + "inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<1>*)\n" + "{\n" + " return static_cast(a[0]);\n" + "}\n\n" + "template \n" + "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<2>*)\n" + "{\n" + " return static_cast(a[1] * x + a[0]);\n" + "}\n\n" + "template \n" + "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<3>*)\n" + "{\n" + " return static_cast((a[2] * x + a[1]) * x + a[0]);\n" + "}\n\n" + "template \n" + "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<4>*)\n" + "{\n" + " return static_cast(((a[3] * x + a[2]) * x + a[1]) * x + a[0]);\n" + "}\n\n"; + + for(int order = 5; order <= i; ++order) + { + ofs << + "template \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("; + + 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 \n" + "inline V evaluate_polynomial_c_imp(const T*, const V&, const mpl::int_<0>*)\n" + "{\n" + " return static_cast(0);\n" + "}\n" + "\n" + "template \n" + "inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<1>*)\n" + "{\n" + " return static_cast(a[0]);\n" + "}\n\n" + "template \n" + "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<2>*)\n" + "{\n" + " return static_cast(a[1] * x + a[0]);\n" + "}\n\n" + "template \n" + "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<3>*)\n" + "{\n" + " return static_cast((a[2] * x + a[1]) * x + a[0]);\n" + "}\n\n" + "template \n" + "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<4>*)\n" + "{\n" + " return static_cast(((a[3] * x + a[2]) * x + a[1]) * x + a[0]);\n" + "}\n\n"; + + for(int order = 5; order <= i; ++order) + { + ofs << + "template \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(a[" << order - 1 << "] * x2 + a[" << order - 3 << "]);\n" ; + ofs << " t[1] = static_cast(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(a[" << item << "]);\n"; + if(item - 1 >= 0) + ofs << " t[1] += static_cast(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(a[" << item << "]);\n"; + if(item - 1 >= 0) + ofs << " t[1] += static_cast(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 \n" + "inline V evaluate_rational_c_imp(const T*, const U*, const V&, const mpl::int_<0>*)\n" + "{\n" + " return static_cast(0);\n" + "}\n" + "\n" + "template \n" + "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const mpl::int_<1>*)\n" + "{\n" + " return static_cast(a[0]) / static_cast(b[0]);\n" + "}\n\n"; + + for(int order = 2; order <= i; ++order) + { + ofs << + "template \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(("; + + 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(("; + + 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 \n" + "inline V evaluate_rational_c_imp(const T*, const U*, const V&, const mpl::int_<0>*)\n" + "{\n" + " return static_cast(0);\n" + "}\n" + "\n" + "template \n" + "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const mpl::int_<1>*)\n" + "{\n" + " return static_cast(a[0]) / static_cast(b[0]);\n" + "}\n\n" + "template \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((a[1] * x + a[0]) / (b[1] * x + b[0]));\n" + "}\n\n" + "template \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(((a[2] * x + a[1]) * x + a[0]) / ((b[2] * x + b[1]) * x + b[0]));\n" + "}\n\n" + "template \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((((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 \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(("; + + 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(("; + + 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 \n" + "inline V evaluate_rational_c_imp(const T*, const U*, const V&, const mpl::int_<0>*)\n" + "{\n" + " return static_cast(0);\n" + "}\n" + "\n" + "template \n" + "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const mpl::int_<1>*)\n" + "{\n" + " return static_cast(a[0]) / static_cast(b[0]);\n" + "}\n\n" + "template \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((a[1] * x + a[0]) / (b[1] * x + b[0]));\n" + "}\n\n" + "template \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(((a[2] * x + a[1]) * x + a[0]) / ((b[2] * x + b[1]) * x + b[0]));\n" + "}\n\n" + "template \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((((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 \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(a[" << item << "]);\n"; + if(item - 1 >= 0) + ofs << " t[1] += static_cast(a[" << item - 1 << "]);\n"; + ofs << " t[2] += static_cast(b[" << item << "]);\n"; + if(item - 1 >= 0) + ofs << " t[3] += static_cast(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(a[" << item << "]);\n"; + if(item - 1 >= 0) + ofs << " t[1] += static_cast(a[" << item - 1 << "]);\n"; + ofs << " t[2] += static_cast(b[" << item << "]);\n"; + if(item - 1 >= 0) + ofs << " t[3] += static_cast(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(a[" << item << "]);\n"; + if(item + 1 < order) + ofs << " t[1] += static_cast(a[" << item + 1 << "]);\n"; + ofs << " t[2] += static_cast(b[" << item << "]);\n"; + if(item + 1 < order) + ofs << " t[3] += static_cast(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(a[" << item << "]);\n"; + if(item + 1 < order) + ofs << " t[1] += static_cast(a[" << item + 1 << "]);\n"; + ofs << " t[2] += static_cast(b[" << item << "]);\n"; + if(item + 1 < order) + ofs << " t[3] += static_cast(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; +} + + + -- cgit v1.2.3