mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2025-10-30 20:40:28 -04:00
739 lines
28 KiB
C++
739 lines
28 KiB
C++
// (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;
|
|
}
|
|
|
|
|
|
|