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)
6 #include "boost/multiprecision/cpp_int.hpp"
9 template <class T
, unsigned Order
>
10 struct const_polynomial
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
{}
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
];
25 constexpr T
& operator[](std::size_t N
)
29 constexpr const T
& operator[](std::size_t N
) const
34 constexpr T
operator()(U val
) const
36 T result
= data
[Order
];
37 for (unsigned i
= Order
; i
> 0; --i
)
40 result
+= data
[i
- 1];
44 constexpr const_polynomial
<T
, Order
- 1> derivative() const
46 const_polynomial
<T
, Order
- 1> result
;
47 for (unsigned i
= 1; i
<= Order
; ++i
)
49 result
[i
- 1] = (*this)[i
] * i
;
53 constexpr const_polynomial
operator-()
55 const_polynomial
t(*this);
56 for (unsigned i
= 0; i
<= Order
; ++i
)
61 constexpr const_polynomial
& operator*=(U val
)
63 for (unsigned i
= 0; i
<= Order
; ++i
)
64 data
[i
] = data
[i
] * val
;
68 constexpr const_polynomial
& operator/=(U val
)
70 for (unsigned i
= 0; i
<= Order
; ++i
)
71 data
[i
] = data
[i
] / val
;
75 constexpr const_polynomial
& operator+=(U val
)
81 constexpr const_polynomial
& operator-=(U val
)
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
)
92 constexpr(Order1
> Order2
)
94 const_polynomial
<T
, Order1
> result(a
);
95 for (unsigned i
= 0; i
<= Order2
; ++i
)
101 const_polynomial
<T
, Order2
> result(b
);
102 for (unsigned i
= 0; i
<= Order1
; ++i
)
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
)
111 constexpr(Order1
> Order2
)
113 const_polynomial
<T
, Order1
> result(a
);
114 for (unsigned i
= 0; i
<= Order2
; ++i
)
120 const_polynomial
<T
, Order2
> result(b
);
121 for (unsigned i
= 0; i
<= Order1
; ++i
)
122 result
[i
] = a
[i
] - b
[i
];
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
)
129 const_polynomial
<T
, Order1
+ Order2
> result
;
130 for (unsigned i
= 0; i
<= Order1
; ++i
)
132 for (unsigned j
= 0; j
<= Order2
; ++j
)
134 result
[i
+ j
] += a
[i
] * b
[j
];
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
)
142 const_polynomial
<T
, Order
> result(a
);
143 for (unsigned i
= 0; i
<= Order
; ++i
)
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
)
152 const_polynomial
<T
, Order
> result(a
);
153 for (unsigned i
= 0; i
<= Order
; ++i
)
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
)
162 const_polynomial
<T
, Order
> result
;
163 for (unsigned i
= 0; i
<= Order
; ++i
)
170 template <class T
, unsigned Order
>
171 class hermite_polynomial
173 const_polynomial
<T
, Order
> m_data
;
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())
179 constexpr const const_polynomial
<T
, Order
>& data() const
183 constexpr const T
& operator[](std::size_t N
) const
188 constexpr T
operator()(U val
) const
195 class hermite_polynomial
<T
, 0>
197 const_polynomial
<T
, 0> m_data
;
200 constexpr hermite_polynomial() : m_data
{1} {}
201 constexpr const const_polynomial
<T
, 0>& data() const
205 constexpr const T
& operator[](std::size_t N
) const
210 constexpr T
operator()(U val
)
217 class hermite_polynomial
<T
, 1>
219 const_polynomial
<T
, 1> m_data
;
222 constexpr hermite_polynomial() : m_data
{0, 2} {}
223 constexpr const const_polynomial
<T
, 1>& data() const
227 constexpr const T
& operator[](std::size_t N
) const
232 constexpr T
operator()(U val
)
242 using namespace boost::multiprecision::literals
;
244 typedef boost::multiprecision::checked_int1024_t int_backend
;
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
;
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);
264 return boost::report_errors();