rational_horner3_19.hpp 37 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458
  1. // (C) Copyright John Maddock 2007.
  2. // Use, modification and distribution are subject to the
  3. // Boost Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  5. //
  6. // This file is machine generated, do not edit by hand
  7. // Polynomial evaluation using second order Horners rule
  8. #ifndef BOOST_MATH_TOOLS_RAT_EVAL_19_HPP
  9. #define BOOST_MATH_TOOLS_RAT_EVAL_19_HPP
  10. namespace boost{ namespace math{ namespace tools{ namespace detail{
  11. template <class T, class U, class V>
  12. inline V evaluate_rational_c_imp(const T*, const U*, const V&, const mpl::int_<0>*) BOOST_MATH_NOEXCEPT(V)
  13. {
  14. return static_cast<V>(0);
  15. }
  16. template <class T, class U, class V>
  17. inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const mpl::int_<1>*) BOOST_MATH_NOEXCEPT(V)
  18. {
  19. return static_cast<V>(a[0]) / static_cast<V>(b[0]);
  20. }
  21. template <class T, class U, class V>
  22. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<2>*) BOOST_MATH_NOEXCEPT(V)
  23. {
  24. return static_cast<V>((a[1] * x + a[0]) / (b[1] * x + b[0]));
  25. }
  26. template <class T, class U, class V>
  27. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<3>*) BOOST_MATH_NOEXCEPT(V)
  28. {
  29. return static_cast<V>(((a[2] * x + a[1]) * x + a[0]) / ((b[2] * x + b[1]) * x + b[0]));
  30. }
  31. template <class T, class U, class V>
  32. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<4>*) BOOST_MATH_NOEXCEPT(V)
  33. {
  34. 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]));
  35. }
  36. template <class T, class U, class V>
  37. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<5>*) BOOST_MATH_NOEXCEPT(V)
  38. {
  39. if(x <= 1)
  40. {
  41. V x2 = x * x;
  42. V t[4];
  43. t[0] = a[4] * x2 + a[2];
  44. t[1] = a[3] * x2 + a[1];
  45. t[2] = b[4] * x2 + b[2];
  46. t[3] = b[3] * x2 + b[1];
  47. t[0] *= x2;
  48. t[2] *= x2;
  49. t[0] += static_cast<V>(a[0]);
  50. t[2] += static_cast<V>(b[0]);
  51. t[1] *= x;
  52. t[3] *= x;
  53. return (t[0] + t[1]) / (t[2] + t[3]);
  54. }
  55. else
  56. {
  57. V z = 1 / x;
  58. V z2 = 1 / (x * x);
  59. V t[4];
  60. t[0] = a[0] * z2 + a[2];
  61. t[1] = a[1] * z2 + a[3];
  62. t[2] = b[0] * z2 + b[2];
  63. t[3] = b[1] * z2 + b[3];
  64. t[0] *= z2;
  65. t[2] *= z2;
  66. t[0] += static_cast<V>(a[4]);
  67. t[2] += static_cast<V>(b[4]);
  68. t[1] *= z;
  69. t[3] *= z;
  70. return (t[0] + t[1]) / (t[2] + t[3]);
  71. }
  72. }
  73. template <class T, class U, class V>
  74. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<6>*) BOOST_MATH_NOEXCEPT(V)
  75. {
  76. if(x <= 1)
  77. {
  78. V x2 = x * x;
  79. V t[4];
  80. t[0] = a[5] * x2 + a[3];
  81. t[1] = a[4] * x2 + a[2];
  82. t[2] = b[5] * x2 + b[3];
  83. t[3] = b[4] * x2 + b[2];
  84. t[0] *= x2;
  85. t[1] *= x2;
  86. t[2] *= x2;
  87. t[3] *= x2;
  88. t[0] += static_cast<V>(a[1]);
  89. t[1] += static_cast<V>(a[0]);
  90. t[2] += static_cast<V>(b[1]);
  91. t[3] += static_cast<V>(b[0]);
  92. t[0] *= x;
  93. t[2] *= x;
  94. return (t[0] + t[1]) / (t[2] + t[3]);
  95. }
  96. else
  97. {
  98. V z = 1 / x;
  99. V z2 = 1 / (x * x);
  100. V t[4];
  101. t[0] = a[0] * z2 + a[2];
  102. t[1] = a[1] * z2 + a[3];
  103. t[2] = b[0] * z2 + b[2];
  104. t[3] = b[1] * z2 + b[3];
  105. t[0] *= z2;
  106. t[1] *= z2;
  107. t[2] *= z2;
  108. t[3] *= z2;
  109. t[0] += static_cast<V>(a[4]);
  110. t[1] += static_cast<V>(a[5]);
  111. t[2] += static_cast<V>(b[4]);
  112. t[3] += static_cast<V>(b[5]);
  113. t[0] *= z;
  114. t[2] *= z;
  115. return (t[0] + t[1]) / (t[2] + t[3]);
  116. }
  117. }
  118. template <class T, class U, class V>
  119. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<7>*) BOOST_MATH_NOEXCEPT(V)
  120. {
  121. if(x <= 1)
  122. {
  123. V x2 = x * x;
  124. V t[4];
  125. t[0] = a[6] * x2 + a[4];
  126. t[1] = a[5] * x2 + a[3];
  127. t[2] = b[6] * x2 + b[4];
  128. t[3] = b[5] * x2 + b[3];
  129. t[0] *= x2;
  130. t[1] *= x2;
  131. t[2] *= x2;
  132. t[3] *= x2;
  133. t[0] += static_cast<V>(a[2]);
  134. t[1] += static_cast<V>(a[1]);
  135. t[2] += static_cast<V>(b[2]);
  136. t[3] += static_cast<V>(b[1]);
  137. t[0] *= x2;
  138. t[2] *= x2;
  139. t[0] += static_cast<V>(a[0]);
  140. t[2] += static_cast<V>(b[0]);
  141. t[1] *= x;
  142. t[3] *= x;
  143. return (t[0] + t[1]) / (t[2] + t[3]);
  144. }
  145. else
  146. {
  147. V z = 1 / x;
  148. V z2 = 1 / (x * x);
  149. V t[4];
  150. t[0] = a[0] * z2 + a[2];
  151. t[1] = a[1] * z2 + a[3];
  152. t[2] = b[0] * z2 + b[2];
  153. t[3] = b[1] * z2 + b[3];
  154. t[0] *= z2;
  155. t[1] *= z2;
  156. t[2] *= z2;
  157. t[3] *= z2;
  158. t[0] += static_cast<V>(a[4]);
  159. t[1] += static_cast<V>(a[5]);
  160. t[2] += static_cast<V>(b[4]);
  161. t[3] += static_cast<V>(b[5]);
  162. t[0] *= z2;
  163. t[2] *= z2;
  164. t[0] += static_cast<V>(a[6]);
  165. t[2] += static_cast<V>(b[6]);
  166. t[1] *= z;
  167. t[3] *= z;
  168. return (t[0] + t[1]) / (t[2] + t[3]);
  169. }
  170. }
  171. template <class T, class U, class V>
  172. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<8>*) BOOST_MATH_NOEXCEPT(V)
  173. {
  174. if(x <= 1)
  175. {
  176. V x2 = x * x;
  177. V t[4];
  178. t[0] = a[7] * x2 + a[5];
  179. t[1] = a[6] * x2 + a[4];
  180. t[2] = b[7] * x2 + b[5];
  181. t[3] = b[6] * x2 + b[4];
  182. t[0] *= x2;
  183. t[1] *= x2;
  184. t[2] *= x2;
  185. t[3] *= x2;
  186. t[0] += static_cast<V>(a[3]);
  187. t[1] += static_cast<V>(a[2]);
  188. t[2] += static_cast<V>(b[3]);
  189. t[3] += static_cast<V>(b[2]);
  190. t[0] *= x2;
  191. t[1] *= x2;
  192. t[2] *= x2;
  193. t[3] *= x2;
  194. t[0] += static_cast<V>(a[1]);
  195. t[1] += static_cast<V>(a[0]);
  196. t[2] += static_cast<V>(b[1]);
  197. t[3] += static_cast<V>(b[0]);
  198. t[0] *= x;
  199. t[2] *= x;
  200. return (t[0] + t[1]) / (t[2] + t[3]);
  201. }
  202. else
  203. {
  204. V z = 1 / x;
  205. V z2 = 1 / (x * x);
  206. V t[4];
  207. t[0] = a[0] * z2 + a[2];
  208. t[1] = a[1] * z2 + a[3];
  209. t[2] = b[0] * z2 + b[2];
  210. t[3] = b[1] * z2 + b[3];
  211. t[0] *= z2;
  212. t[1] *= z2;
  213. t[2] *= z2;
  214. t[3] *= z2;
  215. t[0] += static_cast<V>(a[4]);
  216. t[1] += static_cast<V>(a[5]);
  217. t[2] += static_cast<V>(b[4]);
  218. t[3] += static_cast<V>(b[5]);
  219. t[0] *= z2;
  220. t[1] *= z2;
  221. t[2] *= z2;
  222. t[3] *= z2;
  223. t[0] += static_cast<V>(a[6]);
  224. t[1] += static_cast<V>(a[7]);
  225. t[2] += static_cast<V>(b[6]);
  226. t[3] += static_cast<V>(b[7]);
  227. t[0] *= z;
  228. t[2] *= z;
  229. return (t[0] + t[1]) / (t[2] + t[3]);
  230. }
  231. }
  232. template <class T, class U, class V>
  233. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<9>*) BOOST_MATH_NOEXCEPT(V)
  234. {
  235. if(x <= 1)
  236. {
  237. V x2 = x * x;
  238. V t[4];
  239. t[0] = a[8] * x2 + a[6];
  240. t[1] = a[7] * x2 + a[5];
  241. t[2] = b[8] * x2 + b[6];
  242. t[3] = b[7] * x2 + b[5];
  243. t[0] *= x2;
  244. t[1] *= x2;
  245. t[2] *= x2;
  246. t[3] *= x2;
  247. t[0] += static_cast<V>(a[4]);
  248. t[1] += static_cast<V>(a[3]);
  249. t[2] += static_cast<V>(b[4]);
  250. t[3] += static_cast<V>(b[3]);
  251. t[0] *= x2;
  252. t[1] *= x2;
  253. t[2] *= x2;
  254. t[3] *= x2;
  255. t[0] += static_cast<V>(a[2]);
  256. t[1] += static_cast<V>(a[1]);
  257. t[2] += static_cast<V>(b[2]);
  258. t[3] += static_cast<V>(b[1]);
  259. t[0] *= x2;
  260. t[2] *= x2;
  261. t[0] += static_cast<V>(a[0]);
  262. t[2] += static_cast<V>(b[0]);
  263. t[1] *= x;
  264. t[3] *= x;
  265. return (t[0] + t[1]) / (t[2] + t[3]);
  266. }
  267. else
  268. {
  269. V z = 1 / x;
  270. V z2 = 1 / (x * x);
  271. V t[4];
  272. t[0] = a[0] * z2 + a[2];
  273. t[1] = a[1] * z2 + a[3];
  274. t[2] = b[0] * z2 + b[2];
  275. t[3] = b[1] * z2 + b[3];
  276. t[0] *= z2;
  277. t[1] *= z2;
  278. t[2] *= z2;
  279. t[3] *= z2;
  280. t[0] += static_cast<V>(a[4]);
  281. t[1] += static_cast<V>(a[5]);
  282. t[2] += static_cast<V>(b[4]);
  283. t[3] += static_cast<V>(b[5]);
  284. t[0] *= z2;
  285. t[1] *= z2;
  286. t[2] *= z2;
  287. t[3] *= z2;
  288. t[0] += static_cast<V>(a[6]);
  289. t[1] += static_cast<V>(a[7]);
  290. t[2] += static_cast<V>(b[6]);
  291. t[3] += static_cast<V>(b[7]);
  292. t[0] *= z2;
  293. t[2] *= z2;
  294. t[0] += static_cast<V>(a[8]);
  295. t[2] += static_cast<V>(b[8]);
  296. t[1] *= z;
  297. t[3] *= z;
  298. return (t[0] + t[1]) / (t[2] + t[3]);
  299. }
  300. }
  301. template <class T, class U, class V>
  302. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<10>*) BOOST_MATH_NOEXCEPT(V)
  303. {
  304. if(x <= 1)
  305. {
  306. V x2 = x * x;
  307. V t[4];
  308. t[0] = a[9] * x2 + a[7];
  309. t[1] = a[8] * x2 + a[6];
  310. t[2] = b[9] * x2 + b[7];
  311. t[3] = b[8] * x2 + b[6];
  312. t[0] *= x2;
  313. t[1] *= x2;
  314. t[2] *= x2;
  315. t[3] *= x2;
  316. t[0] += static_cast<V>(a[5]);
  317. t[1] += static_cast<V>(a[4]);
  318. t[2] += static_cast<V>(b[5]);
  319. t[3] += static_cast<V>(b[4]);
  320. t[0] *= x2;
  321. t[1] *= x2;
  322. t[2] *= x2;
  323. t[3] *= x2;
  324. t[0] += static_cast<V>(a[3]);
  325. t[1] += static_cast<V>(a[2]);
  326. t[2] += static_cast<V>(b[3]);
  327. t[3] += static_cast<V>(b[2]);
  328. t[0] *= x2;
  329. t[1] *= x2;
  330. t[2] *= x2;
  331. t[3] *= x2;
  332. t[0] += static_cast<V>(a[1]);
  333. t[1] += static_cast<V>(a[0]);
  334. t[2] += static_cast<V>(b[1]);
  335. t[3] += static_cast<V>(b[0]);
  336. t[0] *= x;
  337. t[2] *= x;
  338. return (t[0] + t[1]) / (t[2] + t[3]);
  339. }
  340. else
  341. {
  342. V z = 1 / x;
  343. V z2 = 1 / (x * x);
  344. V t[4];
  345. t[0] = a[0] * z2 + a[2];
  346. t[1] = a[1] * z2 + a[3];
  347. t[2] = b[0] * z2 + b[2];
  348. t[3] = b[1] * z2 + b[3];
  349. t[0] *= z2;
  350. t[1] *= z2;
  351. t[2] *= z2;
  352. t[3] *= z2;
  353. t[0] += static_cast<V>(a[4]);
  354. t[1] += static_cast<V>(a[5]);
  355. t[2] += static_cast<V>(b[4]);
  356. t[3] += static_cast<V>(b[5]);
  357. t[0] *= z2;
  358. t[1] *= z2;
  359. t[2] *= z2;
  360. t[3] *= z2;
  361. t[0] += static_cast<V>(a[6]);
  362. t[1] += static_cast<V>(a[7]);
  363. t[2] += static_cast<V>(b[6]);
  364. t[3] += static_cast<V>(b[7]);
  365. t[0] *= z2;
  366. t[1] *= z2;
  367. t[2] *= z2;
  368. t[3] *= z2;
  369. t[0] += static_cast<V>(a[8]);
  370. t[1] += static_cast<V>(a[9]);
  371. t[2] += static_cast<V>(b[8]);
  372. t[3] += static_cast<V>(b[9]);
  373. t[0] *= z;
  374. t[2] *= z;
  375. return (t[0] + t[1]) / (t[2] + t[3]);
  376. }
  377. }
  378. template <class T, class U, class V>
  379. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<11>*) BOOST_MATH_NOEXCEPT(V)
  380. {
  381. if(x <= 1)
  382. {
  383. V x2 = x * x;
  384. V t[4];
  385. t[0] = a[10] * x2 + a[8];
  386. t[1] = a[9] * x2 + a[7];
  387. t[2] = b[10] * x2 + b[8];
  388. t[3] = b[9] * x2 + b[7];
  389. t[0] *= x2;
  390. t[1] *= x2;
  391. t[2] *= x2;
  392. t[3] *= x2;
  393. t[0] += static_cast<V>(a[6]);
  394. t[1] += static_cast<V>(a[5]);
  395. t[2] += static_cast<V>(b[6]);
  396. t[3] += static_cast<V>(b[5]);
  397. t[0] *= x2;
  398. t[1] *= x2;
  399. t[2] *= x2;
  400. t[3] *= x2;
  401. t[0] += static_cast<V>(a[4]);
  402. t[1] += static_cast<V>(a[3]);
  403. t[2] += static_cast<V>(b[4]);
  404. t[3] += static_cast<V>(b[3]);
  405. t[0] *= x2;
  406. t[1] *= x2;
  407. t[2] *= x2;
  408. t[3] *= x2;
  409. t[0] += static_cast<V>(a[2]);
  410. t[1] += static_cast<V>(a[1]);
  411. t[2] += static_cast<V>(b[2]);
  412. t[3] += static_cast<V>(b[1]);
  413. t[0] *= x2;
  414. t[2] *= x2;
  415. t[0] += static_cast<V>(a[0]);
  416. t[2] += static_cast<V>(b[0]);
  417. t[1] *= x;
  418. t[3] *= x;
  419. return (t[0] + t[1]) / (t[2] + t[3]);
  420. }
  421. else
  422. {
  423. V z = 1 / x;
  424. V z2 = 1 / (x * x);
  425. V t[4];
  426. t[0] = a[0] * z2 + a[2];
  427. t[1] = a[1] * z2 + a[3];
  428. t[2] = b[0] * z2 + b[2];
  429. t[3] = b[1] * z2 + b[3];
  430. t[0] *= z2;
  431. t[1] *= z2;
  432. t[2] *= z2;
  433. t[3] *= z2;
  434. t[0] += static_cast<V>(a[4]);
  435. t[1] += static_cast<V>(a[5]);
  436. t[2] += static_cast<V>(b[4]);
  437. t[3] += static_cast<V>(b[5]);
  438. t[0] *= z2;
  439. t[1] *= z2;
  440. t[2] *= z2;
  441. t[3] *= z2;
  442. t[0] += static_cast<V>(a[6]);
  443. t[1] += static_cast<V>(a[7]);
  444. t[2] += static_cast<V>(b[6]);
  445. t[3] += static_cast<V>(b[7]);
  446. t[0] *= z2;
  447. t[1] *= z2;
  448. t[2] *= z2;
  449. t[3] *= z2;
  450. t[0] += static_cast<V>(a[8]);
  451. t[1] += static_cast<V>(a[9]);
  452. t[2] += static_cast<V>(b[8]);
  453. t[3] += static_cast<V>(b[9]);
  454. t[0] *= z2;
  455. t[2] *= z2;
  456. t[0] += static_cast<V>(a[10]);
  457. t[2] += static_cast<V>(b[10]);
  458. t[1] *= z;
  459. t[3] *= z;
  460. return (t[0] + t[1]) / (t[2] + t[3]);
  461. }
  462. }
  463. template <class T, class U, class V>
  464. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<12>*) BOOST_MATH_NOEXCEPT(V)
  465. {
  466. if(x <= 1)
  467. {
  468. V x2 = x * x;
  469. V t[4];
  470. t[0] = a[11] * x2 + a[9];
  471. t[1] = a[10] * x2 + a[8];
  472. t[2] = b[11] * x2 + b[9];
  473. t[3] = b[10] * x2 + b[8];
  474. t[0] *= x2;
  475. t[1] *= x2;
  476. t[2] *= x2;
  477. t[3] *= x2;
  478. t[0] += static_cast<V>(a[7]);
  479. t[1] += static_cast<V>(a[6]);
  480. t[2] += static_cast<V>(b[7]);
  481. t[3] += static_cast<V>(b[6]);
  482. t[0] *= x2;
  483. t[1] *= x2;
  484. t[2] *= x2;
  485. t[3] *= x2;
  486. t[0] += static_cast<V>(a[5]);
  487. t[1] += static_cast<V>(a[4]);
  488. t[2] += static_cast<V>(b[5]);
  489. t[3] += static_cast<V>(b[4]);
  490. t[0] *= x2;
  491. t[1] *= x2;
  492. t[2] *= x2;
  493. t[3] *= x2;
  494. t[0] += static_cast<V>(a[3]);
  495. t[1] += static_cast<V>(a[2]);
  496. t[2] += static_cast<V>(b[3]);
  497. t[3] += static_cast<V>(b[2]);
  498. t[0] *= x2;
  499. t[1] *= x2;
  500. t[2] *= x2;
  501. t[3] *= x2;
  502. t[0] += static_cast<V>(a[1]);
  503. t[1] += static_cast<V>(a[0]);
  504. t[2] += static_cast<V>(b[1]);
  505. t[3] += static_cast<V>(b[0]);
  506. t[0] *= x;
  507. t[2] *= x;
  508. return (t[0] + t[1]) / (t[2] + t[3]);
  509. }
  510. else
  511. {
  512. V z = 1 / x;
  513. V z2 = 1 / (x * x);
  514. V t[4];
  515. t[0] = a[0] * z2 + a[2];
  516. t[1] = a[1] * z2 + a[3];
  517. t[2] = b[0] * z2 + b[2];
  518. t[3] = b[1] * z2 + b[3];
  519. t[0] *= z2;
  520. t[1] *= z2;
  521. t[2] *= z2;
  522. t[3] *= z2;
  523. t[0] += static_cast<V>(a[4]);
  524. t[1] += static_cast<V>(a[5]);
  525. t[2] += static_cast<V>(b[4]);
  526. t[3] += static_cast<V>(b[5]);
  527. t[0] *= z2;
  528. t[1] *= z2;
  529. t[2] *= z2;
  530. t[3] *= z2;
  531. t[0] += static_cast<V>(a[6]);
  532. t[1] += static_cast<V>(a[7]);
  533. t[2] += static_cast<V>(b[6]);
  534. t[3] += static_cast<V>(b[7]);
  535. t[0] *= z2;
  536. t[1] *= z2;
  537. t[2] *= z2;
  538. t[3] *= z2;
  539. t[0] += static_cast<V>(a[8]);
  540. t[1] += static_cast<V>(a[9]);
  541. t[2] += static_cast<V>(b[8]);
  542. t[3] += static_cast<V>(b[9]);
  543. t[0] *= z2;
  544. t[1] *= z2;
  545. t[2] *= z2;
  546. t[3] *= z2;
  547. t[0] += static_cast<V>(a[10]);
  548. t[1] += static_cast<V>(a[11]);
  549. t[2] += static_cast<V>(b[10]);
  550. t[3] += static_cast<V>(b[11]);
  551. t[0] *= z;
  552. t[2] *= z;
  553. return (t[0] + t[1]) / (t[2] + t[3]);
  554. }
  555. }
  556. template <class T, class U, class V>
  557. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<13>*) BOOST_MATH_NOEXCEPT(V)
  558. {
  559. if(x <= 1)
  560. {
  561. V x2 = x * x;
  562. V t[4];
  563. t[0] = a[12] * x2 + a[10];
  564. t[1] = a[11] * x2 + a[9];
  565. t[2] = b[12] * x2 + b[10];
  566. t[3] = b[11] * x2 + b[9];
  567. t[0] *= x2;
  568. t[1] *= x2;
  569. t[2] *= x2;
  570. t[3] *= x2;
  571. t[0] += static_cast<V>(a[8]);
  572. t[1] += static_cast<V>(a[7]);
  573. t[2] += static_cast<V>(b[8]);
  574. t[3] += static_cast<V>(b[7]);
  575. t[0] *= x2;
  576. t[1] *= x2;
  577. t[2] *= x2;
  578. t[3] *= x2;
  579. t[0] += static_cast<V>(a[6]);
  580. t[1] += static_cast<V>(a[5]);
  581. t[2] += static_cast<V>(b[6]);
  582. t[3] += static_cast<V>(b[5]);
  583. t[0] *= x2;
  584. t[1] *= x2;
  585. t[2] *= x2;
  586. t[3] *= x2;
  587. t[0] += static_cast<V>(a[4]);
  588. t[1] += static_cast<V>(a[3]);
  589. t[2] += static_cast<V>(b[4]);
  590. t[3] += static_cast<V>(b[3]);
  591. t[0] *= x2;
  592. t[1] *= x2;
  593. t[2] *= x2;
  594. t[3] *= x2;
  595. t[0] += static_cast<V>(a[2]);
  596. t[1] += static_cast<V>(a[1]);
  597. t[2] += static_cast<V>(b[2]);
  598. t[3] += static_cast<V>(b[1]);
  599. t[0] *= x2;
  600. t[2] *= x2;
  601. t[0] += static_cast<V>(a[0]);
  602. t[2] += static_cast<V>(b[0]);
  603. t[1] *= x;
  604. t[3] *= x;
  605. return (t[0] + t[1]) / (t[2] + t[3]);
  606. }
  607. else
  608. {
  609. V z = 1 / x;
  610. V z2 = 1 / (x * x);
  611. V t[4];
  612. t[0] = a[0] * z2 + a[2];
  613. t[1] = a[1] * z2 + a[3];
  614. t[2] = b[0] * z2 + b[2];
  615. t[3] = b[1] * z2 + b[3];
  616. t[0] *= z2;
  617. t[1] *= z2;
  618. t[2] *= z2;
  619. t[3] *= z2;
  620. t[0] += static_cast<V>(a[4]);
  621. t[1] += static_cast<V>(a[5]);
  622. t[2] += static_cast<V>(b[4]);
  623. t[3] += static_cast<V>(b[5]);
  624. t[0] *= z2;
  625. t[1] *= z2;
  626. t[2] *= z2;
  627. t[3] *= z2;
  628. t[0] += static_cast<V>(a[6]);
  629. t[1] += static_cast<V>(a[7]);
  630. t[2] += static_cast<V>(b[6]);
  631. t[3] += static_cast<V>(b[7]);
  632. t[0] *= z2;
  633. t[1] *= z2;
  634. t[2] *= z2;
  635. t[3] *= z2;
  636. t[0] += static_cast<V>(a[8]);
  637. t[1] += static_cast<V>(a[9]);
  638. t[2] += static_cast<V>(b[8]);
  639. t[3] += static_cast<V>(b[9]);
  640. t[0] *= z2;
  641. t[1] *= z2;
  642. t[2] *= z2;
  643. t[3] *= z2;
  644. t[0] += static_cast<V>(a[10]);
  645. t[1] += static_cast<V>(a[11]);
  646. t[2] += static_cast<V>(b[10]);
  647. t[3] += static_cast<V>(b[11]);
  648. t[0] *= z2;
  649. t[2] *= z2;
  650. t[0] += static_cast<V>(a[12]);
  651. t[2] += static_cast<V>(b[12]);
  652. t[1] *= z;
  653. t[3] *= z;
  654. return (t[0] + t[1]) / (t[2] + t[3]);
  655. }
  656. }
  657. template <class T, class U, class V>
  658. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<14>*) BOOST_MATH_NOEXCEPT(V)
  659. {
  660. if(x <= 1)
  661. {
  662. V x2 = x * x;
  663. V t[4];
  664. t[0] = a[13] * x2 + a[11];
  665. t[1] = a[12] * x2 + a[10];
  666. t[2] = b[13] * x2 + b[11];
  667. t[3] = b[12] * x2 + b[10];
  668. t[0] *= x2;
  669. t[1] *= x2;
  670. t[2] *= x2;
  671. t[3] *= x2;
  672. t[0] += static_cast<V>(a[9]);
  673. t[1] += static_cast<V>(a[8]);
  674. t[2] += static_cast<V>(b[9]);
  675. t[3] += static_cast<V>(b[8]);
  676. t[0] *= x2;
  677. t[1] *= x2;
  678. t[2] *= x2;
  679. t[3] *= x2;
  680. t[0] += static_cast<V>(a[7]);
  681. t[1] += static_cast<V>(a[6]);
  682. t[2] += static_cast<V>(b[7]);
  683. t[3] += static_cast<V>(b[6]);
  684. t[0] *= x2;
  685. t[1] *= x2;
  686. t[2] *= x2;
  687. t[3] *= x2;
  688. t[0] += static_cast<V>(a[5]);
  689. t[1] += static_cast<V>(a[4]);
  690. t[2] += static_cast<V>(b[5]);
  691. t[3] += static_cast<V>(b[4]);
  692. t[0] *= x2;
  693. t[1] *= x2;
  694. t[2] *= x2;
  695. t[3] *= x2;
  696. t[0] += static_cast<V>(a[3]);
  697. t[1] += static_cast<V>(a[2]);
  698. t[2] += static_cast<V>(b[3]);
  699. t[3] += static_cast<V>(b[2]);
  700. t[0] *= x2;
  701. t[1] *= x2;
  702. t[2] *= x2;
  703. t[3] *= x2;
  704. t[0] += static_cast<V>(a[1]);
  705. t[1] += static_cast<V>(a[0]);
  706. t[2] += static_cast<V>(b[1]);
  707. t[3] += static_cast<V>(b[0]);
  708. t[0] *= x;
  709. t[2] *= x;
  710. return (t[0] + t[1]) / (t[2] + t[3]);
  711. }
  712. else
  713. {
  714. V z = 1 / x;
  715. V z2 = 1 / (x * x);
  716. V t[4];
  717. t[0] = a[0] * z2 + a[2];
  718. t[1] = a[1] * z2 + a[3];
  719. t[2] = b[0] * z2 + b[2];
  720. t[3] = b[1] * z2 + b[3];
  721. t[0] *= z2;
  722. t[1] *= z2;
  723. t[2] *= z2;
  724. t[3] *= z2;
  725. t[0] += static_cast<V>(a[4]);
  726. t[1] += static_cast<V>(a[5]);
  727. t[2] += static_cast<V>(b[4]);
  728. t[3] += static_cast<V>(b[5]);
  729. t[0] *= z2;
  730. t[1] *= z2;
  731. t[2] *= z2;
  732. t[3] *= z2;
  733. t[0] += static_cast<V>(a[6]);
  734. t[1] += static_cast<V>(a[7]);
  735. t[2] += static_cast<V>(b[6]);
  736. t[3] += static_cast<V>(b[7]);
  737. t[0] *= z2;
  738. t[1] *= z2;
  739. t[2] *= z2;
  740. t[3] *= z2;
  741. t[0] += static_cast<V>(a[8]);
  742. t[1] += static_cast<V>(a[9]);
  743. t[2] += static_cast<V>(b[8]);
  744. t[3] += static_cast<V>(b[9]);
  745. t[0] *= z2;
  746. t[1] *= z2;
  747. t[2] *= z2;
  748. t[3] *= z2;
  749. t[0] += static_cast<V>(a[10]);
  750. t[1] += static_cast<V>(a[11]);
  751. t[2] += static_cast<V>(b[10]);
  752. t[3] += static_cast<V>(b[11]);
  753. t[0] *= z2;
  754. t[1] *= z2;
  755. t[2] *= z2;
  756. t[3] *= z2;
  757. t[0] += static_cast<V>(a[12]);
  758. t[1] += static_cast<V>(a[13]);
  759. t[2] += static_cast<V>(b[12]);
  760. t[3] += static_cast<V>(b[13]);
  761. t[0] *= z;
  762. t[2] *= z;
  763. return (t[0] + t[1]) / (t[2] + t[3]);
  764. }
  765. }
  766. template <class T, class U, class V>
  767. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<15>*) BOOST_MATH_NOEXCEPT(V)
  768. {
  769. if(x <= 1)
  770. {
  771. V x2 = x * x;
  772. V t[4];
  773. t[0] = a[14] * x2 + a[12];
  774. t[1] = a[13] * x2 + a[11];
  775. t[2] = b[14] * x2 + b[12];
  776. t[3] = b[13] * x2 + b[11];
  777. t[0] *= x2;
  778. t[1] *= x2;
  779. t[2] *= x2;
  780. t[3] *= x2;
  781. t[0] += static_cast<V>(a[10]);
  782. t[1] += static_cast<V>(a[9]);
  783. t[2] += static_cast<V>(b[10]);
  784. t[3] += static_cast<V>(b[9]);
  785. t[0] *= x2;
  786. t[1] *= x2;
  787. t[2] *= x2;
  788. t[3] *= x2;
  789. t[0] += static_cast<V>(a[8]);
  790. t[1] += static_cast<V>(a[7]);
  791. t[2] += static_cast<V>(b[8]);
  792. t[3] += static_cast<V>(b[7]);
  793. t[0] *= x2;
  794. t[1] *= x2;
  795. t[2] *= x2;
  796. t[3] *= x2;
  797. t[0] += static_cast<V>(a[6]);
  798. t[1] += static_cast<V>(a[5]);
  799. t[2] += static_cast<V>(b[6]);
  800. t[3] += static_cast<V>(b[5]);
  801. t[0] *= x2;
  802. t[1] *= x2;
  803. t[2] *= x2;
  804. t[3] *= x2;
  805. t[0] += static_cast<V>(a[4]);
  806. t[1] += static_cast<V>(a[3]);
  807. t[2] += static_cast<V>(b[4]);
  808. t[3] += static_cast<V>(b[3]);
  809. t[0] *= x2;
  810. t[1] *= x2;
  811. t[2] *= x2;
  812. t[3] *= x2;
  813. t[0] += static_cast<V>(a[2]);
  814. t[1] += static_cast<V>(a[1]);
  815. t[2] += static_cast<V>(b[2]);
  816. t[3] += static_cast<V>(b[1]);
  817. t[0] *= x2;
  818. t[2] *= x2;
  819. t[0] += static_cast<V>(a[0]);
  820. t[2] += static_cast<V>(b[0]);
  821. t[1] *= x;
  822. t[3] *= x;
  823. return (t[0] + t[1]) / (t[2] + t[3]);
  824. }
  825. else
  826. {
  827. V z = 1 / x;
  828. V z2 = 1 / (x * x);
  829. V t[4];
  830. t[0] = a[0] * z2 + a[2];
  831. t[1] = a[1] * z2 + a[3];
  832. t[2] = b[0] * z2 + b[2];
  833. t[3] = b[1] * z2 + b[3];
  834. t[0] *= z2;
  835. t[1] *= z2;
  836. t[2] *= z2;
  837. t[3] *= z2;
  838. t[0] += static_cast<V>(a[4]);
  839. t[1] += static_cast<V>(a[5]);
  840. t[2] += static_cast<V>(b[4]);
  841. t[3] += static_cast<V>(b[5]);
  842. t[0] *= z2;
  843. t[1] *= z2;
  844. t[2] *= z2;
  845. t[3] *= z2;
  846. t[0] += static_cast<V>(a[6]);
  847. t[1] += static_cast<V>(a[7]);
  848. t[2] += static_cast<V>(b[6]);
  849. t[3] += static_cast<V>(b[7]);
  850. t[0] *= z2;
  851. t[1] *= z2;
  852. t[2] *= z2;
  853. t[3] *= z2;
  854. t[0] += static_cast<V>(a[8]);
  855. t[1] += static_cast<V>(a[9]);
  856. t[2] += static_cast<V>(b[8]);
  857. t[3] += static_cast<V>(b[9]);
  858. t[0] *= z2;
  859. t[1] *= z2;
  860. t[2] *= z2;
  861. t[3] *= z2;
  862. t[0] += static_cast<V>(a[10]);
  863. t[1] += static_cast<V>(a[11]);
  864. t[2] += static_cast<V>(b[10]);
  865. t[3] += static_cast<V>(b[11]);
  866. t[0] *= z2;
  867. t[1] *= z2;
  868. t[2] *= z2;
  869. t[3] *= z2;
  870. t[0] += static_cast<V>(a[12]);
  871. t[1] += static_cast<V>(a[13]);
  872. t[2] += static_cast<V>(b[12]);
  873. t[3] += static_cast<V>(b[13]);
  874. t[0] *= z2;
  875. t[2] *= z2;
  876. t[0] += static_cast<V>(a[14]);
  877. t[2] += static_cast<V>(b[14]);
  878. t[1] *= z;
  879. t[3] *= z;
  880. return (t[0] + t[1]) / (t[2] + t[3]);
  881. }
  882. }
  883. template <class T, class U, class V>
  884. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<16>*) BOOST_MATH_NOEXCEPT(V)
  885. {
  886. if(x <= 1)
  887. {
  888. V x2 = x * x;
  889. V t[4];
  890. t[0] = a[15] * x2 + a[13];
  891. t[1] = a[14] * x2 + a[12];
  892. t[2] = b[15] * x2 + b[13];
  893. t[3] = b[14] * x2 + b[12];
  894. t[0] *= x2;
  895. t[1] *= x2;
  896. t[2] *= x2;
  897. t[3] *= x2;
  898. t[0] += static_cast<V>(a[11]);
  899. t[1] += static_cast<V>(a[10]);
  900. t[2] += static_cast<V>(b[11]);
  901. t[3] += static_cast<V>(b[10]);
  902. t[0] *= x2;
  903. t[1] *= x2;
  904. t[2] *= x2;
  905. t[3] *= x2;
  906. t[0] += static_cast<V>(a[9]);
  907. t[1] += static_cast<V>(a[8]);
  908. t[2] += static_cast<V>(b[9]);
  909. t[3] += static_cast<V>(b[8]);
  910. t[0] *= x2;
  911. t[1] *= x2;
  912. t[2] *= x2;
  913. t[3] *= x2;
  914. t[0] += static_cast<V>(a[7]);
  915. t[1] += static_cast<V>(a[6]);
  916. t[2] += static_cast<V>(b[7]);
  917. t[3] += static_cast<V>(b[6]);
  918. t[0] *= x2;
  919. t[1] *= x2;
  920. t[2] *= x2;
  921. t[3] *= x2;
  922. t[0] += static_cast<V>(a[5]);
  923. t[1] += static_cast<V>(a[4]);
  924. t[2] += static_cast<V>(b[5]);
  925. t[3] += static_cast<V>(b[4]);
  926. t[0] *= x2;
  927. t[1] *= x2;
  928. t[2] *= x2;
  929. t[3] *= x2;
  930. t[0] += static_cast<V>(a[3]);
  931. t[1] += static_cast<V>(a[2]);
  932. t[2] += static_cast<V>(b[3]);
  933. t[3] += static_cast<V>(b[2]);
  934. t[0] *= x2;
  935. t[1] *= x2;
  936. t[2] *= x2;
  937. t[3] *= x2;
  938. t[0] += static_cast<V>(a[1]);
  939. t[1] += static_cast<V>(a[0]);
  940. t[2] += static_cast<V>(b[1]);
  941. t[3] += static_cast<V>(b[0]);
  942. t[0] *= x;
  943. t[2] *= x;
  944. return (t[0] + t[1]) / (t[2] + t[3]);
  945. }
  946. else
  947. {
  948. V z = 1 / x;
  949. V z2 = 1 / (x * x);
  950. V t[4];
  951. t[0] = a[0] * z2 + a[2];
  952. t[1] = a[1] * z2 + a[3];
  953. t[2] = b[0] * z2 + b[2];
  954. t[3] = b[1] * z2 + b[3];
  955. t[0] *= z2;
  956. t[1] *= z2;
  957. t[2] *= z2;
  958. t[3] *= z2;
  959. t[0] += static_cast<V>(a[4]);
  960. t[1] += static_cast<V>(a[5]);
  961. t[2] += static_cast<V>(b[4]);
  962. t[3] += static_cast<V>(b[5]);
  963. t[0] *= z2;
  964. t[1] *= z2;
  965. t[2] *= z2;
  966. t[3] *= z2;
  967. t[0] += static_cast<V>(a[6]);
  968. t[1] += static_cast<V>(a[7]);
  969. t[2] += static_cast<V>(b[6]);
  970. t[3] += static_cast<V>(b[7]);
  971. t[0] *= z2;
  972. t[1] *= z2;
  973. t[2] *= z2;
  974. t[3] *= z2;
  975. t[0] += static_cast<V>(a[8]);
  976. t[1] += static_cast<V>(a[9]);
  977. t[2] += static_cast<V>(b[8]);
  978. t[3] += static_cast<V>(b[9]);
  979. t[0] *= z2;
  980. t[1] *= z2;
  981. t[2] *= z2;
  982. t[3] *= z2;
  983. t[0] += static_cast<V>(a[10]);
  984. t[1] += static_cast<V>(a[11]);
  985. t[2] += static_cast<V>(b[10]);
  986. t[3] += static_cast<V>(b[11]);
  987. t[0] *= z2;
  988. t[1] *= z2;
  989. t[2] *= z2;
  990. t[3] *= z2;
  991. t[0] += static_cast<V>(a[12]);
  992. t[1] += static_cast<V>(a[13]);
  993. t[2] += static_cast<V>(b[12]);
  994. t[3] += static_cast<V>(b[13]);
  995. t[0] *= z2;
  996. t[1] *= z2;
  997. t[2] *= z2;
  998. t[3] *= z2;
  999. t[0] += static_cast<V>(a[14]);
  1000. t[1] += static_cast<V>(a[15]);
  1001. t[2] += static_cast<V>(b[14]);
  1002. t[3] += static_cast<V>(b[15]);
  1003. t[0] *= z;
  1004. t[2] *= z;
  1005. return (t[0] + t[1]) / (t[2] + t[3]);
  1006. }
  1007. }
  1008. template <class T, class U, class V>
  1009. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<17>*) BOOST_MATH_NOEXCEPT(V)
  1010. {
  1011. if(x <= 1)
  1012. {
  1013. V x2 = x * x;
  1014. V t[4];
  1015. t[0] = a[16] * x2 + a[14];
  1016. t[1] = a[15] * x2 + a[13];
  1017. t[2] = b[16] * x2 + b[14];
  1018. t[3] = b[15] * x2 + b[13];
  1019. t[0] *= x2;
  1020. t[1] *= x2;
  1021. t[2] *= x2;
  1022. t[3] *= x2;
  1023. t[0] += static_cast<V>(a[12]);
  1024. t[1] += static_cast<V>(a[11]);
  1025. t[2] += static_cast<V>(b[12]);
  1026. t[3] += static_cast<V>(b[11]);
  1027. t[0] *= x2;
  1028. t[1] *= x2;
  1029. t[2] *= x2;
  1030. t[3] *= x2;
  1031. t[0] += static_cast<V>(a[10]);
  1032. t[1] += static_cast<V>(a[9]);
  1033. t[2] += static_cast<V>(b[10]);
  1034. t[3] += static_cast<V>(b[9]);
  1035. t[0] *= x2;
  1036. t[1] *= x2;
  1037. t[2] *= x2;
  1038. t[3] *= x2;
  1039. t[0] += static_cast<V>(a[8]);
  1040. t[1] += static_cast<V>(a[7]);
  1041. t[2] += static_cast<V>(b[8]);
  1042. t[3] += static_cast<V>(b[7]);
  1043. t[0] *= x2;
  1044. t[1] *= x2;
  1045. t[2] *= x2;
  1046. t[3] *= x2;
  1047. t[0] += static_cast<V>(a[6]);
  1048. t[1] += static_cast<V>(a[5]);
  1049. t[2] += static_cast<V>(b[6]);
  1050. t[3] += static_cast<V>(b[5]);
  1051. t[0] *= x2;
  1052. t[1] *= x2;
  1053. t[2] *= x2;
  1054. t[3] *= x2;
  1055. t[0] += static_cast<V>(a[4]);
  1056. t[1] += static_cast<V>(a[3]);
  1057. t[2] += static_cast<V>(b[4]);
  1058. t[3] += static_cast<V>(b[3]);
  1059. t[0] *= x2;
  1060. t[1] *= x2;
  1061. t[2] *= x2;
  1062. t[3] *= x2;
  1063. t[0] += static_cast<V>(a[2]);
  1064. t[1] += static_cast<V>(a[1]);
  1065. t[2] += static_cast<V>(b[2]);
  1066. t[3] += static_cast<V>(b[1]);
  1067. t[0] *= x2;
  1068. t[2] *= x2;
  1069. t[0] += static_cast<V>(a[0]);
  1070. t[2] += static_cast<V>(b[0]);
  1071. t[1] *= x;
  1072. t[3] *= x;
  1073. return (t[0] + t[1]) / (t[2] + t[3]);
  1074. }
  1075. else
  1076. {
  1077. V z = 1 / x;
  1078. V z2 = 1 / (x * x);
  1079. V t[4];
  1080. t[0] = a[0] * z2 + a[2];
  1081. t[1] = a[1] * z2 + a[3];
  1082. t[2] = b[0] * z2 + b[2];
  1083. t[3] = b[1] * z2 + b[3];
  1084. t[0] *= z2;
  1085. t[1] *= z2;
  1086. t[2] *= z2;
  1087. t[3] *= z2;
  1088. t[0] += static_cast<V>(a[4]);
  1089. t[1] += static_cast<V>(a[5]);
  1090. t[2] += static_cast<V>(b[4]);
  1091. t[3] += static_cast<V>(b[5]);
  1092. t[0] *= z2;
  1093. t[1] *= z2;
  1094. t[2] *= z2;
  1095. t[3] *= z2;
  1096. t[0] += static_cast<V>(a[6]);
  1097. t[1] += static_cast<V>(a[7]);
  1098. t[2] += static_cast<V>(b[6]);
  1099. t[3] += static_cast<V>(b[7]);
  1100. t[0] *= z2;
  1101. t[1] *= z2;
  1102. t[2] *= z2;
  1103. t[3] *= z2;
  1104. t[0] += static_cast<V>(a[8]);
  1105. t[1] += static_cast<V>(a[9]);
  1106. t[2] += static_cast<V>(b[8]);
  1107. t[3] += static_cast<V>(b[9]);
  1108. t[0] *= z2;
  1109. t[1] *= z2;
  1110. t[2] *= z2;
  1111. t[3] *= z2;
  1112. t[0] += static_cast<V>(a[10]);
  1113. t[1] += static_cast<V>(a[11]);
  1114. t[2] += static_cast<V>(b[10]);
  1115. t[3] += static_cast<V>(b[11]);
  1116. t[0] *= z2;
  1117. t[1] *= z2;
  1118. t[2] *= z2;
  1119. t[3] *= z2;
  1120. t[0] += static_cast<V>(a[12]);
  1121. t[1] += static_cast<V>(a[13]);
  1122. t[2] += static_cast<V>(b[12]);
  1123. t[3] += static_cast<V>(b[13]);
  1124. t[0] *= z2;
  1125. t[1] *= z2;
  1126. t[2] *= z2;
  1127. t[3] *= z2;
  1128. t[0] += static_cast<V>(a[14]);
  1129. t[1] += static_cast<V>(a[15]);
  1130. t[2] += static_cast<V>(b[14]);
  1131. t[3] += static_cast<V>(b[15]);
  1132. t[0] *= z2;
  1133. t[2] *= z2;
  1134. t[0] += static_cast<V>(a[16]);
  1135. t[2] += static_cast<V>(b[16]);
  1136. t[1] *= z;
  1137. t[3] *= z;
  1138. return (t[0] + t[1]) / (t[2] + t[3]);
  1139. }
  1140. }
  1141. template <class T, class U, class V>
  1142. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<18>*) BOOST_MATH_NOEXCEPT(V)
  1143. {
  1144. if(x <= 1)
  1145. {
  1146. V x2 = x * x;
  1147. V t[4];
  1148. t[0] = a[17] * x2 + a[15];
  1149. t[1] = a[16] * x2 + a[14];
  1150. t[2] = b[17] * x2 + b[15];
  1151. t[3] = b[16] * x2 + b[14];
  1152. t[0] *= x2;
  1153. t[1] *= x2;
  1154. t[2] *= x2;
  1155. t[3] *= x2;
  1156. t[0] += static_cast<V>(a[13]);
  1157. t[1] += static_cast<V>(a[12]);
  1158. t[2] += static_cast<V>(b[13]);
  1159. t[3] += static_cast<V>(b[12]);
  1160. t[0] *= x2;
  1161. t[1] *= x2;
  1162. t[2] *= x2;
  1163. t[3] *= x2;
  1164. t[0] += static_cast<V>(a[11]);
  1165. t[1] += static_cast<V>(a[10]);
  1166. t[2] += static_cast<V>(b[11]);
  1167. t[3] += static_cast<V>(b[10]);
  1168. t[0] *= x2;
  1169. t[1] *= x2;
  1170. t[2] *= x2;
  1171. t[3] *= x2;
  1172. t[0] += static_cast<V>(a[9]);
  1173. t[1] += static_cast<V>(a[8]);
  1174. t[2] += static_cast<V>(b[9]);
  1175. t[3] += static_cast<V>(b[8]);
  1176. t[0] *= x2;
  1177. t[1] *= x2;
  1178. t[2] *= x2;
  1179. t[3] *= x2;
  1180. t[0] += static_cast<V>(a[7]);
  1181. t[1] += static_cast<V>(a[6]);
  1182. t[2] += static_cast<V>(b[7]);
  1183. t[3] += static_cast<V>(b[6]);
  1184. t[0] *= x2;
  1185. t[1] *= x2;
  1186. t[2] *= x2;
  1187. t[3] *= x2;
  1188. t[0] += static_cast<V>(a[5]);
  1189. t[1] += static_cast<V>(a[4]);
  1190. t[2] += static_cast<V>(b[5]);
  1191. t[3] += static_cast<V>(b[4]);
  1192. t[0] *= x2;
  1193. t[1] *= x2;
  1194. t[2] *= x2;
  1195. t[3] *= x2;
  1196. t[0] += static_cast<V>(a[3]);
  1197. t[1] += static_cast<V>(a[2]);
  1198. t[2] += static_cast<V>(b[3]);
  1199. t[3] += static_cast<V>(b[2]);
  1200. t[0] *= x2;
  1201. t[1] *= x2;
  1202. t[2] *= x2;
  1203. t[3] *= x2;
  1204. t[0] += static_cast<V>(a[1]);
  1205. t[1] += static_cast<V>(a[0]);
  1206. t[2] += static_cast<V>(b[1]);
  1207. t[3] += static_cast<V>(b[0]);
  1208. t[0] *= x;
  1209. t[2] *= x;
  1210. return (t[0] + t[1]) / (t[2] + t[3]);
  1211. }
  1212. else
  1213. {
  1214. V z = 1 / x;
  1215. V z2 = 1 / (x * x);
  1216. V t[4];
  1217. t[0] = a[0] * z2 + a[2];
  1218. t[1] = a[1] * z2 + a[3];
  1219. t[2] = b[0] * z2 + b[2];
  1220. t[3] = b[1] * z2 + b[3];
  1221. t[0] *= z2;
  1222. t[1] *= z2;
  1223. t[2] *= z2;
  1224. t[3] *= z2;
  1225. t[0] += static_cast<V>(a[4]);
  1226. t[1] += static_cast<V>(a[5]);
  1227. t[2] += static_cast<V>(b[4]);
  1228. t[3] += static_cast<V>(b[5]);
  1229. t[0] *= z2;
  1230. t[1] *= z2;
  1231. t[2] *= z2;
  1232. t[3] *= z2;
  1233. t[0] += static_cast<V>(a[6]);
  1234. t[1] += static_cast<V>(a[7]);
  1235. t[2] += static_cast<V>(b[6]);
  1236. t[3] += static_cast<V>(b[7]);
  1237. t[0] *= z2;
  1238. t[1] *= z2;
  1239. t[2] *= z2;
  1240. t[3] *= z2;
  1241. t[0] += static_cast<V>(a[8]);
  1242. t[1] += static_cast<V>(a[9]);
  1243. t[2] += static_cast<V>(b[8]);
  1244. t[3] += static_cast<V>(b[9]);
  1245. t[0] *= z2;
  1246. t[1] *= z2;
  1247. t[2] *= z2;
  1248. t[3] *= z2;
  1249. t[0] += static_cast<V>(a[10]);
  1250. t[1] += static_cast<V>(a[11]);
  1251. t[2] += static_cast<V>(b[10]);
  1252. t[3] += static_cast<V>(b[11]);
  1253. t[0] *= z2;
  1254. t[1] *= z2;
  1255. t[2] *= z2;
  1256. t[3] *= z2;
  1257. t[0] += static_cast<V>(a[12]);
  1258. t[1] += static_cast<V>(a[13]);
  1259. t[2] += static_cast<V>(b[12]);
  1260. t[3] += static_cast<V>(b[13]);
  1261. t[0] *= z2;
  1262. t[1] *= z2;
  1263. t[2] *= z2;
  1264. t[3] *= z2;
  1265. t[0] += static_cast<V>(a[14]);
  1266. t[1] += static_cast<V>(a[15]);
  1267. t[2] += static_cast<V>(b[14]);
  1268. t[3] += static_cast<V>(b[15]);
  1269. t[0] *= z2;
  1270. t[1] *= z2;
  1271. t[2] *= z2;
  1272. t[3] *= z2;
  1273. t[0] += static_cast<V>(a[16]);
  1274. t[1] += static_cast<V>(a[17]);
  1275. t[2] += static_cast<V>(b[16]);
  1276. t[3] += static_cast<V>(b[17]);
  1277. t[0] *= z;
  1278. t[2] *= z;
  1279. return (t[0] + t[1]) / (t[2] + t[3]);
  1280. }
  1281. }
  1282. template <class T, class U, class V>
  1283. inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const mpl::int_<19>*) BOOST_MATH_NOEXCEPT(V)
  1284. {
  1285. if(x <= 1)
  1286. {
  1287. V x2 = x * x;
  1288. V t[4];
  1289. t[0] = a[18] * x2 + a[16];
  1290. t[1] = a[17] * x2 + a[15];
  1291. t[2] = b[18] * x2 + b[16];
  1292. t[3] = b[17] * x2 + b[15];
  1293. t[0] *= x2;
  1294. t[1] *= x2;
  1295. t[2] *= x2;
  1296. t[3] *= x2;
  1297. t[0] += static_cast<V>(a[14]);
  1298. t[1] += static_cast<V>(a[13]);
  1299. t[2] += static_cast<V>(b[14]);
  1300. t[3] += static_cast<V>(b[13]);
  1301. t[0] *= x2;
  1302. t[1] *= x2;
  1303. t[2] *= x2;
  1304. t[3] *= x2;
  1305. t[0] += static_cast<V>(a[12]);
  1306. t[1] += static_cast<V>(a[11]);
  1307. t[2] += static_cast<V>(b[12]);
  1308. t[3] += static_cast<V>(b[11]);
  1309. t[0] *= x2;
  1310. t[1] *= x2;
  1311. t[2] *= x2;
  1312. t[3] *= x2;
  1313. t[0] += static_cast<V>(a[10]);
  1314. t[1] += static_cast<V>(a[9]);
  1315. t[2] += static_cast<V>(b[10]);
  1316. t[3] += static_cast<V>(b[9]);
  1317. t[0] *= x2;
  1318. t[1] *= x2;
  1319. t[2] *= x2;
  1320. t[3] *= x2;
  1321. t[0] += static_cast<V>(a[8]);
  1322. t[1] += static_cast<V>(a[7]);
  1323. t[2] += static_cast<V>(b[8]);
  1324. t[3] += static_cast<V>(b[7]);
  1325. t[0] *= x2;
  1326. t[1] *= x2;
  1327. t[2] *= x2;
  1328. t[3] *= x2;
  1329. t[0] += static_cast<V>(a[6]);
  1330. t[1] += static_cast<V>(a[5]);
  1331. t[2] += static_cast<V>(b[6]);
  1332. t[3] += static_cast<V>(b[5]);
  1333. t[0] *= x2;
  1334. t[1] *= x2;
  1335. t[2] *= x2;
  1336. t[3] *= x2;
  1337. t[0] += static_cast<V>(a[4]);
  1338. t[1] += static_cast<V>(a[3]);
  1339. t[2] += static_cast<V>(b[4]);
  1340. t[3] += static_cast<V>(b[3]);
  1341. t[0] *= x2;
  1342. t[1] *= x2;
  1343. t[2] *= x2;
  1344. t[3] *= x2;
  1345. t[0] += static_cast<V>(a[2]);
  1346. t[1] += static_cast<V>(a[1]);
  1347. t[2] += static_cast<V>(b[2]);
  1348. t[3] += static_cast<V>(b[1]);
  1349. t[0] *= x2;
  1350. t[2] *= x2;
  1351. t[0] += static_cast<V>(a[0]);
  1352. t[2] += static_cast<V>(b[0]);
  1353. t[1] *= x;
  1354. t[3] *= x;
  1355. return (t[0] + t[1]) / (t[2] + t[3]);
  1356. }
  1357. else
  1358. {
  1359. V z = 1 / x;
  1360. V z2 = 1 / (x * x);
  1361. V t[4];
  1362. t[0] = a[0] * z2 + a[2];
  1363. t[1] = a[1] * z2 + a[3];
  1364. t[2] = b[0] * z2 + b[2];
  1365. t[3] = b[1] * z2 + b[3];
  1366. t[0] *= z2;
  1367. t[1] *= z2;
  1368. t[2] *= z2;
  1369. t[3] *= z2;
  1370. t[0] += static_cast<V>(a[4]);
  1371. t[1] += static_cast<V>(a[5]);
  1372. t[2] += static_cast<V>(b[4]);
  1373. t[3] += static_cast<V>(b[5]);
  1374. t[0] *= z2;
  1375. t[1] *= z2;
  1376. t[2] *= z2;
  1377. t[3] *= z2;
  1378. t[0] += static_cast<V>(a[6]);
  1379. t[1] += static_cast<V>(a[7]);
  1380. t[2] += static_cast<V>(b[6]);
  1381. t[3] += static_cast<V>(b[7]);
  1382. t[0] *= z2;
  1383. t[1] *= z2;
  1384. t[2] *= z2;
  1385. t[3] *= z2;
  1386. t[0] += static_cast<V>(a[8]);
  1387. t[1] += static_cast<V>(a[9]);
  1388. t[2] += static_cast<V>(b[8]);
  1389. t[3] += static_cast<V>(b[9]);
  1390. t[0] *= z2;
  1391. t[1] *= z2;
  1392. t[2] *= z2;
  1393. t[3] *= z2;
  1394. t[0] += static_cast<V>(a[10]);
  1395. t[1] += static_cast<V>(a[11]);
  1396. t[2] += static_cast<V>(b[10]);
  1397. t[3] += static_cast<V>(b[11]);
  1398. t[0] *= z2;
  1399. t[1] *= z2;
  1400. t[2] *= z2;
  1401. t[3] *= z2;
  1402. t[0] += static_cast<V>(a[12]);
  1403. t[1] += static_cast<V>(a[13]);
  1404. t[2] += static_cast<V>(b[12]);
  1405. t[3] += static_cast<V>(b[13]);
  1406. t[0] *= z2;
  1407. t[1] *= z2;
  1408. t[2] *= z2;
  1409. t[3] *= z2;
  1410. t[0] += static_cast<V>(a[14]);
  1411. t[1] += static_cast<V>(a[15]);
  1412. t[2] += static_cast<V>(b[14]);
  1413. t[3] += static_cast<V>(b[15]);
  1414. t[0] *= z2;
  1415. t[1] *= z2;
  1416. t[2] *= z2;
  1417. t[3] *= z2;
  1418. t[0] += static_cast<V>(a[16]);
  1419. t[1] += static_cast<V>(a[17]);
  1420. t[2] += static_cast<V>(b[16]);
  1421. t[3] += static_cast<V>(b[17]);
  1422. t[0] *= z2;
  1423. t[2] *= z2;
  1424. t[0] += static_cast<V>(a[18]);
  1425. t[2] += static_cast<V>(b[18]);
  1426. t[1] *= z;
  1427. t[3] *= z;
  1428. return (t[0] + t[1]) / (t[2] + t[3]);
  1429. }
  1430. }
  1431. }}}} // namespaces
  1432. #endif // include guard