tools/detail/rational_horner3_17.hpp

0001 //  (C) Copyright John Maddock 2007.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005 //
0006 //  This file is machine generated, do not edit by hand
0007
0008 // Polynomial evaluation using second order Horners rule
0009 #ifndef BOOST_MATH_TOOLS_RAT_EVAL_17_HPP
0010 #define BOOST_MATH_TOOLS_RAT_EVAL_17_HPP
0011
0012 namespace boost{ namespace math{ namespace tools{ namespace detail{
0013
0014 template <class T, class U, class V>
0015 inline V evaluate_rational_c_imp(const T*, const U*, const V&, const std::integral_constant<int, 0>*) BOOST_MATH_NOEXCEPT(V)
0016 {
0017    return static_cast<V>(0);
0018 }
0019
0020 template <class T, class U, class V>
0021 inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const std::integral_constant<int, 1>*) BOOST_MATH_NOEXCEPT(V)
0022 {
0023    return static_cast<V>(a[0]) / static_cast<V>(b[0]);
0024 }
0025
0026 template <class T, class U, class V>
0027 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 2>*) BOOST_MATH_NOEXCEPT(V)
0028 {
0029    return static_cast<V>((a[1] * x + a[0]) / (b[1] * x + b[0]));
0030 }
0031
0032 template <class T, class U, class V>
0033 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 3>*) BOOST_MATH_NOEXCEPT(V)
0034 {
0035    return static_cast<V>(((a[2] * x + a[1]) * x + a[0]) / ((b[2] * x + b[1]) * x + b[0]));
0036 }
0037
0038 template <class T, class U, class V>
0039 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 4>*) BOOST_MATH_NOEXCEPT(V)
0040 {
0041    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]));
0042 }
0043
0044 template <class T, class U, class V>
0045 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 5>*) BOOST_MATH_NOEXCEPT(V)
0046 {
0047    if(x <= 1)
0048    {
0049       V x2 = x * x;
0050       V t[4];
0051       t[0] = a[4] * x2 + a[2];
0052       t[1] = a[3] * x2 + a[1];
0053       t[2] = b[4] * x2 + b[2];
0054       t[3] = b[3] * x2 + b[1];
0055       t[0] *= x2;
0056       t[2] *= x2;
0057       t[0] += static_cast<V>(a[0]);
0058       t[2] += static_cast<V>(b[0]);
0059       t[1] *= x;
0060       t[3] *= x;
0061       return (t[0] + t[1]) / (t[2] + t[3]);
0062    }
0063    else
0064    {
0065       V z = 1 / x;
0066       V z2 = 1 / (x * x);
0067       V t[4];
0068       t[0] = a[0] * z2 + a[2];
0069       t[1] = a[1] * z2 + a[3];
0070       t[2] = b[0] * z2 + b[2];
0071       t[3] = b[1] * z2 + b[3];
0072       t[0] *= z2;
0073       t[2] *= z2;
0074       t[0] += static_cast<V>(a[4]);
0075       t[2] += static_cast<V>(b[4]);
0076       t[1] *= z;
0077       t[3] *= z;
0078       return (t[0] + t[1]) / (t[2] + t[3]);
0079    }
0080 }
0081
0082 template <class T, class U, class V>
0083 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 6>*) BOOST_MATH_NOEXCEPT(V)
0084 {
0085    if(x <= 1)
0086    {
0087       V x2 = x * x;
0088       V t[4];
0089       t[0] = a[5] * x2 + a[3];
0090       t[1] = a[4] * x2 + a[2];
0091       t[2] = b[5] * x2 + b[3];
0092       t[3] = b[4] * x2 + b[2];
0093       t[0] *= x2;
0094       t[1] *= x2;
0095       t[2] *= x2;
0096       t[3] *= x2;
0097       t[0] += static_cast<V>(a[1]);
0098       t[1] += static_cast<V>(a[0]);
0099       t[2] += static_cast<V>(b[1]);
0100       t[3] += static_cast<V>(b[0]);
0101       t[0] *= x;
0102       t[2] *= x;
0103       return (t[0] + t[1]) / (t[2] + t[3]);
0104    }
0105    else
0106    {
0107       V z = 1 / x;
0108       V z2 = 1 / (x * x);
0109       V t[4];
0110       t[0] = a[0] * z2 + a[2];
0111       t[1] = a[1] * z2 + a[3];
0112       t[2] = b[0] * z2 + b[2];
0113       t[3] = b[1] * z2 + b[3];
0114       t[0] *= z2;
0115       t[1] *= z2;
0116       t[2] *= z2;
0117       t[3] *= z2;
0118       t[0] += static_cast<V>(a[4]);
0119       t[1] += static_cast<V>(a[5]);
0120       t[2] += static_cast<V>(b[4]);
0121       t[3] += static_cast<V>(b[5]);
0122       t[0] *= z;
0123       t[2] *= z;
0124       return (t[0] + t[1]) / (t[2] + t[3]);
0125    }
0126 }
0127
0128 template <class T, class U, class V>
0129 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 7>*) BOOST_MATH_NOEXCEPT(V)
0130 {
0131    if(x <= 1)
0132    {
0133       V x2 = x * x;
0134       V t[4];
0135       t[0] = a[6] * x2 + a[4];
0136       t[1] = a[5] * x2 + a[3];
0137       t[2] = b[6] * x2 + b[4];
0138       t[3] = b[5] * x2 + b[3];
0139       t[0] *= x2;
0140       t[1] *= x2;
0141       t[2] *= x2;
0142       t[3] *= x2;
0143       t[0] += static_cast<V>(a[2]);
0144       t[1] += static_cast<V>(a[1]);
0145       t[2] += static_cast<V>(b[2]);
0146       t[3] += static_cast<V>(b[1]);
0147       t[0] *= x2;
0148       t[2] *= x2;
0149       t[0] += static_cast<V>(a[0]);
0150       t[2] += static_cast<V>(b[0]);
0151       t[1] *= x;
0152       t[3] *= x;
0153       return (t[0] + t[1]) / (t[2] + t[3]);
0154    }
0155    else
0156    {
0157       V z = 1 / x;
0158       V z2 = 1 / (x * x);
0159       V t[4];
0160       t[0] = a[0] * z2 + a[2];
0161       t[1] = a[1] * z2 + a[3];
0162       t[2] = b[0] * z2 + b[2];
0163       t[3] = b[1] * z2 + b[3];
0164       t[0] *= z2;
0165       t[1] *= z2;
0166       t[2] *= z2;
0167       t[3] *= z2;
0168       t[0] += static_cast<V>(a[4]);
0169       t[1] += static_cast<V>(a[5]);
0170       t[2] += static_cast<V>(b[4]);
0171       t[3] += static_cast<V>(b[5]);
0172       t[0] *= z2;
0173       t[2] *= z2;
0174       t[0] += static_cast<V>(a[6]);
0175       t[2] += static_cast<V>(b[6]);
0176       t[1] *= z;
0177       t[3] *= z;
0178       return (t[0] + t[1]) / (t[2] + t[3]);
0179    }
0180 }
0181
0182 template <class T, class U, class V>
0183 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 8>*) BOOST_MATH_NOEXCEPT(V)
0184 {
0185    if(x <= 1)
0186    {
0187       V x2 = x * x;
0188       V t[4];
0189       t[0] = a[7] * x2 + a[5];
0190       t[1] = a[6] * x2 + a[4];
0191       t[2] = b[7] * x2 + b[5];
0192       t[3] = b[6] * x2 + b[4];
0193       t[0] *= x2;
0194       t[1] *= x2;
0195       t[2] *= x2;
0196       t[3] *= x2;
0197       t[0] += static_cast<V>(a[3]);
0198       t[1] += static_cast<V>(a[2]);
0199       t[2] += static_cast<V>(b[3]);
0200       t[3] += static_cast<V>(b[2]);
0201       t[0] *= x2;
0202       t[1] *= x2;
0203       t[2] *= x2;
0204       t[3] *= x2;
0205       t[0] += static_cast<V>(a[1]);
0206       t[1] += static_cast<V>(a[0]);
0207       t[2] += static_cast<V>(b[1]);
0208       t[3] += static_cast<V>(b[0]);
0209       t[0] *= x;
0210       t[2] *= x;
0211       return (t[0] + t[1]) / (t[2] + t[3]);
0212    }
0213    else
0214    {
0215       V z = 1 / x;
0216       V z2 = 1 / (x * x);
0217       V t[4];
0218       t[0] = a[0] * z2 + a[2];
0219       t[1] = a[1] * z2 + a[3];
0220       t[2] = b[0] * z2 + b[2];
0221       t[3] = b[1] * z2 + b[3];
0222       t[0] *= z2;
0223       t[1] *= z2;
0224       t[2] *= z2;
0225       t[3] *= z2;
0226       t[0] += static_cast<V>(a[4]);
0227       t[1] += static_cast<V>(a[5]);
0228       t[2] += static_cast<V>(b[4]);
0229       t[3] += static_cast<V>(b[5]);
0230       t[0] *= z2;
0231       t[1] *= z2;
0232       t[2] *= z2;
0233       t[3] *= z2;
0234       t[0] += static_cast<V>(a[6]);
0235       t[1] += static_cast<V>(a[7]);
0236       t[2] += static_cast<V>(b[6]);
0237       t[3] += static_cast<V>(b[7]);
0238       t[0] *= z;
0239       t[2] *= z;
0240       return (t[0] + t[1]) / (t[2] + t[3]);
0241    }
0242 }
0243
0244 template <class T, class U, class V>
0245 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 9>*) BOOST_MATH_NOEXCEPT(V)
0246 {
0247    if(x <= 1)
0248    {
0249       V x2 = x * x;
0250       V t[4];
0251       t[0] = a[8] * x2 + a[6];
0252       t[1] = a[7] * x2 + a[5];
0253       t[2] = b[8] * x2 + b[6];
0254       t[3] = b[7] * x2 + b[5];
0255       t[0] *= x2;
0256       t[1] *= x2;
0257       t[2] *= x2;
0258       t[3] *= x2;
0259       t[0] += static_cast<V>(a[4]);
0260       t[1] += static_cast<V>(a[3]);
0261       t[2] += static_cast<V>(b[4]);
0262       t[3] += static_cast<V>(b[3]);
0263       t[0] *= x2;
0264       t[1] *= x2;
0265       t[2] *= x2;
0266       t[3] *= x2;
0267       t[0] += static_cast<V>(a[2]);
0268       t[1] += static_cast<V>(a[1]);
0269       t[2] += static_cast<V>(b[2]);
0270       t[3] += static_cast<V>(b[1]);
0271       t[0] *= x2;
0272       t[2] *= x2;
0273       t[0] += static_cast<V>(a[0]);
0274       t[2] += static_cast<V>(b[0]);
0275       t[1] *= x;
0276       t[3] *= x;
0277       return (t[0] + t[1]) / (t[2] + t[3]);
0278    }
0279    else
0280    {
0281       V z = 1 / x;
0282       V z2 = 1 / (x * x);
0283       V t[4];
0284       t[0] = a[0] * z2 + a[2];
0285       t[1] = a[1] * z2 + a[3];
0286       t[2] = b[0] * z2 + b[2];
0287       t[3] = b[1] * z2 + b[3];
0288       t[0] *= z2;
0289       t[1] *= z2;
0290       t[2] *= z2;
0291       t[3] *= z2;
0292       t[0] += static_cast<V>(a[4]);
0293       t[1] += static_cast<V>(a[5]);
0294       t[2] += static_cast<V>(b[4]);
0295       t[3] += static_cast<V>(b[5]);
0296       t[0] *= z2;
0297       t[1] *= z2;
0298       t[2] *= z2;
0299       t[3] *= z2;
0300       t[0] += static_cast<V>(a[6]);
0301       t[1] += static_cast<V>(a[7]);
0302       t[2] += static_cast<V>(b[6]);
0303       t[3] += static_cast<V>(b[7]);
0304       t[0] *= z2;
0305       t[2] *= z2;
0306       t[0] += static_cast<V>(a[8]);
0307       t[2] += static_cast<V>(b[8]);
0308       t[1] *= z;
0309       t[3] *= z;
0310       return (t[0] + t[1]) / (t[2] + t[3]);
0311    }
0312 }
0313
0314 template <class T, class U, class V>
0315 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 10>*) BOOST_MATH_NOEXCEPT(V)
0316 {
0317    if(x <= 1)
0318    {
0319       V x2 = x * x;
0320       V t[4];
0321       t[0] = a[9] * x2 + a[7];
0322       t[1] = a[8] * x2 + a[6];
0323       t[2] = b[9] * x2 + b[7];
0324       t[3] = b[8] * x2 + b[6];
0325       t[0] *= x2;
0326       t[1] *= x2;
0327       t[2] *= x2;
0328       t[3] *= x2;
0329       t[0] += static_cast<V>(a[5]);
0330       t[1] += static_cast<V>(a[4]);
0331       t[2] += static_cast<V>(b[5]);
0332       t[3] += static_cast<V>(b[4]);
0333       t[0] *= x2;
0334       t[1] *= x2;
0335       t[2] *= x2;
0336       t[3] *= x2;
0337       t[0] += static_cast<V>(a[3]);
0338       t[1] += static_cast<V>(a[2]);
0339       t[2] += static_cast<V>(b[3]);
0340       t[3] += static_cast<V>(b[2]);
0341       t[0] *= x2;
0342       t[1] *= x2;
0343       t[2] *= x2;
0344       t[3] *= x2;
0345       t[0] += static_cast<V>(a[1]);
0346       t[1] += static_cast<V>(a[0]);
0347       t[2] += static_cast<V>(b[1]);
0348       t[3] += static_cast<V>(b[0]);
0349       t[0] *= x;
0350       t[2] *= x;
0351       return (t[0] + t[1]) / (t[2] + t[3]);
0352    }
0353    else
0354    {
0355       V z = 1 / x;
0356       V z2 = 1 / (x * x);
0357       V t[4];
0358       t[0] = a[0] * z2 + a[2];
0359       t[1] = a[1] * z2 + a[3];
0360       t[2] = b[0] * z2 + b[2];
0361       t[3] = b[1] * z2 + b[3];
0362       t[0] *= z2;
0363       t[1] *= z2;
0364       t[2] *= z2;
0365       t[3] *= z2;
0366       t[0] += static_cast<V>(a[4]);
0367       t[1] += static_cast<V>(a[5]);
0368       t[2] += static_cast<V>(b[4]);
0369       t[3] += static_cast<V>(b[5]);
0370       t[0] *= z2;
0371       t[1] *= z2;
0372       t[2] *= z2;
0373       t[3] *= z2;
0374       t[0] += static_cast<V>(a[6]);
0375       t[1] += static_cast<V>(a[7]);
0376       t[2] += static_cast<V>(b[6]);
0377       t[3] += static_cast<V>(b[7]);
0378       t[0] *= z2;
0379       t[1] *= z2;
0380       t[2] *= z2;
0381       t[3] *= z2;
0382       t[0] += static_cast<V>(a[8]);
0383       t[1] += static_cast<V>(a[9]);
0384       t[2] += static_cast<V>(b[8]);
0385       t[3] += static_cast<V>(b[9]);
0386       t[0] *= z;
0387       t[2] *= z;
0388       return (t[0] + t[1]) / (t[2] + t[3]);
0389    }
0390 }
0391
0392 template <class T, class U, class V>
0393 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 11>*) BOOST_MATH_NOEXCEPT(V)
0394 {
0395    if(x <= 1)
0396    {
0397       V x2 = x * x;
0398       V t[4];
0399       t[0] = a[10] * x2 + a[8];
0400       t[1] = a[9] * x2 + a[7];
0401       t[2] = b[10] * x2 + b[8];
0402       t[3] = b[9] * x2 + b[7];
0403       t[0] *= x2;
0404       t[1] *= x2;
0405       t[2] *= x2;
0406       t[3] *= x2;
0407       t[0] += static_cast<V>(a[6]);
0408       t[1] += static_cast<V>(a[5]);
0409       t[2] += static_cast<V>(b[6]);
0410       t[3] += static_cast<V>(b[5]);
0411       t[0] *= x2;
0412       t[1] *= x2;
0413       t[2] *= x2;
0414       t[3] *= x2;
0415       t[0] += static_cast<V>(a[4]);
0416       t[1] += static_cast<V>(a[3]);
0417       t[2] += static_cast<V>(b[4]);
0418       t[3] += static_cast<V>(b[3]);
0419       t[0] *= x2;
0420       t[1] *= x2;
0421       t[2] *= x2;
0422       t[3] *= x2;
0423       t[0] += static_cast<V>(a[2]);
0424       t[1] += static_cast<V>(a[1]);
0425       t[2] += static_cast<V>(b[2]);
0426       t[3] += static_cast<V>(b[1]);
0427       t[0] *= x2;
0428       t[2] *= x2;
0429       t[0] += static_cast<V>(a[0]);
0430       t[2] += static_cast<V>(b[0]);
0431       t[1] *= x;
0432       t[3] *= x;
0433       return (t[0] + t[1]) / (t[2] + t[3]);
0434    }
0435    else
0436    {
0437       V z = 1 / x;
0438       V z2 = 1 / (x * x);
0439       V t[4];
0440       t[0] = a[0] * z2 + a[2];
0441       t[1] = a[1] * z2 + a[3];
0442       t[2] = b[0] * z2 + b[2];
0443       t[3] = b[1] * z2 + b[3];
0444       t[0] *= z2;
0445       t[1] *= z2;
0446       t[2] *= z2;
0447       t[3] *= z2;
0448       t[0] += static_cast<V>(a[4]);
0449       t[1] += static_cast<V>(a[5]);
0450       t[2] += static_cast<V>(b[4]);
0451       t[3] += static_cast<V>(b[5]);
0452       t[0] *= z2;
0453       t[1] *= z2;
0454       t[2] *= z2;
0455       t[3] *= z2;
0456       t[0] += static_cast<V>(a[6]);
0457       t[1] += static_cast<V>(a[7]);
0458       t[2] += static_cast<V>(b[6]);
0459       t[3] += static_cast<V>(b[7]);
0460       t[0] *= z2;
0461       t[1] *= z2;
0462       t[2] *= z2;
0463       t[3] *= z2;
0464       t[0] += static_cast<V>(a[8]);
0465       t[1] += static_cast<V>(a[9]);
0466       t[2] += static_cast<V>(b[8]);
0467       t[3] += static_cast<V>(b[9]);
0468       t[0] *= z2;
0469       t[2] *= z2;
0470       t[0] += static_cast<V>(a[10]);
0471       t[2] += static_cast<V>(b[10]);
0472       t[1] *= z;
0473       t[3] *= z;
0474       return (t[0] + t[1]) / (t[2] + t[3]);
0475    }
0476 }
0477
0478 template <class T, class U, class V>
0479 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 12>*) BOOST_MATH_NOEXCEPT(V)
0480 {
0481    if(x <= 1)
0482    {
0483       V x2 = x * x;
0484       V t[4];
0485       t[0] = a[11] * x2 + a[9];
0486       t[1] = a[10] * x2 + a[8];
0487       t[2] = b[11] * x2 + b[9];
0488       t[3] = b[10] * x2 + b[8];
0489       t[0] *= x2;
0490       t[1] *= x2;
0491       t[2] *= x2;
0492       t[3] *= x2;
0493       t[0] += static_cast<V>(a[7]);
0494       t[1] += static_cast<V>(a[6]);
0495       t[2] += static_cast<V>(b[7]);
0496       t[3] += static_cast<V>(b[6]);
0497       t[0] *= x2;
0498       t[1] *= x2;
0499       t[2] *= x2;
0500       t[3] *= x2;
0501       t[0] += static_cast<V>(a[5]);
0502       t[1] += static_cast<V>(a[4]);
0503       t[2] += static_cast<V>(b[5]);
0504       t[3] += static_cast<V>(b[4]);
0505       t[0] *= x2;
0506       t[1] *= x2;
0507       t[2] *= x2;
0508       t[3] *= x2;
0509       t[0] += static_cast<V>(a[3]);
0510       t[1] += static_cast<V>(a[2]);
0511       t[2] += static_cast<V>(b[3]);
0512       t[3] += static_cast<V>(b[2]);
0513       t[0] *= x2;
0514       t[1] *= x2;
0515       t[2] *= x2;
0516       t[3] *= x2;
0517       t[0] += static_cast<V>(a[1]);
0518       t[1] += static_cast<V>(a[0]);
0519       t[2] += static_cast<V>(b[1]);
0520       t[3] += static_cast<V>(b[0]);
0521       t[0] *= x;
0522       t[2] *= x;
0523       return (t[0] + t[1]) / (t[2] + t[3]);
0524    }
0525    else
0526    {
0527       V z = 1 / x;
0528       V z2 = 1 / (x * x);
0529       V t[4];
0530       t[0] = a[0] * z2 + a[2];
0531       t[1] = a[1] * z2 + a[3];
0532       t[2] = b[0] * z2 + b[2];
0533       t[3] = b[1] * z2 + b[3];
0534       t[0] *= z2;
0535       t[1] *= z2;
0536       t[2] *= z2;
0537       t[3] *= z2;
0538       t[0] += static_cast<V>(a[4]);
0539       t[1] += static_cast<V>(a[5]);
0540       t[2] += static_cast<V>(b[4]);
0541       t[3] += static_cast<V>(b[5]);
0542       t[0] *= z2;
0543       t[1] *= z2;
0544       t[2] *= z2;
0545       t[3] *= z2;
0546       t[0] += static_cast<V>(a[6]);
0547       t[1] += static_cast<V>(a[7]);
0548       t[2] += static_cast<V>(b[6]);
0549       t[3] += static_cast<V>(b[7]);
0550       t[0] *= z2;
0551       t[1] *= z2;
0552       t[2] *= z2;
0553       t[3] *= z2;
0554       t[0] += static_cast<V>(a[8]);
0555       t[1] += static_cast<V>(a[9]);
0556       t[2] += static_cast<V>(b[8]);
0557       t[3] += static_cast<V>(b[9]);
0558       t[0] *= z2;
0559       t[1] *= z2;
0560       t[2] *= z2;
0561       t[3] *= z2;
0562       t[0] += static_cast<V>(a[10]);
0563       t[1] += static_cast<V>(a[11]);
0564       t[2] += static_cast<V>(b[10]);
0565       t[3] += static_cast<V>(b[11]);
0566       t[0] *= z;
0567       t[2] *= z;
0568       return (t[0] + t[1]) / (t[2] + t[3]);
0569    }
0570 }
0571
0572 template <class T, class U, class V>
0573 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 13>*) BOOST_MATH_NOEXCEPT(V)
0574 {
0575    if(x <= 1)
0576    {
0577       V x2 = x * x;
0578       V t[4];
0579       t[0] = a[12] * x2 + a[10];
0580       t[1] = a[11] * x2 + a[9];
0581       t[2] = b[12] * x2 + b[10];
0582       t[3] = b[11] * x2 + b[9];
0583       t[0] *= x2;
0584       t[1] *= x2;
0585       t[2] *= x2;
0586       t[3] *= x2;
0587       t[0] += static_cast<V>(a[8]);
0588       t[1] += static_cast<V>(a[7]);
0589       t[2] += static_cast<V>(b[8]);
0590       t[3] += static_cast<V>(b[7]);
0591       t[0] *= x2;
0592       t[1] *= x2;
0593       t[2] *= x2;
0594       t[3] *= x2;
0595       t[0] += static_cast<V>(a[6]);
0596       t[1] += static_cast<V>(a[5]);
0597       t[2] += static_cast<V>(b[6]);
0598       t[3] += static_cast<V>(b[5]);
0599       t[0] *= x2;
0600       t[1] *= x2;
0601       t[2] *= x2;
0602       t[3] *= x2;
0603       t[0] += static_cast<V>(a[4]);
0604       t[1] += static_cast<V>(a[3]);
0605       t[2] += static_cast<V>(b[4]);
0606       t[3] += static_cast<V>(b[3]);
0607       t[0] *= x2;
0608       t[1] *= x2;
0609       t[2] *= x2;
0610       t[3] *= x2;
0611       t[0] += static_cast<V>(a[2]);
0612       t[1] += static_cast<V>(a[1]);
0613       t[2] += static_cast<V>(b[2]);
0614       t[3] += static_cast<V>(b[1]);
0615       t[0] *= x2;
0616       t[2] *= x2;
0617       t[0] += static_cast<V>(a[0]);
0618       t[2] += static_cast<V>(b[0]);
0619       t[1] *= x;
0620       t[3] *= x;
0621       return (t[0] + t[1]) / (t[2] + t[3]);
0622    }
0623    else
0624    {
0625       V z = 1 / x;
0626       V z2 = 1 / (x * x);
0627       V t[4];
0628       t[0] = a[0] * z2 + a[2];
0629       t[1] = a[1] * z2 + a[3];
0630       t[2] = b[0] * z2 + b[2];
0631       t[3] = b[1] * z2 + b[3];
0632       t[0] *= z2;
0633       t[1] *= z2;
0634       t[2] *= z2;
0635       t[3] *= z2;
0636       t[0] += static_cast<V>(a[4]);
0637       t[1] += static_cast<V>(a[5]);
0638       t[2] += static_cast<V>(b[4]);
0639       t[3] += static_cast<V>(b[5]);
0640       t[0] *= z2;
0641       t[1] *= z2;
0642       t[2] *= z2;
0643       t[3] *= z2;
0644       t[0] += static_cast<V>(a[6]);
0645       t[1] += static_cast<V>(a[7]);
0646       t[2] += static_cast<V>(b[6]);
0647       t[3] += static_cast<V>(b[7]);
0648       t[0] *= z2;
0649       t[1] *= z2;
0650       t[2] *= z2;
0651       t[3] *= z2;
0652       t[0] += static_cast<V>(a[8]);
0653       t[1] += static_cast<V>(a[9]);
0654       t[2] += static_cast<V>(b[8]);
0655       t[3] += static_cast<V>(b[9]);
0656       t[0] *= z2;
0657       t[1] *= z2;
0658       t[2] *= z2;
0659       t[3] *= z2;
0660       t[0] += static_cast<V>(a[10]);
0661       t[1] += static_cast<V>(a[11]);
0662       t[2] += static_cast<V>(b[10]);
0663       t[3] += static_cast<V>(b[11]);
0664       t[0] *= z2;
0665       t[2] *= z2;
0666       t[0] += static_cast<V>(a[12]);
0667       t[2] += static_cast<V>(b[12]);
0668       t[1] *= z;
0669       t[3] *= z;
0670       return (t[0] + t[1]) / (t[2] + t[3]);
0671    }
0672 }
0673
0674 template <class T, class U, class V>
0675 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 14>*) BOOST_MATH_NOEXCEPT(V)
0676 {
0677    if(x <= 1)
0678    {
0679       V x2 = x * x;
0680       V t[4];
0681       t[0] = a[13] * x2 + a[11];
0682       t[1] = a[12] * x2 + a[10];
0683       t[2] = b[13] * x2 + b[11];
0684       t[3] = b[12] * x2 + b[10];
0685       t[0] *= x2;
0686       t[1] *= x2;
0687       t[2] *= x2;
0688       t[3] *= x2;
0689       t[0] += static_cast<V>(a[9]);
0690       t[1] += static_cast<V>(a[8]);
0691       t[2] += static_cast<V>(b[9]);
0692       t[3] += static_cast<V>(b[8]);
0693       t[0] *= x2;
0694       t[1] *= x2;
0695       t[2] *= x2;
0696       t[3] *= x2;
0697       t[0] += static_cast<V>(a[7]);
0698       t[1] += static_cast<V>(a[6]);
0699       t[2] += static_cast<V>(b[7]);
0700       t[3] += static_cast<V>(b[6]);
0701       t[0] *= x2;
0702       t[1] *= x2;
0703       t[2] *= x2;
0704       t[3] *= x2;
0705       t[0] += static_cast<V>(a[5]);
0706       t[1] += static_cast<V>(a[4]);
0707       t[2] += static_cast<V>(b[5]);
0708       t[3] += static_cast<V>(b[4]);
0709       t[0] *= x2;
0710       t[1] *= x2;
0711       t[2] *= x2;
0712       t[3] *= x2;
0713       t[0] += static_cast<V>(a[3]);
0714       t[1] += static_cast<V>(a[2]);
0715       t[2] += static_cast<V>(b[3]);
0716       t[3] += static_cast<V>(b[2]);
0717       t[0] *= x2;
0718       t[1] *= x2;
0719       t[2] *= x2;
0720       t[3] *= x2;
0721       t[0] += static_cast<V>(a[1]);
0722       t[1] += static_cast<V>(a[0]);
0723       t[2] += static_cast<V>(b[1]);
0724       t[3] += static_cast<V>(b[0]);
0725       t[0] *= x;
0726       t[2] *= x;
0727       return (t[0] + t[1]) / (t[2] + t[3]);
0728    }
0729    else
0730    {
0731       V z = 1 / x;
0732       V z2 = 1 / (x * x);
0733       V t[4];
0734       t[0] = a[0] * z2 + a[2];
0735       t[1] = a[1] * z2 + a[3];
0736       t[2] = b[0] * z2 + b[2];
0737       t[3] = b[1] * z2 + b[3];
0738       t[0] *= z2;
0739       t[1] *= z2;
0740       t[2] *= z2;
0741       t[3] *= z2;
0742       t[0] += static_cast<V>(a[4]);
0743       t[1] += static_cast<V>(a[5]);
0744       t[2] += static_cast<V>(b[4]);
0745       t[3] += static_cast<V>(b[5]);
0746       t[0] *= z2;
0747       t[1] *= z2;
0748       t[2] *= z2;
0749       t[3] *= z2;
0750       t[0] += static_cast<V>(a[6]);
0751       t[1] += static_cast<V>(a[7]);
0752       t[2] += static_cast<V>(b[6]);
0753       t[3] += static_cast<V>(b[7]);
0754       t[0] *= z2;
0755       t[1] *= z2;
0756       t[2] *= z2;
0757       t[3] *= z2;
0758       t[0] += static_cast<V>(a[8]);
0759       t[1] += static_cast<V>(a[9]);
0760       t[2] += static_cast<V>(b[8]);
0761       t[3] += static_cast<V>(b[9]);
0762       t[0] *= z2;
0763       t[1] *= z2;
0764       t[2] *= z2;
0765       t[3] *= z2;
0766       t[0] += static_cast<V>(a[10]);
0767       t[1] += static_cast<V>(a[11]);
0768       t[2] += static_cast<V>(b[10]);
0769       t[3] += static_cast<V>(b[11]);
0770       t[0] *= z2;
0771       t[1] *= z2;
0772       t[2] *= z2;
0773       t[3] *= z2;
0774       t[0] += static_cast<V>(a[12]);
0775       t[1] += static_cast<V>(a[13]);
0776       t[2] += static_cast<V>(b[12]);
0777       t[3] += static_cast<V>(b[13]);
0778       t[0] *= z;
0779       t[2] *= z;
0780       return (t[0] + t[1]) / (t[2] + t[3]);
0781    }
0782 }
0783
0784 template <class T, class U, class V>
0785 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 15>*) BOOST_MATH_NOEXCEPT(V)
0786 {
0787    if(x <= 1)
0788    {
0789       V x2 = x * x;
0790       V t[4];
0791       t[0] = a[14] * x2 + a[12];
0792       t[1] = a[13] * x2 + a[11];
0793       t[2] = b[14] * x2 + b[12];
0794       t[3] = b[13] * x2 + b[11];
0795       t[0] *= x2;
0796       t[1] *= x2;
0797       t[2] *= x2;
0798       t[3] *= x2;
0799       t[0] += static_cast<V>(a[10]);
0800       t[1] += static_cast<V>(a[9]);
0801       t[2] += static_cast<V>(b[10]);
0802       t[3] += static_cast<V>(b[9]);
0803       t[0] *= x2;
0804       t[1] *= x2;
0805       t[2] *= x2;
0806       t[3] *= x2;
0807       t[0] += static_cast<V>(a[8]);
0808       t[1] += static_cast<V>(a[7]);
0809       t[2] += static_cast<V>(b[8]);
0810       t[3] += static_cast<V>(b[7]);
0811       t[0] *= x2;
0812       t[1] *= x2;
0813       t[2] *= x2;
0814       t[3] *= x2;
0815       t[0] += static_cast<V>(a[6]);
0816       t[1] += static_cast<V>(a[5]);
0817       t[2] += static_cast<V>(b[6]);
0818       t[3] += static_cast<V>(b[5]);
0819       t[0] *= x2;
0820       t[1] *= x2;
0821       t[2] *= x2;
0822       t[3] *= x2;
0823       t[0] += static_cast<V>(a[4]);
0824       t[1] += static_cast<V>(a[3]);
0825       t[2] += static_cast<V>(b[4]);
0826       t[3] += static_cast<V>(b[3]);
0827       t[0] *= x2;
0828       t[1] *= x2;
0829       t[2] *= x2;
0830       t[3] *= x2;
0831       t[0] += static_cast<V>(a[2]);
0832       t[1] += static_cast<V>(a[1]);
0833       t[2] += static_cast<V>(b[2]);
0834       t[3] += static_cast<V>(b[1]);
0835       t[0] *= x2;
0836       t[2] *= x2;
0837       t[0] += static_cast<V>(a[0]);
0838       t[2] += static_cast<V>(b[0]);
0839       t[1] *= x;
0840       t[3] *= x;
0841       return (t[0] + t[1]) / (t[2] + t[3]);
0842    }
0843    else
0844    {
0845       V z = 1 / x;
0846       V z2 = 1 / (x * x);
0847       V t[4];
0848       t[0] = a[0] * z2 + a[2];
0849       t[1] = a[1] * z2 + a[3];
0850       t[2] = b[0] * z2 + b[2];
0851       t[3] = b[1] * z2 + b[3];
0852       t[0] *= z2;
0853       t[1] *= z2;
0854       t[2] *= z2;
0855       t[3] *= z2;
0856       t[0] += static_cast<V>(a[4]);
0857       t[1] += static_cast<V>(a[5]);
0858       t[2] += static_cast<V>(b[4]);
0859       t[3] += static_cast<V>(b[5]);
0860       t[0] *= z2;
0861       t[1] *= z2;
0862       t[2] *= z2;
0863       t[3] *= z2;
0864       t[0] += static_cast<V>(a[6]);
0865       t[1] += static_cast<V>(a[7]);
0866       t[2] += static_cast<V>(b[6]);
0867       t[3] += static_cast<V>(b[7]);
0868       t[0] *= z2;
0869       t[1] *= z2;
0870       t[2] *= z2;
0871       t[3] *= z2;
0872       t[0] += static_cast<V>(a[8]);
0873       t[1] += static_cast<V>(a[9]);
0874       t[2] += static_cast<V>(b[8]);
0875       t[3] += static_cast<V>(b[9]);
0876       t[0] *= z2;
0877       t[1] *= z2;
0878       t[2] *= z2;
0879       t[3] *= z2;
0880       t[0] += static_cast<V>(a[10]);
0881       t[1] += static_cast<V>(a[11]);
0882       t[2] += static_cast<V>(b[10]);
0883       t[3] += static_cast<V>(b[11]);
0884       t[0] *= z2;
0885       t[1] *= z2;
0886       t[2] *= z2;
0887       t[3] *= z2;
0888       t[0] += static_cast<V>(a[12]);
0889       t[1] += static_cast<V>(a[13]);
0890       t[2] += static_cast<V>(b[12]);
0891       t[3] += static_cast<V>(b[13]);
0892       t[0] *= z2;
0893       t[2] *= z2;
0894       t[0] += static_cast<V>(a[14]);
0895       t[2] += static_cast<V>(b[14]);
0896       t[1] *= z;
0897       t[3] *= z;
0898       return (t[0] + t[1]) / (t[2] + t[3]);
0899    }
0900 }
0901
0902 template <class T, class U, class V>
0903 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 16>*) BOOST_MATH_NOEXCEPT(V)
0904 {
0905    if(x <= 1)
0906    {
0907       V x2 = x * x;
0908       V t[4];
0909       t[0] = a[15] * x2 + a[13];
0910       t[1] = a[14] * x2 + a[12];
0911       t[2] = b[15] * x2 + b[13];
0912       t[3] = b[14] * x2 + b[12];
0913       t[0] *= x2;
0914       t[1] *= x2;
0915       t[2] *= x2;
0916       t[3] *= x2;
0917       t[0] += static_cast<V>(a[11]);
0918       t[1] += static_cast<V>(a[10]);
0919       t[2] += static_cast<V>(b[11]);
0920       t[3] += static_cast<V>(b[10]);
0921       t[0] *= x2;
0922       t[1] *= x2;
0923       t[2] *= x2;
0924       t[3] *= x2;
0925       t[0] += static_cast<V>(a[9]);
0926       t[1] += static_cast<V>(a[8]);
0927       t[2] += static_cast<V>(b[9]);
0928       t[3] += static_cast<V>(b[8]);
0929       t[0] *= x2;
0930       t[1] *= x2;
0931       t[2] *= x2;
0932       t[3] *= x2;
0933       t[0] += static_cast<V>(a[7]);
0934       t[1] += static_cast<V>(a[6]);
0935       t[2] += static_cast<V>(b[7]);
0936       t[3] += static_cast<V>(b[6]);
0937       t[0] *= x2;
0938       t[1] *= x2;
0939       t[2] *= x2;
0940       t[3] *= x2;
0941       t[0] += static_cast<V>(a[5]);
0942       t[1] += static_cast<V>(a[4]);
0943       t[2] += static_cast<V>(b[5]);
0944       t[3] += static_cast<V>(b[4]);
0945       t[0] *= x2;
0946       t[1] *= x2;
0947       t[2] *= x2;
0948       t[3] *= x2;
0949       t[0] += static_cast<V>(a[3]);
0950       t[1] += static_cast<V>(a[2]);
0951       t[2] += static_cast<V>(b[3]);
0952       t[3] += static_cast<V>(b[2]);
0953       t[0] *= x2;
0954       t[1] *= x2;
0955       t[2] *= x2;
0956       t[3] *= x2;
0957       t[0] += static_cast<V>(a[1]);
0958       t[1] += static_cast<V>(a[0]);
0959       t[2] += static_cast<V>(b[1]);
0960       t[3] += static_cast<V>(b[0]);
0961       t[0] *= x;
0962       t[2] *= x;
0963       return (t[0] + t[1]) / (t[2] + t[3]);
0964    }
0965    else
0966    {
0967       V z = 1 / x;
0968       V z2 = 1 / (x * x);
0969       V t[4];
0970       t[0] = a[0] * z2 + a[2];
0971       t[1] = a[1] * z2 + a[3];
0972       t[2] = b[0] * z2 + b[2];
0973       t[3] = b[1] * z2 + b[3];
0974       t[0] *= z2;
0975       t[1] *= z2;
0976       t[2] *= z2;
0977       t[3] *= z2;
0978       t[0] += static_cast<V>(a[4]);
0979       t[1] += static_cast<V>(a[5]);
0980       t[2] += static_cast<V>(b[4]);
0981       t[3] += static_cast<V>(b[5]);
0982       t[0] *= z2;
0983       t[1] *= z2;
0984       t[2] *= z2;
0985       t[3] *= z2;
0986       t[0] += static_cast<V>(a[6]);
0987       t[1] += static_cast<V>(a[7]);
0988       t[2] += static_cast<V>(b[6]);
0989       t[3] += static_cast<V>(b[7]);
0990       t[0] *= z2;
0991       t[1] *= z2;
0992       t[2] *= z2;
0993       t[3] *= z2;
0994       t[0] += static_cast<V>(a[8]);
0995       t[1] += static_cast<V>(a[9]);
0996       t[2] += static_cast<V>(b[8]);
0997       t[3] += static_cast<V>(b[9]);
0998       t[0] *= z2;
0999       t[1] *= z2;
1000       t[2] *= z2;
1001       t[3] *= z2;
1002       t[0] += static_cast<V>(a[10]);
1003       t[1] += static_cast<V>(a[11]);
1004       t[2] += static_cast<V>(b[10]);
1005       t[3] += static_cast<V>(b[11]);
1006       t[0] *= z2;
1007       t[1] *= z2;
1008       t[2] *= z2;
1009       t[3] *= z2;
1010       t[0] += static_cast<V>(a[12]);
1011       t[1] += static_cast<V>(a[13]);
1012       t[2] += static_cast<V>(b[12]);
1013       t[3] += static_cast<V>(b[13]);
1014       t[0] *= z2;
1015       t[1] *= z2;
1016       t[2] *= z2;
1017       t[3] *= z2;
1018       t[0] += static_cast<V>(a[14]);
1019       t[1] += static_cast<V>(a[15]);
1020       t[2] += static_cast<V>(b[14]);
1021       t[3] += static_cast<V>(b[15]);
1022       t[0] *= z;
1023       t[2] *= z;
1024       return (t[0] + t[1]) / (t[2] + t[3]);
1025    }
1026 }
1027
1028 template <class T, class U, class V>
1029 inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const std::integral_constant<int, 17>*) BOOST_MATH_NOEXCEPT(V)
1030 {
1031    if(x <= 1)
1032    {
1033       V x2 = x * x;
1034       V t[4];
1035       t[0] = a[16] * x2 + a[14];
1036       t[1] = a[15] * x2 + a[13];
1037       t[2] = b[16] * x2 + b[14];
1038       t[3] = b[15] * x2 + b[13];
1039       t[0] *= x2;
1040       t[1] *= x2;
1041       t[2] *= x2;
1042       t[3] *= x2;
1043       t[0] += static_cast<V>(a[12]);
1044       t[1] += static_cast<V>(a[11]);
1045       t[2] += static_cast<V>(b[12]);
1046       t[3] += static_cast<V>(b[11]);
1047       t[0] *= x2;
1048       t[1] *= x2;
1049       t[2] *= x2;
1050       t[3] *= x2;
1051       t[0] += static_cast<V>(a[10]);
1052       t[1] += static_cast<V>(a[9]);
1053       t[2] += static_cast<V>(b[10]);
1054       t[3] += static_cast<V>(b[9]);
1055       t[0] *= x2;
1056       t[1] *= x2;
1057       t[2] *= x2;
1058       t[3] *= x2;
1059       t[0] += static_cast<V>(a[8]);
1060       t[1] += static_cast<V>(a[7]);
1061       t[2] += static_cast<V>(b[8]);
1062       t[3] += static_cast<V>(b[7]);
1063       t[0] *= x2;
1064       t[1] *= x2;
1065       t[2] *= x2;
1066       t[3] *= x2;
1067       t[0] += static_cast<V>(a[6]);
1068       t[1] += static_cast<V>(a[5]);
1069       t[2] += static_cast<V>(b[6]);
1070       t[3] += static_cast<V>(b[5]);
1071       t[0] *= x2;
1072       t[1] *= x2;
1073       t[2] *= x2;
1074       t[3] *= x2;
1075       t[0] += static_cast<V>(a[4]);
1076       t[1] += static_cast<V>(a[3]);
1077       t[2] += static_cast<V>(b[4]);
1078       t[3] += static_cast<V>(b[3]);
1079       t[0] *= x2;
1080       t[1] *= x2;
1081       t[2] *= x2;
1082       t[3] *= x2;
1083       t[0] += static_cast<V>(a[2]);
1084       t[1] += static_cast<V>(a[1]);
1085       t[2] += static_cast<V>(b[2]);
1086       t[3] += static_cast<V>(b[1]);
1087       t[0] *= x2;
1088       t[2] *= x2;
1089       t[0] += static_cast<V>(a[0]);
1090       t[2] += static_cast<V>(b[0]);
1091       t[1] *= x;
1092       t[3] *= x;
1093       return (t[0] + t[1]) / (t[2] + t[3]);
1094    }
1095    else
1096    {
1097       V z = 1 / x;
1098       V z2 = 1 / (x * x);
1099       V t[4];
1100       t[0] = a[0] * z2 + a[2];
1101       t[1] = a[1] * z2 + a[3];
1102       t[2] = b[0] * z2 + b[2];
1103       t[3] = b[1] * z2 + b[3];
1104       t[0] *= z2;
1105       t[1] *= z2;
1106       t[2] *= z2;
1107       t[3] *= z2;
1108       t[0] += static_cast<V>(a[4]);
1109       t[1] += static_cast<V>(a[5]);
1110       t[2] += static_cast<V>(b[4]);
1111       t[3] += static_cast<V>(b[5]);
1112       t[0] *= z2;
1113       t[1] *= z2;
1114       t[2] *= z2;
1115       t[3] *= z2;
1116       t[0] += static_cast<V>(a[6]);
1117       t[1] += static_cast<V>(a[7]);
1118       t[2] += static_cast<V>(b[6]);
1119       t[3] += static_cast<V>(b[7]);
1120       t[0] *= z2;
1121       t[1] *= z2;
1122       t[2] *= z2;
1123       t[3] *= z2;
1124       t[0] += static_cast<V>(a[8]);
1125       t[1] += static_cast<V>(a[9]);
1126       t[2] += static_cast<V>(b[8]);
1127       t[3] += static_cast<V>(b[9]);
1128       t[0] *= z2;
1129       t[1] *= z2;
1130       t[2] *= z2;
1131       t[3] *= z2;
1132       t[0] += static_cast<V>(a[10]);
1133       t[1] += static_cast<V>(a[11]);
1134       t[2] += static_cast<V>(b[10]);
1135       t[3] += static_cast<V>(b[11]);
1136       t[0] *= z2;
1137       t[1] *= z2;
1138       t[2] *= z2;
1139       t[3] *= z2;
1140       t[0] += static_cast<V>(a[12]);
1141       t[1] += static_cast<V>(a[13]);
1142       t[2] += static_cast<V>(b[12]);
1143       t[3] += static_cast<V>(b[13]);
1144       t[0] *= z2;
1145       t[1] *= z2;
1146       t[2] *= z2;
1147       t[3] *= z2;
1148       t[0] += static_cast<V>(a[14]);
1149       t[1] += static_cast<V>(a[15]);
1150       t[2] += static_cast<V>(b[14]);
1151       t[3] += static_cast<V>(b[15]);
1152       t[0] *= z2;
1153       t[2] *= z2;
1154       t[0] += static_cast<V>(a[16]);
1155       t[2] += static_cast<V>(b[16]);
1156       t[1] *= z;
1157       t[3] *= z;
1158       return (t[0] + t[1]) / (t[2] + t[3]);
1159    }
1160 }
1161
1162
1163 }}}} // namespaces
1164
1165 #endif // include guard
1166