1 // Use, modification and distribution are subject to the
2 // Boost Software License, Version 1.0.
3 // (See accompanying file LICENSE_1_0.txt
4 // or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 // Copyright Jeremy W. Murphy 2015.
7
8 // This file is written to be included from a Quickbook .qbk document.
9 // It can be compiled by the C++ compiler, and run. Any output can
10 // also be added here as comment or included or pasted in elsewhere.
11 // Caution: this file contains Quickbook markup as well as code
12 // and comments: don't change any of the special comment markups!
13
14 //[polynomial_arithmetic_0
15 /*`First include the essential polynomial header (and others) to make the example:
16 */
17 #include <boost/math/tools/polynomial.hpp>
18 //] [polynomial_arithmetic_0
19
20 #include <boost/array.hpp>
21 #include <boost/lexical_cast.hpp>
22 #include <boost/assert.hpp>
23
24 #include <iostream>
25 #include <stdexcept>
26 #include <cmath>
27 #include <string>
28 #include <utility>
29
30 //[polynomial_arithmetic_1
31 /*`and some using statements are convenient:
32 */
33
34 using std::string;
35 using std::exception;
36 using std::cout;
37 using std::abs;
38 using std::pair;
39
40 using namespace boost::math;
41 using namespace boost::math::tools; // for polynomial
42 using boost::lexical_cast;
43
44 //] [/polynomial_arithmetic_1]
45
46 template <typename T>
sign_str(T const & x)47 string sign_str(T const &x)
48 {
49 return x < 0 ? "-" : "+";
50 }
51
52 template <typename T>
inner_coefficient(T const & x)53 string inner_coefficient(T const &x)
54 {
55 string result(" " + sign_str(x) + " ");
56 if (abs(x) != T(1))
57 result += lexical_cast<string>(abs(x));
58 return result;
59 }
60
61 /*! Output in formula format.
62 For example: from a polynomial in Boost container storage [ 10, -6, -4, 3 ]
63 show as human-friendly formula notation: 3x^3 - 4x^2 - 6x + 10.
64 */
65 template <typename T>
formula_format(polynomial<T> const & a)66 string formula_format(polynomial<T> const &a)
67 {
68 string result;
69 if (a.size() == 0)
70 result += lexical_cast<string>(T(0));
71 else
72 {
73 // First one is a special case as it may need unary negate.
74 unsigned i = a.size() - 1;
75 if (a[i] < 0)
76 result += "-";
77 if (abs(a[i]) != T(1))
78 result += lexical_cast<string>(abs(a[i]));
79
80 if (i > 0)
81 {
82 result += "x";
83 if (i > 1)
84 {
85 result += "^" + lexical_cast<string>(i);
86 i--;
87 for (; i != 1; i--)
88 if (a[i])
89 result += inner_coefficient(a[i]) + "x^" + lexical_cast<string>(i);
90
91 if (a[i])
92 result += inner_coefficient(a[i]) + "x";
93 }
94 i--;
95
96 if (a[i])
97 result += " " + sign_str(a[i]) + " " + lexical_cast<string>(abs(a[i]));
98 }
99 }
100 return result;
101 } // string formula_format(polynomial<T> const &a)
102
103
main()104 int main()
105 {
106 cout << "Example: Polynomial arithmetic.\n\n";
107
108 try
109 {
110 //[polynomial_arithmetic_2
111 /*`Store the coefficients in a convenient way to access them,
112 then create some polynomials using construction from an iterator range,
113 and finally output in a 'pretty' formula format.
114
115 [tip Although we might conventionally write a polynomial from left to right
116 in descending order of degree, Boost.Math stores in [*ascending order of degree].]
117
118 Read/write for humans: 3x^3 - 4x^2 - 6x + 10
119 Boost polynomial storage: [ 10, -6, -4, 3 ]
120 */
121 boost::array<double, 4> const d3a = {{10, -6, -4, 3}};
122 polynomial<double> const a(d3a.begin(), d3a.end());
123
124 // With C++11 and later, you can also use initializer_list construction.
125 polynomial<double> const b{{-2.0, 1.0}};
126
127 // formula_format() converts from Boost storage to human notation.
128 cout << "a = " << formula_format(a)
129 << "\nb = " << formula_format(b) << "\n\n";
130
131 //] [/polynomial_arithmetic_2]
132
133 //[polynomial_arithmetic_3
134 // Now we can do arithmetic with the usual infix operators: + - * / and %.
135 polynomial<double> s = a + b;
136 cout << "a + b = " << formula_format(s) << "\n";
137 polynomial<double> d = a - b;
138 cout << "a - b = " << formula_format(d) << "\n";
139 polynomial<double> p = a * b;
140 cout << "a * b = " << formula_format(p) << "\n";
141 polynomial<double> q = a / b;
142 cout << "a / b = " << formula_format(q) << "\n";
143 polynomial<double> r = a % b;
144 cout << "a % b = " << formula_format(r) << "\n";
145 //] [/polynomial_arithmetic_3]
146
147 //[polynomial_arithmetic_4
148 /*`
149 Division is a special case where you can calculate two for the price of one.
150
151 Actually, quotient and remainder are always calculated together due to the nature
152 of the algorithm: the infix operators return one result and throw the other
153 away.
154
155 If you are doing a lot of division and want both the quotient and remainder, then
156 you don't want to do twice the work necessary.
157
158 In that case you can call the underlying function, [^quotient_remainder],
159 to get both results together as a pair.
160 */
161 pair< polynomial<double>, polynomial<double> > result;
162 result = quotient_remainder(a, b);
163 // Reassure ourselves that the result is the same.
164 BOOST_ASSERT(result.first == q);
165 BOOST_ASSERT(result.second == r);
166 //] [/polynomial_arithmetic_4]
167 //[polynomial_arithmetic_5
168 /*
169 We can use the right and left shift operators to add and remove a factor of x.
170 This has the same semantics as left and right shift for integers where it is a
171 factor of 2. x is the smallest prime factor of a polynomial as is 2 for integers.
172 */
173 cout << "Right and left shift operators.\n";
174 cout << "\n" << formula_format(p) << "\n";
175 cout << "... right shift by 1 ...\n";
176 p >>= 1;
177 cout << formula_format(p) << "\n";
178 cout << "... left shift by 2 ...\n";
179 p <<= 2;
180 cout << formula_format(p) << "\n";
181
182 /*
183 We can also give a meaning to odd and even for a polynomial that is consistent
184 with these operations: a polynomial is odd if it has a non-zero constant value,
185 even otherwise. That is:
186 x^2 + 1 odd
187 x^2 even
188 */
189 cout << std::boolalpha;
190 cout << "\nPrint whether a polynomial is odd.\n";
191 cout << formula_format(s) << " odd? " << odd(s) << "\n";
192 // We cheekily use the internal details to subtract the constant, making it even.
193 s -= s.data().front();
194 cout << formula_format(s) << " odd? " << odd(s) << "\n";
195 // And of course you can check if it is even:
196 cout << formula_format(s) << " even? " << even(s) << "\n";
197
198
199 //] [/polynomial_arithmetic_5]
200 //[polynomial_arithmetic_6]
201 /* For performance and convenience, we can test whether a polynomial is zero
202 * by implicitly converting to bool with the same semantics as int. */
203 polynomial<double> zero; // Default construction is 0.
204 cout << "zero: " << (zero ? "not zero" : "zero") << "\n";
205 cout << "r: " << (r ? "not zero" : "zero") << "\n";
206 /* We can also set a polynomial to zero without needing a another zero
207 * polynomial to assign to it. */
208 r.set_zero();
209 cout << "r: " << (r ? "not zero" : "zero") << "\n";
210 //] [/polynomial_arithmetic_6]
211 }
212 catch (exception const &e)
213 {
214 cout << "\nMessage from thrown exception was:\n " << e.what() << "\n";
215 }
216 } // int main()
217
218 /*
219 //[polynomial_output_1
220
221 a = 3x^3 - 4x^2 - 6x + 10
222 b = x - 2
223
224 //] [/polynomial_output_1]
225
226
227 //[polynomial_output_2
228
229 a + b = 3x^3 - 4x^2 - 5x + 8
230 a - b = 3x^3 - 4x^2 - 7x + 12
231 a * b = 3x^4 - 10x^3 + 2x^2 + 22x - 20
232 a / b = 3x^2 + 2x - 2
233 a % b = 6
234
235 //] [/polynomial_output_2]
236
237 */
238