]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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 | #include "boost/multiprecision/cpp_int.hpp" | |
7 | #include "test.hpp" | |
8 | ||
9 | template <class T, unsigned Order> | |
10 | struct const_polynomial | |
11 | { | |
12 | public: | |
13 | T data[Order + 1]; | |
14 | ||
15 | public: | |
16 | constexpr const_polynomial(T val = 0) : data{val} {} | |
17 | constexpr const_polynomial(const const_polynomial&) = default; | |
18 | constexpr const_polynomial(const std::initializer_list<T>& init) : data{} | |
19 | { | |
20 | if (init.size() > Order + 1) | |
21 | throw std::range_error("Too many initializers in list"); | |
22 | for (unsigned i = 0; i < init.size(); ++i) | |
23 | data[i] = init.begin()[i]; | |
24 | } | |
25 | constexpr T& operator[](std::size_t N) | |
26 | { | |
27 | return data[N]; | |
28 | } | |
29 | constexpr const T& operator[](std::size_t N) const | |
30 | { | |
31 | return data[N]; | |
32 | } | |
33 | template <class U> | |
34 | constexpr T operator()(U val) const | |
35 | { | |
36 | T result = data[Order]; | |
37 | for (unsigned i = Order; i > 0; --i) | |
38 | { | |
39 | result *= val; | |
40 | result += data[i - 1]; | |
41 | } | |
42 | return result; | |
43 | } | |
44 | constexpr const_polynomial<T, Order - 1> derivative() const | |
45 | { | |
46 | const_polynomial<T, Order - 1> result; | |
47 | for (unsigned i = 1; i <= Order; ++i) | |
48 | { | |
49 | result[i - 1] = (*this)[i] * i; | |
50 | } | |
51 | return result; | |
52 | } | |
53 | constexpr const_polynomial operator-() | |
54 | { | |
55 | const_polynomial t(*this); | |
56 | for (unsigned i = 0; i <= Order; ++i) | |
57 | t[i] = -t[i]; | |
58 | return t; | |
59 | } | |
60 | template <class U> | |
61 | constexpr const_polynomial& operator*=(U val) | |
62 | { | |
63 | for (unsigned i = 0; i <= Order; ++i) | |
64 | data[i] = data[i] * val; | |
65 | return *this; | |
66 | } | |
67 | template <class U> | |
68 | constexpr const_polynomial& operator/=(U val) | |
69 | { | |
70 | for (unsigned i = 0; i <= Order; ++i) | |
71 | data[i] = data[i] / val; | |
72 | return *this; | |
73 | } | |
74 | template <class U> | |
75 | constexpr const_polynomial& operator+=(U val) | |
76 | { | |
77 | data[0] += val; | |
78 | return *this; | |
79 | } | |
80 | template <class U> | |
81 | constexpr const_polynomial& operator-=(U val) | |
82 | { | |
83 | data[0] -= val; | |
84 | return *this; | |
85 | } | |
86 | }; | |
87 | ||
88 | template <class T, unsigned Order1, unsigned Order2> | |
89 | inline constexpr const_polynomial<T, (Order1 > Order2 ? Order1 : Order2)> operator+(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b) | |
90 | { | |
91 | if | |
92 | constexpr(Order1 > Order2) | |
93 | { | |
94 | const_polynomial<T, Order1> result(a); | |
95 | for (unsigned i = 0; i <= Order2; ++i) | |
96 | result[i] += b[i]; | |
97 | return result; | |
98 | } | |
99 | else | |
100 | { | |
101 | const_polynomial<T, Order2> result(b); | |
102 | for (unsigned i = 0; i <= Order1; ++i) | |
103 | result[i] += a[i]; | |
104 | return result; | |
105 | } | |
106 | } | |
107 | template <class T, unsigned Order1, unsigned Order2> | |
108 | inline constexpr const_polynomial<T, (Order1 > Order2 ? Order1 : Order2)> operator-(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b) | |
109 | { | |
110 | if | |
111 | constexpr(Order1 > Order2) | |
112 | { | |
113 | const_polynomial<T, Order1> result(a); | |
114 | for (unsigned i = 0; i <= Order2; ++i) | |
115 | result[i] -= b[i]; | |
116 | return result; | |
117 | } | |
118 | else | |
119 | { | |
120 | const_polynomial<T, Order2> result(b); | |
121 | for (unsigned i = 0; i <= Order1; ++i) | |
122 | result[i] = a[i] - b[i]; | |
123 | return result; | |
124 | } | |
125 | } | |
126 | template <class T, unsigned Order1, unsigned Order2> | |
127 | inline constexpr const_polynomial<T, Order1 + Order2> operator*(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b) | |
128 | { | |
129 | const_polynomial<T, Order1 + Order2> result; | |
130 | for (unsigned i = 0; i <= Order1; ++i) | |
131 | { | |
132 | for (unsigned j = 0; j <= Order2; ++j) | |
133 | { | |
134 | result[i + j] += a[i] * b[j]; | |
135 | } | |
136 | } | |
137 | return result; | |
138 | } | |
139 | template <class T, unsigned Order, class U> | |
140 | inline constexpr const_polynomial<T, Order> operator*(const const_polynomial<T, Order>& a, const U& b) | |
141 | { | |
142 | const_polynomial<T, Order> result(a); | |
143 | for (unsigned i = 0; i <= Order; ++i) | |
144 | { | |
145 | result[i] *= b; | |
146 | } | |
147 | return result; | |
148 | } | |
149 | template <class U, class T, unsigned Order> | |
150 | inline constexpr const_polynomial<T, Order> operator*(const U& b, const const_polynomial<T, Order>& a) | |
151 | { | |
152 | const_polynomial<T, Order> result(a); | |
153 | for (unsigned i = 0; i <= Order; ++i) | |
154 | { | |
155 | result[i] *= b; | |
156 | } | |
157 | return result; | |
158 | } | |
159 | template <class T, unsigned Order, class U> | |
160 | inline constexpr const_polynomial<T, Order> operator/(const const_polynomial<T, Order>& a, const U& b) | |
161 | { | |
162 | const_polynomial<T, Order> result; | |
163 | for (unsigned i = 0; i <= Order; ++i) | |
164 | { | |
165 | result[i] /= b; | |
166 | } | |
167 | return result; | |
168 | } | |
169 | ||
170 | template <class T, unsigned Order> | |
171 | class hermite_polynomial | |
172 | { | |
173 | const_polynomial<T, Order> m_data; | |
174 | ||
175 | public: | |
176 | constexpr hermite_polynomial() : m_data(hermite_polynomial<T, Order - 1>().data() * const_polynomial<T, 1>{0, 2} - hermite_polynomial<T, Order - 1>().data().derivative()) | |
177 | { | |
178 | } | |
179 | constexpr const const_polynomial<T, Order>& data() const | |
180 | { | |
181 | return m_data; | |
182 | } | |
183 | constexpr const T& operator[](std::size_t N) const | |
184 | { | |
185 | return m_data[N]; | |
186 | } | |
187 | template <class U> | |
188 | constexpr T operator()(U val) const | |
189 | { | |
190 | return m_data(val); | |
191 | } | |
192 | }; | |
193 | ||
194 | template <class T> | |
195 | class hermite_polynomial<T, 0> | |
196 | { | |
197 | const_polynomial<T, 0> m_data; | |
198 | ||
199 | public: | |
200 | constexpr hermite_polynomial() : m_data{1} {} | |
201 | constexpr const const_polynomial<T, 0>& data() const | |
202 | { | |
203 | return m_data; | |
204 | } | |
205 | constexpr const T& operator[](std::size_t N) const | |
206 | { | |
207 | return m_data[N]; | |
208 | } | |
209 | template <class U> | |
210 | constexpr T operator()(U val) | |
211 | { | |
212 | return m_data(val); | |
213 | } | |
214 | }; | |
215 | ||
216 | template <class T> | |
217 | class hermite_polynomial<T, 1> | |
218 | { | |
219 | const_polynomial<T, 1> m_data; | |
220 | ||
221 | public: | |
222 | constexpr hermite_polynomial() : m_data{0, 2} {} | |
223 | constexpr const const_polynomial<T, 1>& data() const | |
224 | { | |
225 | return m_data; | |
226 | } | |
227 | constexpr const T& operator[](std::size_t N) const | |
228 | { | |
229 | return m_data[N]; | |
230 | } | |
231 | template <class U> | |
232 | constexpr T operator()(U val) | |
233 | { | |
234 | return m_data(val); | |
235 | } | |
236 | }; | |
237 | ||
238 | ||
239 | ||
240 | int main() | |
241 | { | |
242 | using namespace boost::multiprecision::literals; | |
243 | ||
244 | typedef boost::multiprecision::checked_int1024_t int_backend; | |
245 | ||
246 | // 8192 x^13 - 319488 x^11 + 4392960 x^9 - 26357760 x^7 + 69189120 x^5 - 69189120 x^3 + 17297280 x | |
247 | constexpr hermite_polynomial<int_backend, 13> h; | |
248 | ||
249 | static_assert(h[0] == 0); | |
250 | static_assert(h[1] == 17297280); | |
251 | static_assert(h[2] == 0); | |
252 | static_assert(h[3] == -69189120); | |
253 | static_assert(h[4] == 0); | |
254 | static_assert(h[5] == 69189120); | |
255 | static_assert(h[6] == 0); | |
256 | static_assert(h[7] == -26357760); | |
257 | static_assert(h[8] == 0); | |
258 | static_assert(h[9] == 4392960); | |
259 | static_assert(h[10] == 0); | |
260 | static_assert(h[11] == -319488); | |
261 | static_assert(h[12] == 0); | |
262 | static_assert(h[13] == 8192); | |
263 | ||
264 | return boost::report_errors(); | |
265 | } | |
266 |