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 #include <iostream>
7 #include <iomanip>
8 #include <fstream>
9 #include <string>
10 #include <sstream>
11
12 int max_order = 20;
13 const char* path_prefix = "..\\..\\..\\boost\\math\\tools\\detail\\polynomial_";
14 const char* path_prefix2 = "..\\..\\..\\boost\\math\\tools\\detail\\rational_";
15
16 const char* copyright_string =
17 "// (C) Copyright John Maddock 2007.\n"
18 "// Use, modification and distribution are subject to the\n"
19 "// Boost Software License, Version 1.0. (See accompanying file\n"
20 "// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)\n"
21 "//\n"
22 "// This file is machine generated, do not edit by hand\n\n";
23
24
print_polynomials(int max_order)25 void print_polynomials(int max_order)
26 {
27 for(int i = 2; i <= max_order; ++i)
28 {
29 std::stringstream filename;
30 filename << path_prefix << "horner1_" << i << ".hpp";
31 std::ofstream ofs(filename.str().c_str());
32 if(ofs.bad())
33 break;
34 //
35 // Output the boilerplate at the top of the header:
36 //
37 ofs << copyright_string <<
38 "// Polynomial evaluation using Horners rule\n"
39 "#ifndef BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n"
40 "#define BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n\n"
41 "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
42 "template <class T, class V>\n"
43 "inline V evaluate_polynomial_c_imp(const T*, const V&, const boost::integral_constant<int, 0>*)\n"
44 "{\n"
45 " return static_cast<V>(0);\n"
46 "}\n"
47 "\n"
48 "template <class T, class V>\n"
49 "inline V evaluate_polynomial_c_imp(const T* a, const V&, const boost::integral_constant<int, 1>*)\n"
50 "{\n"
51 " return static_cast<V>(a[0]);\n"
52 "}\n\n";
53
54 for(int order = 2; order <= i; ++order)
55 {
56 ofs <<
57 "template <class T, class V>\n"
58 "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const boost::integral_constant<int, " << order << ">*)\n"
59 "{\n"
60 " return static_cast<V>(";
61
62 for(int bracket = 2; bracket < order; ++bracket)
63 ofs << "(";
64 ofs << "a[" << order - 1 << "] * x + a[" << order - 2 << "]" ;
65 for(int item = order - 3; item >= 0; --item)
66 {
67 ofs << ") * x + a[" << item << "]";
68 }
69
70 ofs << ");\n"
71 "}\n\n";
72 }
73 //
74 // And finally the boilerplate at the end of the header:
75 //
76 ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";
77
78 filename.str("");
79 filename << path_prefix << "horner2_" << i << ".hpp";
80 ofs.close();
81 ofs.open(filename.str().c_str());
82 if(ofs.bad())
83 break;
84 //
85 // Output the boilerplate at the top of the header:
86 //
87 ofs << copyright_string <<
88 "// Polynomial evaluation using second order Horners rule\n"
89 "#ifndef BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n"
90 "#define BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n\n"
91 "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
92 "template <class T, class V>\n"
93 "inline V evaluate_polynomial_c_imp(const T*, const V&, const boost::integral_constant<int, 0>*)\n"
94 "{\n"
95 " return static_cast<V>(0);\n"
96 "}\n"
97 "\n"
98 "template <class T, class V>\n"
99 "inline V evaluate_polynomial_c_imp(const T* a, const V&, const boost::integral_constant<int, 1>*)\n"
100 "{\n"
101 " return static_cast<V>(a[0]);\n"
102 "}\n\n"
103 "template <class T, class V>\n"
104 "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const boost::integral_constant<int, 2>*)\n"
105 "{\n"
106 " return static_cast<V>(a[1] * x + a[0]);\n"
107 "}\n\n"
108 "template <class T, class V>\n"
109 "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const boost::integral_constant<int, 3>*)\n"
110 "{\n"
111 " return static_cast<V>((a[2] * x + a[1]) * x + a[0]);\n"
112 "}\n\n"
113 "template <class T, class V>\n"
114 "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const boost::integral_constant<int, 4>*)\n"
115 "{\n"
116 " return static_cast<V>(((a[3] * x + a[2]) * x + a[1]) * x + a[0]);\n"
117 "}\n\n";
118
119 for(int order = 5; order <= i; ++order)
120 {
121 ofs <<
122 "template <class T, class V>\n"
123 "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const boost::integral_constant<int, " << order << ">*)\n"
124 "{\n"
125 " V x2 = x * x;\n"
126 " return static_cast<V>(";
127
128 if(order & 1)
129 {
130 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
131 ofs << "(";
132 ofs << "a[" << order - 1 << "] * x2 + a[" << order - 3 << "]" ;
133 for(int item = order - 5; item >= 0; item -= 2)
134 {
135 ofs << ") * x2 + a[" << item << "]";
136 }
137 ofs << " + ";
138 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
139 ofs << "(";
140 ofs << "a[" << order - 2 << "] * x2 + a[" << order - 4 << "]" ;
141 for(int item = order - 6; item >= 0; item -= 2)
142 {
143 ofs << ") * x2 + a[" << item << "]";
144 }
145 ofs << ") * x";
146 }
147 else
148 {
149 for(int bracket = 0; bracket < (order - 1) / 2; ++bracket)
150 ofs << "(";
151 ofs << "a[" << order - 1 << "] * x2 + a[" << order - 3 << "]" ;
152 for(int item = order - 5; item >= 0; item -= 2)
153 {
154 ofs << ") * x2 + a[" << item << "]";
155 }
156 ofs << ") * x + ";
157 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
158 ofs << "(";
159 ofs << "a[" << order - 2 << "] * x2 + a[" << order - 4 << "]" ;
160 for(int item = order - 6; item >= 0; item -= 2)
161 {
162 ofs << ") * x2 + a[" << item << "]";
163 }
164 }
165 ofs << ");\n"
166 "}\n\n";
167 }
168 //
169 // And finally the boilerplate at the end of the header:
170 //
171 ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";
172
173
174 filename.str("");
175 filename << path_prefix << "horner3_" << i << ".hpp";
176 ofs.close();
177 ofs.open(filename.str().c_str());
178 if(ofs.bad())
179 break;
180 //
181 // Output the boilerplate at the top of the header:
182 //
183 ofs << copyright_string <<
184 "// Unrolled polynomial evaluation using second order Horners rule\n"
185 "#ifndef BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n"
186 "#define BOOST_MATH_TOOLS_POLY_EVAL_" << i << "_HPP\n\n"
187 "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
188 "template <class T, class V>\n"
189 "inline V evaluate_polynomial_c_imp(const T*, const V&, const boost::integral_constant<int, 0>*)\n"
190 "{\n"
191 " return static_cast<V>(0);\n"
192 "}\n"
193 "\n"
194 "template <class T, class V>\n"
195 "inline V evaluate_polynomial_c_imp(const T* a, const V&, const boost::integral_constant<int, 1>*)\n"
196 "{\n"
197 " return static_cast<V>(a[0]);\n"
198 "}\n\n"
199 "template <class T, class V>\n"
200 "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const boost::integral_constant<int, 2>*)\n"
201 "{\n"
202 " return static_cast<V>(a[1] * x + a[0]);\n"
203 "}\n\n"
204 "template <class T, class V>\n"
205 "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const boost::integral_constant<int, 3>*)\n"
206 "{\n"
207 " return static_cast<V>((a[2] * x + a[1]) * x + a[0]);\n"
208 "}\n\n"
209 "template <class T, class V>\n"
210 "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const boost::integral_constant<int, 4>*)\n"
211 "{\n"
212 " return static_cast<V>(((a[3] * x + a[2]) * x + a[1]) * x + a[0]);\n"
213 "}\n\n";
214
215 for(int order = 5; order <= i; ++order)
216 {
217 ofs <<
218 "template <class T, class V>\n"
219 "inline V evaluate_polynomial_c_imp(const T* a, const V& x, const boost::integral_constant<int, " << order << ">*)\n"
220 "{\n"
221 " V x2 = x * x;\n"
222 " V t[2];\n";
223
224 if(order & 1)
225 {
226 ofs << " t[0] = static_cast<V>(a[" << order - 1 << "] * x2 + a[" << order - 3 << "]);\n" ;
227 ofs << " t[1] = static_cast<V>(a[" << order - 2 << "] * x2 + a[" << order - 4 << "]);\n" ;
228 for(int item = order - 5; item >= 0; item -= 2)
229 {
230 ofs << " t[0] *= x2;\n";
231 if(item - 1 >= 0)
232 ofs << " t[1] *= x2;\n";
233 ofs << " t[0] += static_cast<V>(a[" << item << "]);\n";
234 if(item - 1 >= 0)
235 ofs << " t[1] += static_cast<V>(a[" << item - 1 << "]);\n";
236 }
237 ofs <<
238 " t[1] *= x;\n"
239 " return t[0] + t[1];\n";
240 }
241 else
242 {
243 ofs << " t[0] = a[" << order - 1 << "] * x2 + a[" << order - 3 << "];\n" ;
244 ofs << " t[1] = a[" << order - 2 << "] * x2 + a[" << order - 4 << "];\n" ;
245 for(int item = order - 5; item >= 0; item -= 2)
246 {
247 ofs << " t[0] *= x2;\n";
248 if(item - 1 >= 0)
249 ofs << " t[1] *= x2;\n";
250 ofs << " t[0] += static_cast<V>(a[" << item << "]);\n";
251 if(item - 1 >= 0)
252 ofs << " t[1] += static_cast<V>(a[" << item - 1 << "]);\n";
253 }
254 ofs << " t[0] *= x;\n";
255 ofs << " return t[0] + t[1];\n";
256 }
257 ofs << "}\n\n";
258 }
259 //
260 // And finally the boilerplate at the end of the header:
261 //
262 ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";
263 }
264 }
265
print_rationals(int max_order)266 void print_rationals(int max_order)
267 {
268 for(int i = 2; i <= max_order; ++i)
269 {
270 std::stringstream filename;
271 filename << path_prefix2 << "horner1_" << i << ".hpp";
272 std::ofstream ofs(filename.str().c_str());
273 if(ofs.bad())
274 break;
275 //
276 // Output the boilerplate at the top of the header:
277 //
278 ofs << copyright_string <<
279 "// Polynomial evaluation using Horners rule\n"
280 "#ifndef BOOST_MATH_TOOLS_POLY_RAT_" << i << "_HPP\n"
281 "#define BOOST_MATH_TOOLS_POLY_RAT_" << i << "_HPP\n\n"
282 "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
283 "template <class T, class U, class V>\n"
284 "inline V evaluate_rational_c_imp(const T*, const U*, const V&, const boost::integral_constant<int, 0>*)\n"
285 "{\n"
286 " return static_cast<V>(0);\n"
287 "}\n"
288 "\n"
289 "template <class T, class U, class V>\n"
290 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const boost::integral_constant<int, 1>*)\n"
291 "{\n"
292 " return static_cast<V>(a[0]) / static_cast<V>(b[0]);\n"
293 "}\n\n";
294
295 for(int order = 2; order <= i; ++order)
296 {
297 ofs <<
298 "template <class T, class U, class V>\n"
299 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const boost::integral_constant<int, " << order << ">*)\n"
300 "{\n"
301 " if(x <= 1)\n"
302 " return static_cast<V>((";
303
304 for(int bracket = 2; bracket < order; ++bracket)
305 ofs << "(";
306 ofs << "a[" << order - 1 << "] * x + a[" << order - 2 << "]" ;
307 for(int item = order - 3; item >= 0; --item)
308 {
309 ofs << ") * x + a[" << item << "]";
310 }
311
312 ofs << ") / (";
313 for(int bracket = 2; bracket < order; ++bracket)
314 ofs << "(";
315 ofs << "b[" << order - 1 << "] * x + b[" << order - 2 << "]" ;
316 for(int item = order - 3; item >= 0; --item)
317 {
318 ofs << ") * x + b[" << item << "]";
319 }
320
321 ofs << "));\n else\n {\n V z = 1 / x;\n return static_cast<V>((";
322
323 for(int bracket = order - 1; bracket > 1; --bracket)
324 ofs << "(";
325 ofs << "a[" << 0 << "] * z + a[" << 1 << "]" ;
326 for(int item = 2; item <= order - 1; ++item)
327 {
328 ofs << ") * z + a[" << item << "]";
329 }
330
331 ofs << ") / (";
332 for(int bracket = 2; bracket < order; ++bracket)
333 ofs << "(";
334 ofs << "b[" << 0 << "] * z + b[" << 1 << "]" ;
335 for(int item = 2; item <= order - 1; ++item)
336 {
337 ofs << ") * z + b[" << item << "]";
338 }
339
340 ofs << "));\n }\n";
341
342 ofs << "}\n\n";
343 }
344 //
345 // And finally the boilerplate at the end of the header:
346 //
347 ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";
348
349 filename.str("");
350 filename << path_prefix2 << "horner2_" << i << ".hpp";
351 ofs.close();
352 ofs.open(filename.str().c_str());
353 if(ofs.bad())
354 break;
355 //
356 // Output the boilerplate at the top of the header:
357 //
358 ofs << copyright_string <<
359 "// Polynomial evaluation using second order Horners rule\n"
360 "#ifndef BOOST_MATH_TOOLS_RAT_EVAL_" << i << "_HPP\n"
361 "#define BOOST_MATH_TOOLS_RAT_EVAL_" << i << "_HPP\n\n"
362 "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
363 "template <class T, class U, class V>\n"
364 "inline V evaluate_rational_c_imp(const T*, const U*, const V&, const boost::integral_constant<int, 0>*)\n"
365 "{\n"
366 " return static_cast<V>(0);\n"
367 "}\n"
368 "\n"
369 "template <class T, class U, class V>\n"
370 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const boost::integral_constant<int, 1>*)\n"
371 "{\n"
372 " return static_cast<V>(a[0]) / static_cast<V>(b[0]);\n"
373 "}\n\n"
374 "template <class T, class U, class V>\n"
375 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const boost::integral_constant<int, 2>*)\n"
376 "{\n"
377 " return static_cast<V>((a[1] * x + a[0]) / (b[1] * x + b[0]));\n"
378 "}\n\n"
379 "template <class T, class U, class V>\n"
380 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const boost::integral_constant<int, 3>*)\n"
381 "{\n"
382 " return static_cast<V>(((a[2] * x + a[1]) * x + a[0]) / ((b[2] * x + b[1]) * x + b[0]));\n"
383 "}\n\n"
384 "template <class T, class U, class V>\n"
385 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const boost::integral_constant<int, 4>*)\n"
386 "{\n"
387 " 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"
388 "}\n\n";
389
390 for(int order = 5; order <= i; ++order)
391 {
392 ofs <<
393 "template <class T, class U, class V>\n"
394 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const boost::integral_constant<int, " << order << ">*)\n"
395 "{\n"
396 " if(x <= 1)\n {\n"
397 " V x2 = x * x;\n"
398 " return static_cast<V>((";
399
400 if(order & 1)
401 {
402 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
403 ofs << "(";
404 ofs << "a[" << order - 1 << "] * x2 + a[" << order - 3 << "]" ;
405 for(int item = order - 5; item >= 0; item -= 2)
406 {
407 ofs << ") * x2 + a[" << item << "]";
408 }
409 ofs << " + ";
410 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
411 ofs << "(";
412 ofs << "a[" << order - 2 << "] * x2 + a[" << order - 4 << "]" ;
413 for(int item = order - 6; item >= 0; item -= 2)
414 {
415 ofs << ") * x2 + a[" << item << "]";
416 }
417 ofs << ") * x";
418 }
419 else
420 {
421 for(int bracket = 0; bracket < (order - 1) / 2; ++bracket)
422 ofs << "(";
423 ofs << "a[" << order - 1 << "] * x2 + a[" << order - 3 << "]" ;
424 for(int item = order - 5; item >= 0; item -= 2)
425 {
426 ofs << ") * x2 + a[" << item << "]";
427 }
428 ofs << ") * x + ";
429 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
430 ofs << "(";
431 ofs << "a[" << order - 2 << "] * x2 + a[" << order - 4 << "]" ;
432 for(int item = order - 6; item >= 0; item -= 2)
433 {
434 ofs << ") * x2 + a[" << item << "]";
435 }
436 }
437 ofs << ") / (";
438 if(order & 1)
439 {
440 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
441 ofs << "(";
442 ofs << "b[" << order - 1 << "] * x2 + b[" << order - 3 << "]" ;
443 for(int item = order - 5; item >= 0; item -= 2)
444 {
445 ofs << ") * x2 + b[" << item << "]";
446 }
447 ofs << " + ";
448 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
449 ofs << "(";
450 ofs << "b[" << order - 2 << "] * x2 + b[" << order - 4 << "]" ;
451 for(int item = order - 6; item >= 0; item -= 2)
452 {
453 ofs << ") * x2 + b[" << item << "]";
454 }
455 ofs << ") * x";
456 }
457 else
458 {
459 for(int bracket = 0; bracket < (order - 1) / 2; ++bracket)
460 ofs << "(";
461 ofs << "b[" << order - 1 << "] * x2 + b[" << order - 3 << "]" ;
462 for(int item = order - 5; item >= 0; item -= 2)
463 {
464 ofs << ") * x2 + b[" << item << "]";
465 }
466 ofs << ") * x + ";
467 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
468 ofs << "(";
469 ofs << "b[" << order - 2 << "] * x2 + b[" << order - 4 << "]" ;
470 for(int item = order - 6; item >= 0; item -= 2)
471 {
472 ofs << ") * x2 + b[" << item << "]";
473 }
474 }
475
476 ofs << "));\n }\n else\n {\n V z = 1 / x;\n V z2 = 1 / (x * x);\n return static_cast<V>((";
477
478 if(order & 1)
479 {
480 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
481 ofs << "(";
482 ofs << "a[" << 0 << "] * z2 + a[" << 2 << "]" ;
483 for(int item = 4; item < order; item += 2)
484 {
485 ofs << ") * z2 + a[" << item << "]";
486 }
487 ofs << " + ";
488 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
489 ofs << "(";
490 ofs << "a[" << 1 << "] * z2 + a[" << 3 << "]" ;
491 for(int item = 5; item < order; item += 2)
492 {
493 ofs << ") * z2 + a[" << item << "]";
494 }
495 ofs << ") * z";
496 }
497 else
498 {
499 for(int bracket = 0; bracket < (order - 1) / 2; ++bracket)
500 ofs << "(";
501 ofs << "a[" << 0 << "] * z2 + a[" << 2 << "]" ;
502 for(int item = 4; item < order; item += 2)
503 {
504 ofs << ") * z2 + a[" << item << "]";
505 }
506 ofs << ") * z + ";
507 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
508 ofs << "(";
509 ofs << "a[" << 1 << "] * z2 + a[" << 3 << "]" ;
510 for(int item = 5; item < order; item += 2)
511 {
512 ofs << ") * z2 + a[" << item << "]";
513 }
514 }
515
516 ofs << ") / (";
517
518 if(order & 1)
519 {
520 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
521 ofs << "(";
522 ofs << "b[" << 0 << "] * z2 + b[" << 2 << "]" ;
523 for(int item = 4; item < order; item += 2)
524 {
525 ofs << ") * z2 + b[" << item << "]";
526 }
527 ofs << " + ";
528 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
529 ofs << "(";
530 ofs << "b[" << 1 << "] * z2 + b[" << 3 << "]" ;
531 for(int item = 5; item < order; item += 2)
532 {
533 ofs << ") * z2 + b[" << item << "]";
534 }
535 ofs << ") * z";
536 }
537 else
538 {
539 for(int bracket = 0; bracket < (order - 1) / 2; ++bracket)
540 ofs << "(";
541 ofs << "b[" << 0 << "] * z2 + b[" << 2 << "]" ;
542 for(int item = 4; item < order; item += 2)
543 {
544 ofs << ") * z2 + b[" << item << "]";
545 }
546 ofs << ") * z + ";
547 for(int bracket = 0; bracket < (order - 1) / 2 - 1; ++bracket)
548 ofs << "(";
549 ofs << "b[" << 1 << "] * z2 + b[" << 3 << "]" ;
550 for(int item = 5; item < order; item += 2)
551 {
552 ofs << ") * z2 + b[" << item << "]";
553 }
554 }
555 ofs << "));\n }\n";
556
557 ofs << "}\n\n";
558 }
559 //
560 // And finally the boilerplate at the end of the header:
561 //
562 ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";
563
564
565 filename.str("");
566 filename << path_prefix2 << "horner3_" << i << ".hpp";
567 ofs.close();
568 ofs.open(filename.str().c_str());
569 if(ofs.bad())
570 break;
571 //
572 // Output the boilerplate at the top of the header:
573 //
574 ofs << copyright_string <<
575 "// Polynomial evaluation using second order Horners rule\n"
576 "#ifndef BOOST_MATH_TOOLS_RAT_EVAL_" << i << "_HPP\n"
577 "#define BOOST_MATH_TOOLS_RAT_EVAL_" << i << "_HPP\n\n"
578 "namespace boost{ namespace math{ namespace tools{ namespace detail{\n\n"
579 "template <class T, class U, class V>\n"
580 "inline V evaluate_rational_c_imp(const T*, const U*, const V&, const boost::integral_constant<int, 0>*)\n"
581 "{\n"
582 " return static_cast<V>(0);\n"
583 "}\n"
584 "\n"
585 "template <class T, class U, class V>\n"
586 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V&, const boost::integral_constant<int, 1>*)\n"
587 "{\n"
588 " return static_cast<V>(a[0]) / static_cast<V>(b[0]);\n"
589 "}\n\n"
590 "template <class T, class U, class V>\n"
591 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const boost::integral_constant<int, 2>*)\n"
592 "{\n"
593 " return static_cast<V>((a[1] * x + a[0]) / (b[1] * x + b[0]));\n"
594 "}\n\n"
595 "template <class T, class U, class V>\n"
596 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const boost::integral_constant<int, 3>*)\n"
597 "{\n"
598 " return static_cast<V>(((a[2] * x + a[1]) * x + a[0]) / ((b[2] * x + b[1]) * x + b[0]));\n"
599 "}\n\n"
600 "template <class T, class U, class V>\n"
601 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const boost::integral_constant<int, 4>*)\n"
602 "{\n"
603 " 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"
604 "}\n\n";
605
606 for(int order = 5; order <= i; ++order)
607 {
608 ofs <<
609 "template <class T, class U, class V>\n"
610 "inline V evaluate_rational_c_imp(const T* a, const U* b, const V& x, const boost::integral_constant<int, " << order << ">*)\n"
611 "{\n"
612 " if(x <= 1)\n {\n"
613 " V x2 = x * x;\n"
614 " V t[4];\n";
615
616 if(order & 1)
617 {
618 ofs << " t[0] = a[" << order - 1 << "] * x2 + a[" << order - 3 << "];\n" ;
619 ofs << " t[1] = a[" << order - 2 << "] * x2 + a[" << order - 4 << "];\n" ;
620 ofs << " t[2] = b[" << order - 1 << "] * x2 + b[" << order - 3 << "];\n" ;
621 ofs << " t[3] = b[" << order - 2 << "] * x2 + b[" << order - 4 << "];\n" ;
622 for(int item = order - 5; item >= 0; item -= 2)
623 {
624 ofs << " t[0] *= x2;\n";
625 if(item - 1 >= 0)
626 ofs << " t[1] *= x2;\n";
627 ofs << " t[2] *= x2;\n";
628 if(item - 1 >= 0)
629 ofs << " t[3] *= x2;\n";
630 ofs << " t[0] += static_cast<V>(a[" << item << "]);\n";
631 if(item - 1 >= 0)
632 ofs << " t[1] += static_cast<V>(a[" << item - 1 << "]);\n";
633 ofs << " t[2] += static_cast<V>(b[" << item << "]);\n";
634 if(item - 1 >= 0)
635 ofs << " t[3] += static_cast<V>(b[" << item - 1 << "]);\n";
636 }
637 ofs << " t[1] *= x;\n";
638 ofs << " t[3] *= x;\n";
639 }
640 else
641 {
642 ofs << " t[0] = a[" << order - 1 << "] * x2 + a[" << order - 3 << "];\n" ;
643 ofs << " t[1] = a[" << order - 2 << "] * x2 + a[" << order - 4 << "];\n" ;
644 ofs << " t[2] = b[" << order - 1 << "] * x2 + b[" << order - 3 << "];\n" ;
645 ofs << " t[3] = b[" << order - 2 << "] * x2 + b[" << order - 4 << "];\n" ;
646 for(int item = order - 5; item >= 0; item -= 2)
647 {
648 ofs << " t[0] *= x2;\n";
649 if(item - 1 >= 0)
650 ofs << " t[1] *= x2;\n";
651 ofs << " t[2] *= x2;\n";
652 if(item - 1 >= 0)
653 ofs << " t[3] *= x2;\n";
654 ofs << " t[0] += static_cast<V>(a[" << item << "]);\n";
655 if(item - 1 >= 0)
656 ofs << " t[1] += static_cast<V>(a[" << item - 1 << "]);\n";
657 ofs << " t[2] += static_cast<V>(b[" << item << "]);\n";
658 if(item - 1 >= 0)
659 ofs << " t[3] += static_cast<V>(b[" << item - 1 << "]);\n";
660 }
661 ofs << " t[0] *= x;\n";
662 ofs << " t[2] *= x;\n";
663 }
664 ofs << " return (t[0] + t[1]) / (t[2] + t[3]);\n";
665
666 ofs << " }\n else\n {\n V z = 1 / x;\n V z2 = 1 / (x * x);\n V t[4];\n";
667
668 if(order & 1)
669 {
670 ofs << " t[0] = a[" << 0 << "] * z2 + a[" << 2 << "];\n" ;
671 ofs << " t[1] = a[" << 1 << "] * z2 + a[" << 3 << "];\n" ;
672 ofs << " t[2] = b[" << 0 << "] * z2 + b[" << 2 << "];\n" ;
673 ofs << " t[3] = b[" << 1 << "] * z2 + b[" << 3 << "];\n" ;
674 for(int item = 4; item < order; item += 2)
675 {
676 ofs << " t[0] *= z2;\n";
677 if(item + 1 < order)
678 ofs << " t[1] *= z2;\n";
679 ofs << " t[2] *= z2;\n";
680 if(item + 1 < order)
681 ofs << " t[3] *= z2;\n";
682 ofs << " t[0] += static_cast<V>(a[" << item << "]);\n";
683 if(item + 1 < order)
684 ofs << " t[1] += static_cast<V>(a[" << item + 1 << "]);\n";
685 ofs << " t[2] += static_cast<V>(b[" << item << "]);\n";
686 if(item + 1 < order)
687 ofs << " t[3] += static_cast<V>(b[" << item + 1 << "]);\n";
688 }
689 ofs << " t[1] *= z;\n";
690 ofs << " t[3] *= z;\n";
691 }
692 else
693 {
694 ofs << " t[0] = a[" << 0 << "] * z2 + a[" << 2 << "];\n" ;
695 ofs << " t[1] = a[" << 1 << "] * z2 + a[" << 3 << "];\n" ;
696 ofs << " t[2] = b[" << 0 << "] * z2 + b[" << 2 << "];\n" ;
697 ofs << " t[3] = b[" << 1 << "] * z2 + b[" << 3 << "];\n" ;
698 for(int item = 4; item < order; item += 2)
699 {
700 ofs << " t[0] *= z2;\n";
701 if(item + 1 < order)
702 ofs << " t[1] *= z2;\n";
703 ofs << " t[2] *= z2;\n";
704 if(item + 1 < order)
705 ofs << " t[3] *= z2;\n";
706 ofs << " t[0] += static_cast<V>(a[" << item << "]);\n";
707 if(item + 1 < order)
708 ofs << " t[1] += static_cast<V>(a[" << item + 1 << "]);\n";
709 ofs << " t[2] += static_cast<V>(b[" << item << "]);\n";
710 if(item + 1 < order)
711 ofs << " t[3] += static_cast<V>(b[" << item + 1 << "]);\n";
712 }
713 ofs << " t[0] *= z;\n";
714 ofs << " t[2] *= z;\n";
715 }
716 ofs << " return (t[0] + t[1]) / (t[2] + t[3]);\n }\n";
717
718 ofs << "}\n\n";
719 }
720 //
721 // And finally the boilerplate at the end of the header:
722 //
723 ofs << "\n}}}} // namespaces\n\n#endif // include guard\n\n";
724 }
725 }
726
main()727 int main()
728 {
729 for(int i = 2; i <= max_order; ++i)
730 {
731 print_polynomials(i);
732 print_rationals(i);
733 }
734 return 0;
735 }
736
737
738
739