123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738 |
- // (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;
- }
|