Skip to content

Commit 86f0c7b

Browse files
MaxSagebaumMaxSagebaum
committed
Taylor polynomial propagation implementation.
1 parent 93e676f commit 86f0c7b

File tree

6 files changed

+946
-0
lines changed

6 files changed

+946
-0
lines changed

‎include/cpp2taylor.h‎

Lines changed: 389 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,389 @@
1+
2+
#ifndef CPP2TAYLOR_H_CPP2
3+
#define CPP2TAYLOR_H_CPP2
4+
5+
6+
//=== Cpp2 type declarations ====================================================
7+
8+
9+
#include "cpp2util.h"
10+
11+
#line 1 "cpp2taylor.h2"
12+
13+
#line 4 "cpp2taylor.h2"
14+
namespace cpp2 {
15+
16+
template<typename R, int dim> class taylor;
17+
18+
19+
#line 192 "cpp2taylor.h2"
20+
}
21+
22+
23+
//=== Cpp2 type definitions and function declarations ===========================
24+
25+
#line 1 "cpp2taylor.h2"
26+
#ifndef CPP2_CPP2TAYLOR_H
27+
#define CPP2_CPP2TAYLOR_H
28+
29+
#line 4 "cpp2taylor.h2"
30+
namespace cpp2 {
31+
32+
template<typename R, int dim> class taylor {
33+
private: std::array<R,dim> v {};
34+
35+
public: explicit taylor();
36+
public: taylor(R const& d1);
37+
#line 10 "cpp2taylor.h2"
38+
public: auto operator=(R const& d1) -> taylor& ;
39+
40+
#line 13 "cpp2taylor.h2"
41+
public: taylor(taylor const& that);
42+
#line 13 "cpp2taylor.h2"
43+
public: auto operator=(taylor const& that) -> taylor& ;
44+
#line 13 "cpp2taylor.h2"
45+
public: taylor(taylor&& that) noexcept;
46+
#line 13 "cpp2taylor.h2"
47+
public: auto operator=(taylor&& that) noexcept -> taylor& ;
48+
49+
// C++ interface
50+
51+
public: [[nodiscard]] auto operator[](cpp2::impl::in<int> k) const& -> R;
52+
53+
#line 27 "cpp2taylor.h2"
54+
public: [[nodiscard]] auto operator[](cpp2::impl::in<int> i) & -> auto&&;
55+
56+
#line 32 "cpp2taylor.h2"
57+
// C++2 interface / AD interface
58+
59+
public: [[nodiscard]] auto get(cpp2::impl::in<int> i, cpp2::impl::in<R> v0) const& -> R;
60+
61+
#line 44 "cpp2taylor.h2"
62+
public: [[nodiscard]] auto add(cpp2::impl::in<taylor> o, cpp2::impl::in<R> v0, cpp2::impl::in<R> o0) const& -> taylor;
63+
64+
#line 55 "cpp2taylor.h2"
65+
public: [[nodiscard]] auto sub(cpp2::impl::in<taylor> o, cpp2::impl::in<R> v0, cpp2::impl::in<R> o0) const& -> taylor;
66+
67+
#line 66 "cpp2taylor.h2"
68+
public: [[nodiscard]] auto mul(cpp2::impl::in<taylor> o, cpp2::impl::in<R> v0, cpp2::impl::in<R> o0) const& -> taylor;
69+
70+
#line 79 "cpp2taylor.h2"
71+
public: [[nodiscard]] auto div(cpp2::impl::in<taylor> o, cpp2::impl::in<R> v0, cpp2::impl::in<R> o0) const& -> taylor;
72+
73+
#line 97 "cpp2taylor.h2"
74+
public: [[nodiscard]] auto sqrt(cpp2::impl::in<R> v0) const& -> taylor;
75+
76+
#line 116 "cpp2taylor.h2"
77+
public: [[nodiscard]] auto log(cpp2::impl::in<R> v0) const& -> taylor;
78+
79+
#line 135 "cpp2taylor.h2"
80+
public: [[nodiscard]] auto exp(cpp2::impl::in<R> v0) const& -> taylor;
81+
82+
#line 154 "cpp2taylor.h2"
83+
public: static auto comp_sin_cos(taylor& s, taylor& c, cpp2::impl::in<taylor> u, cpp2::impl::in<R> u0) -> void;
84+
85+
#line 171 "cpp2taylor.h2"
86+
public: [[nodiscard]] auto sin(cpp2::impl::in<R> v0) const& -> taylor;
87+
88+
#line 181 "cpp2taylor.h2"
89+
public: [[nodiscard]] auto cos(cpp2::impl::in<R> v0) const& -> taylor;
90+
91+
#line 190 "cpp2taylor.h2"
92+
};
93+
94+
} // cpp2 namespace
95+
96+
#endif // CPP2_CPP2TAYLOR_H
97+
98+
99+
//=== Cpp2 function definitions =================================================
100+
101+
#line 1 "cpp2taylor.h2"
102+
103+
#line 4 "cpp2taylor.h2"
104+
namespace cpp2 {
105+
106+
#line 9 "cpp2taylor.h2"
107+
template <typename R, int dim> taylor<R,dim>::taylor(){}
108+
#line 10 "cpp2taylor.h2"
109+
template <typename R, int dim> taylor<R,dim>::taylor(R const& d1)
110+
: v{ d1 }{
111+
112+
#line 12 "cpp2taylor.h2"
113+
}
114+
#line 10 "cpp2taylor.h2"
115+
template <typename R, int dim> auto taylor<R,dim>::operator=(R const& d1) -> taylor& {
116+
v = d1;
117+
return *this;
118+
119+
#line 12 "cpp2taylor.h2"
120+
}
121+
#line 13 "cpp2taylor.h2"
122+
template <typename R, int dim> taylor<R,dim>::taylor(taylor const& that)
123+
: v{ that.v }{}
124+
#line 13 "cpp2taylor.h2"
125+
template <typename R, int dim> auto taylor<R,dim>::operator=(taylor const& that) -> taylor& {
126+
v = that.v;
127+
return *this; }
128+
#line 13 "cpp2taylor.h2"
129+
template <typename R, int dim> taylor<R,dim>::taylor(taylor&& that) noexcept
130+
: v{ std::move(that).v }{}
131+
#line 13 "cpp2taylor.h2"
132+
template <typename R, int dim> auto taylor<R,dim>::operator=(taylor&& that) noexcept -> taylor& {
133+
v = std::move(that).v;
134+
return *this; }
135+
136+
#line 17 "cpp2taylor.h2"
137+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::operator[](cpp2::impl::in<int> k) const& -> R{
138+
if (cpp2::cpp2_default.is_active() && !([_0 = 1, _1 = k, _2 = dim]{ return cpp2::impl::cmp_less_eq(_0,_1) && cpp2::impl::cmp_less_eq(_1,_2); }()) ) { cpp2::cpp2_default.report_violation(""); }
139+
R r {CPP2_ASSERT_IN_BOUNDS(v, k - 1)};
140+
{
141+
auto i{2};
142+
143+
#line 21 "cpp2taylor.h2"
144+
for( ; cpp2::impl::cmp_less_eq(i,k); i += 1 ) {
145+
r *= i;
146+
}
147+
}
148+
#line 24 "cpp2taylor.h2"
149+
return r;
150+
}
151+
152+
#line 27 "cpp2taylor.h2"
153+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::operator[](cpp2::impl::in<int> i) & -> auto&&{
154+
if (cpp2::cpp2_default.is_active() && !([_0 = 1, _1 = 1, _2 = dim]{ return cpp2::impl::cmp_less_eq(_0,_1) && cpp2::impl::cmp_less_eq(_1,_2); }()) ) { cpp2::cpp2_default.report_violation(""); }
155+
return CPP2_ASSERT_IN_BOUNDS(v, i - 1);
156+
}
157+
158+
#line 34 "cpp2taylor.h2"
159+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::get(cpp2::impl::in<int> i, cpp2::impl::in<R> v0) const& -> R{
160+
if (cpp2::cpp2_default.is_active() && !([_0 = 0, _1 = i, _2 = dim]{ return cpp2::impl::cmp_less_eq(_0,_1) && cpp2::impl::cmp_less_eq(_1,_2); }()) ) { cpp2::cpp2_default.report_violation(""); }
161+
162+
if (i == 0) {
163+
return v0;
164+
}
165+
166+
return CPP2_ASSERT_IN_BOUNDS(v, i - 1);
167+
}
168+
169+
#line 44 "cpp2taylor.h2"
170+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::add(cpp2::impl::in<taylor> o, cpp2::impl::in<R> v0, cpp2::impl::in<R> o0) const& -> taylor{
171+
taylor r {};
172+
{
173+
auto k{1};
174+
175+
#line 48 "cpp2taylor.h2"
176+
for( ; cpp2::impl::cmp_less_eq(k,dim); k += 1 ) {
177+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) = get(k, v0) + CPP2_UFCS(get)(o, k, o0);
178+
}
179+
}
180+
181+
#line 52 "cpp2taylor.h2"
182+
return r;
183+
}
184+
185+
#line 55 "cpp2taylor.h2"
186+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::sub(cpp2::impl::in<taylor> o, cpp2::impl::in<R> v0, cpp2::impl::in<R> o0) const& -> taylor{
187+
taylor r {};
188+
{
189+
auto k{1};
190+
191+
#line 59 "cpp2taylor.h2"
192+
for( ; cpp2::impl::cmp_less_eq(k,dim); k += 1 ) {
193+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) = get(k, v0) - CPP2_UFCS(get)(o, k, o0);
194+
}
195+
}
196+
197+
#line 63 "cpp2taylor.h2"
198+
return r;
199+
}
200+
201+
#line 66 "cpp2taylor.h2"
202+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::mul(cpp2::impl::in<taylor> o, cpp2::impl::in<R> v0, cpp2::impl::in<R> o0) const& -> taylor{
203+
taylor r {};
204+
{
205+
auto k{1};
206+
207+
#line 70 "cpp2taylor.h2"
208+
for( ; cpp2::impl::cmp_less_eq(k,dim); k += 1 ) {
209+
{
210+
auto j{0};
211+
212+
#line 72 "cpp2taylor.h2"
213+
for( ; cpp2::impl::cmp_less_eq(j,k); j += 1 ) {
214+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) += get(j, v0) * o.get(k - j, o0);
215+
}
216+
}
217+
#line 75 "cpp2taylor.h2"
218+
}
219+
}
220+
#line 76 "cpp2taylor.h2"
221+
return r;
222+
}
223+
224+
#line 79 "cpp2taylor.h2"
225+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::div(cpp2::impl::in<taylor> o, cpp2::impl::in<R> v0, cpp2::impl::in<R> o0) const& -> taylor{
226+
taylor r {};
227+
R r0 {v0 / CPP2_ASSERT_NOT_ZERO(CPP2_TYPEOF(v0),o0)};
228+
229+
R factor {1.0 / CPP2_ASSERT_NOT_ZERO(CPP2_TYPEOF(1.0),o0)};
230+
{
231+
auto k{1};
232+
233+
#line 86 "cpp2taylor.h2"
234+
for( ; cpp2::impl::cmp_less_eq(k,dim); k += 1 ) {
235+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) = get(k, v0);
236+
{
237+
auto j{0};
238+
239+
#line 89 "cpp2taylor.h2"
240+
for( ; cpp2::impl::cmp_less(j,k); j += 1 ) {
241+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) -= CPP2_UFCS(get)(r, j, r0) * o.get(k - j, o0);
242+
}
243+
}
244+
#line 92 "cpp2taylor.h2"
245+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) *= factor;
246+
}
247+
}
248+
#line 94 "cpp2taylor.h2"
249+
return r;
250+
}
251+
252+
#line 97 "cpp2taylor.h2"
253+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::sqrt(cpp2::impl::in<R> v0) const& -> taylor{
254+
taylor r {};
255+
R r0 {std::sqrt(v0)};
256+
257+
R factor {0.5 / CPP2_ASSERT_NOT_ZERO(CPP2_TYPEOF(0.5),r0)};
258+
{
259+
auto k{1};
260+
261+
#line 104 "cpp2taylor.h2"
262+
for( ; cpp2::impl::cmp_less_eq(k,dim); k += 1 ) {
263+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) = get(k, v0);
264+
{
265+
auto j{1};
266+
267+
#line 107 "cpp2taylor.h2"
268+
for( ; cpp2::impl::cmp_less(j,k); j += 1 ) {
269+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) -= r.get(j, r0) * r.get(k - j, r0);
270+
}
271+
}
272+
#line 110 "cpp2taylor.h2"
273+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) *= factor;
274+
}
275+
}
276+
277+
#line 113 "cpp2taylor.h2"
278+
return r;
279+
}
280+
281+
#line 116 "cpp2taylor.h2"
282+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::log(cpp2::impl::in<R> v0) const& -> taylor{
283+
taylor r {};
284+
R r0 {std::log(v0)};
285+
286+
R factor {1.0 / CPP2_ASSERT_NOT_ZERO(CPP2_TYPEOF(1.0),v0)};
287+
{
288+
auto k{1};
289+
290+
#line 123 "cpp2taylor.h2"
291+
for( ; cpp2::impl::cmp_less_eq(k,dim); k += 1 ) {
292+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) = k * get(k, v0);
293+
{
294+
auto j{1};
295+
296+
#line 126 "cpp2taylor.h2"
297+
for( ; cpp2::impl::cmp_less(j,k); j += 1 ) {
298+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) -= j * get(k - j, v0) * r.get(j, r0);
299+
}
300+
}
301+
#line 129 "cpp2taylor.h2"
302+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) *= factor / CPP2_ASSERT_NOT_ZERO(CPP2_TYPEOF(factor),k);
303+
}
304+
}
305+
306+
#line 132 "cpp2taylor.h2"
307+
return r;
308+
}
309+
310+
#line 135 "cpp2taylor.h2"
311+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::exp(cpp2::impl::in<R> v0) const& -> taylor{
312+
taylor r {};
313+
R r0 {std::exp(v0)};
314+
315+
R factor {1.0 / CPP2_ASSERT_NOT_ZERO(CPP2_TYPEOF(1.0),v0)};
316+
{
317+
auto k{1};
318+
319+
#line 142 "cpp2taylor.h2"
320+
for( ; cpp2::impl::cmp_less_eq(k,dim); k += 1 ) {
321+
{
322+
auto j{1};
323+
324+
#line 144 "cpp2taylor.h2"
325+
for( ; cpp2::impl::cmp_less_eq(j,k); j += 1 ) {
326+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) += j * r.get(k - j, r0) * get(j, v0);
327+
}
328+
}
329+
#line 147 "cpp2taylor.h2"
330+
CPP2_ASSERT_IN_BOUNDS(r.v, k - 1) /= CPP2_ASSERT_NOT_ZERO(CPP2_TYPEOF(CPP2_ASSERT_IN_BOUNDS(r.v, k - 1)),k);
331+
}
332+
}
333+
334+
#line 150 "cpp2taylor.h2"
335+
return r;
336+
}
337+
338+
#line 154 "cpp2taylor.h2"
339+
template <typename R, int dim> auto taylor<R,dim>::comp_sin_cos(taylor& s, taylor& c, cpp2::impl::in<taylor> u, cpp2::impl::in<R> u0) -> void{
340+
R s0 {std::sin(u0)};
341+
R c0 {std::cos(u0)};
342+
{
343+
auto k{1};
344+
345+
#line 159 "cpp2taylor.h2"
346+
for( ; cpp2::impl::cmp_less_eq(k,dim); k += 1 ) {
347+
{
348+
auto j{1};
349+
350+
#line 161 "cpp2taylor.h2"
351+
for( ; cpp2::impl::cmp_less_eq(j,k); j += 1 ) {
352+
CPP2_ASSERT_IN_BOUNDS(s.v, k - 1) += j * u.get(j, u0) * CPP2_UFCS(get)(c, k - j, c0);
353+
CPP2_ASSERT_IN_BOUNDS(c.v, k - 1) -= j * u.get(j, u0) * CPP2_UFCS(get)(s, k - j, s0);
354+
}
355+
}
356+
#line 165 "cpp2taylor.h2"
357+
CPP2_ASSERT_IN_BOUNDS(s.v, k - 1) /= CPP2_ASSERT_NOT_ZERO(CPP2_TYPEOF(CPP2_ASSERT_IN_BOUNDS(s.v, k - 1)),k);
358+
CPP2_ASSERT_IN_BOUNDS(c.v, k - 1) /= CPP2_ASSERT_NOT_ZERO(CPP2_TYPEOF(CPP2_ASSERT_IN_BOUNDS(c.v, k - 1)),k);
359+
}
360+
}
361+
#line 168 "cpp2taylor.h2"
362+
}
363+
364+
#line 171 "cpp2taylor.h2"
365+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::sin(cpp2::impl::in<R> v0) const& -> taylor{
366+
taylor t {};
367+
taylor r {};
368+
369+
comp_sin_cos(r, t, (*this), v0);
370+
static_cast<void>(cpp2::move(t));
371+
372+
return r;
373+
}
374+
375+
#line 181 "cpp2taylor.h2"
376+
template <typename R, int dim> [[nodiscard]] auto taylor<R,dim>::cos(cpp2::impl::in<R> v0) const& -> taylor{
377+
taylor t {};
378+
taylor r {};
379+
380+
comp_sin_cos(t, r, (*this), v0);
381+
static_cast<void>(cpp2::move(t));
382+
383+
return r;
384+
}
385+
386+
#line 192 "cpp2taylor.h2"
387+
}
388+
389+
#endif

0 commit comments

Comments
 (0)