1 // (C) Copyright John Maddock 2019.
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 // Contains Quickbook snippets used by boost/libs/multiprecision/doc/multiprecision.qbk,
7 // section Literal Types and constexpr Support.
8
9 #include <iostream>
10 #include <boost/config.hpp>
11
12 #ifdef BOOST_HAS_FLOAT128
13 #include <boost/multiprecision/float128.hpp>
14 #endif
15
16 //[constexpr_circle
17 #include <boost/math/constants/constants.hpp> // For constant pi with full precision of type T.
18 // using boost::math::constants::pi;
19
20 template <class T>
circumference(T radius)21 inline constexpr T circumference(T radius)
22 {
23 return 2 * boost::math::constants::pi<T>() * radius;
24 }
25
26 template <class T>
area(T radius)27 inline constexpr T area(T radius)
28 {
29 return boost::math::constants::pi<T>() * radius * radius;
30 }
31 //] [/constexpr_circle]
32
33 template <class T, unsigned Order>
34 struct const_polynomial
35 {
36 public:
37 T data[Order + 1];
38
39 public:
const_polynomialconst_polynomial40 constexpr const_polynomial(T val = 0) : data{val} {}
const_polynomialconst_polynomial41 constexpr const_polynomial(const std::initializer_list<T>& init) : data{}
42 {
43 if (init.size() > Order + 1)
44 throw std::range_error("Too many initializers in list!");
45 for (unsigned i = 0; i < init.size(); ++i)
46 data[i] = init.begin()[i];
47 }
operator []const_polynomial48 constexpr T& operator[](std::size_t N)
49 {
50 return data[N];
51 }
operator []const_polynomial52 constexpr const T& operator[](std::size_t N) const
53 {
54 return data[N];
55 }
56 template <class U>
operator ()const_polynomial57 constexpr T operator()(U val)const
58 {
59 T result = data[Order];
60 for (unsigned i = Order; i > 0; --i)
61 {
62 result *= val;
63 result += data[i - 1];
64 }
65 return result;
66 }
derivativeconst_polynomial67 constexpr const_polynomial<T, Order - 1> derivative() const
68 {
69 const_polynomial<T, Order - 1> result;
70 for (unsigned i = 1; i <= Order; ++i)
71 {
72 result[i - 1] = (*this)[i] * i;
73 }
74 return result;
75 }
operator -const_polynomial76 constexpr const_polynomial operator-()
77 {
78 const_polynomial t(*this);
79 for (unsigned i = 0; i <= Order; ++i)
80 t[i] = -t[i];
81 return t;
82 }
83 template <class U>
operator *=const_polynomial84 constexpr const_polynomial& operator*=(U val)
85 {
86 for (unsigned i = 0; i <= Order; ++i)
87 data[i] = data[i] * val;
88 return *this;
89 }
90 template <class U>
operator /=const_polynomial91 constexpr const_polynomial& operator/=(U val)
92 {
93 for (unsigned i = 0; i <= Order; ++i)
94 data[i] = data[i] / val;
95 return *this;
96 }
97 template <class U>
operator +=const_polynomial98 constexpr const_polynomial& operator+=(U val)
99 {
100 data[0] += val;
101 return *this;
102 }
103 template <class U>
operator -=const_polynomial104 constexpr const_polynomial& operator-=(U val)
105 {
106 data[0] -= val;
107 return *this;
108 }
109 };
110
111 template <class T, unsigned Order1, unsigned Order2>
operator +(const const_polynomial<T,Order1> & a,const const_polynomial<T,Order2> & b)112 inline constexpr const_polynomial<T, (Order1 > Order2 ? Order1 : Order2)> operator+(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
113 {
114 if
115 constexpr(Order1 > Order2)
116 {
117 const_polynomial<T, Order1> result(a);
118 for (unsigned i = 0; i <= Order2; ++i)
119 result[i] += b[i];
120 return result;
121 }
122 else
123 {
124 const_polynomial<T, Order2> result(b);
125 for (unsigned i = 0; i <= Order1; ++i)
126 result[i] += a[i];
127 return result;
128 }
129 }
130 template <class T, unsigned Order1, unsigned Order2>
operator -(const const_polynomial<T,Order1> & a,const const_polynomial<T,Order2> & b)131 inline constexpr const_polynomial<T, (Order1 > Order2 ? Order1 : Order2)> operator-(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
132 {
133 if
134 constexpr(Order1 > Order2)
135 {
136 const_polynomial<T, Order1> result(a);
137 for (unsigned i = 0; i <= Order2; ++i)
138 result[i] -= b[i];
139 return result;
140 }
141 else
142 {
143 const_polynomial<T, Order2> result(b);
144 for (unsigned i = 0; i <= Order1; ++i)
145 result[i] = a[i] - b[i];
146 return result;
147 }
148 }
149 template <class T, unsigned Order1, unsigned Order2>
operator *(const const_polynomial<T,Order1> & a,const const_polynomial<T,Order2> & b)150 inline constexpr const_polynomial<T, Order1 + Order2> operator*(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
151 {
152 const_polynomial<T, Order1 + Order2> result;
153 for (unsigned i = 0; i <= Order1; ++i)
154 {
155 for (unsigned j = 0; j <= Order2; ++j)
156 {
157 result[i + j] += a[i] * b[j];
158 }
159 }
160 return result;
161 }
162 template <class T, unsigned Order, class U>
operator *(const const_polynomial<T,Order> & a,const U & b)163 inline constexpr const_polynomial<T, Order> operator*(const const_polynomial<T, Order>& a, const U& b)
164 {
165 const_polynomial<T, Order> result(a);
166 for (unsigned i = 0; i <= Order; ++i)
167 {
168 result[i] *= b;
169 }
170 return result;
171 }
172 template <class U, class T, unsigned Order>
operator *(const U & b,const const_polynomial<T,Order> & a)173 inline constexpr const_polynomial<T, Order> operator*(const U& b, const const_polynomial<T, Order>& a)
174 {
175 const_polynomial<T, Order> result(a);
176 for (unsigned i = 0; i <= Order; ++i)
177 {
178 result[i] *= b;
179 }
180 return result;
181 }
182 template <class T, unsigned Order, class U>
operator /(const const_polynomial<T,Order> & a,const U & b)183 inline constexpr const_polynomial<T, Order> operator/(const const_polynomial<T, Order>& a, const U& b)
184 {
185 const_polynomial<T, Order> result;
186 for (unsigned i = 0; i <= Order; ++i)
187 {
188 result[i] /= b;
189 }
190 return result;
191 }
192
193 //[hermite_example
194 template <class T, unsigned Order>
195 class hermite_polynomial
196 {
197 const_polynomial<T, Order> m_data;
198
199 public:
hermite_polynomial()200 constexpr hermite_polynomial() : m_data(hermite_polynomial<T, Order - 1>().data() * const_polynomial<T, 1>{0, 2} - hermite_polynomial<T, Order - 1>().data().derivative())
201 {
202 }
data() const203 constexpr const const_polynomial<T, Order>& data() const
204 {
205 return m_data;
206 }
operator [](std::size_t N) const207 constexpr const T& operator[](std::size_t N)const
208 {
209 return m_data[N];
210 }
211 template <class U>
operator ()(U val) const212 constexpr T operator()(U val)const
213 {
214 return m_data(val);
215 }
216 };
217 //] [/hermite_example]
218
219 //[hermite_example2
220 template <class T>
221 class hermite_polynomial<T, 0>
222 {
223 const_polynomial<T, 0> m_data;
224
225 public:
hermite_polynomial()226 constexpr hermite_polynomial() : m_data{1} {}
data() const227 constexpr const const_polynomial<T, 0>& data() const
228 {
229 return m_data;
230 }
operator [](std::size_t N) const231 constexpr const T& operator[](std::size_t N) const
232 {
233 return m_data[N];
234 }
235 template <class U>
operator ()(U val)236 constexpr T operator()(U val)
237 {
238 return m_data(val);
239 }
240 };
241
242 template <class T>
243 class hermite_polynomial<T, 1>
244 {
245 const_polynomial<T, 1> m_data;
246
247 public:
hermite_polynomial()248 constexpr hermite_polynomial() : m_data{0, 2} {}
data() const249 constexpr const const_polynomial<T, 1>& data() const
250 {
251 return m_data;
252 }
operator [](std::size_t N) const253 constexpr const T& operator[](std::size_t N) const
254 {
255 return m_data[N];
256 }
257 template <class U>
operator ()(U val)258 constexpr T operator()(U val)
259 {
260 return m_data(val);
261 }
262 };
263 //] [/hermite_example2]
264
265
test_double()266 void test_double()
267 {
268 constexpr double radius = 2.25;
269 constexpr double c = circumference(radius);
270 constexpr double a = area(radius);
271
272 std::cout << "Circumference = " << c << std::endl;
273 std::cout << "Area = " << a << std::endl;
274
275 constexpr const_polynomial<double, 2> pa = {3, 4};
276 constexpr const_polynomial<double, 2> pb = {5, 6};
277 static_assert(pa[0] == 3);
278 static_assert(pa[1] == 4);
279 constexpr auto pc = pa * 2;
280 static_assert(pc[0] == 6);
281 static_assert(pc[1] == 8);
282 constexpr auto pd = 3 * pa;
283 static_assert(pd[0] == 3 * 3);
284 static_assert(pd[1] == 4 * 3);
285 constexpr auto pe = pa + pb;
286 static_assert(pe[0] == 3 + 5);
287 static_assert(pe[1] == 4 + 6);
288 constexpr auto pf = pa - pb;
289 static_assert(pf[0] == 3 - 5);
290 static_assert(pf[1] == 4 - 6);
291 constexpr auto pg = pa * pb;
292 static_assert(pg[0] == 15);
293 static_assert(pg[1] == 38);
294 static_assert(pg[2] == 24);
295
296 constexpr hermite_polynomial<double, 2> h1;
297 static_assert(h1[0] == -2);
298 static_assert(h1[1] == 0);
299 static_assert(h1[2] == 4);
300
301 constexpr hermite_polynomial<double, 3> h3;
302 static_assert(h3[0] == 0);
303 static_assert(h3[1] == -12);
304 static_assert(h3[2] == 0);
305 static_assert(h3[3] == 8);
306
307 constexpr hermite_polynomial<double, 9> h9;
308 static_assert(h9[0] == 0);
309 static_assert(h9[1] == 30240);
310 static_assert(h9[2] == 0);
311 static_assert(h9[3] == -80640);
312 static_assert(h9[4] == 0);
313 static_assert(h9[5] == 48384);
314 static_assert(h9[6] == 0);
315 static_assert(h9[7] == -9216);
316 static_assert(h9[8] == 0);
317 static_assert(h9[9] == 512);
318
319 static_assert(h9(0.5) == 6481);
320
321 }
322
test_float128()323 void test_float128()
324 {
325 #ifdef BOOST_HAS_FLOAT128
326 //[constexpr_circle_usage
327
328 using boost::multiprecision::float128;
329
330 constexpr float128 radius = 2.25;
331 constexpr float128 c = circumference(radius);
332 constexpr float128 a = area(radius);
333
334 std::cout << "Circumference = " << c << std::endl;
335 std::cout << "Area = " << a << std::endl;
336
337 //] [/constexpr_circle_usage]
338
339
340 constexpr hermite_polynomial<float128, 2> h1;
341 static_assert(h1[0] == -2);
342 static_assert(h1[1] == 0);
343 static_assert(h1[2] == 4);
344
345 constexpr hermite_polynomial<float128, 3> h3;
346 static_assert(h3[0] == 0);
347 static_assert(h3[1] == -12);
348 static_assert(h3[2] == 0);
349 static_assert(h3[3] == 8);
350
351 //[hermite_example3
352 constexpr hermite_polynomial<float128, 9> h9;
353 //
354 // Verify that the polynomial's coefficients match the known values:
355 //
356 static_assert(h9[0] == 0);
357 static_assert(h9[1] == 30240);
358 static_assert(h9[2] == 0);
359 static_assert(h9[3] == -80640);
360 static_assert(h9[4] == 0);
361 static_assert(h9[5] == 48384);
362 static_assert(h9[6] == 0);
363 static_assert(h9[7] == -9216);
364 static_assert(h9[8] == 0);
365 static_assert(h9[9] == 512);
366 //
367 // Define an abscissa value to evaluate at:
368 constexpr float128 abscissa(0.5);
369 //
370 // Evaluate H_9(0.5) using all constexpr arithmetic, and check that it has the expected result:
371 static_assert(h9(abscissa) == 6481);
372 //]
373 #endif
374 }
375
main()376 int main()
377 {
378 test_double();
379 test_float128();
380 std::cout << "Done!" << std::endl;
381 }
382