1 // Copyright 2020, 2021 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com)
2 //
3 // This file is part of the heyoka library.
4 //
5 // This Source Code Form is subject to the terms of the Mozilla
6 // Public License v. 2.0. If a copy of the MPL was not distributed
7 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
9 #include <heyoka/config.hpp>
10
11 #include <algorithm>
12 #include <cmath>
13 #include <initializer_list>
14 #include <limits>
15 #include <random>
16 #include <tuple>
17 #include <utility>
18 #include <vector>
19
20 #include <boost/cstdint.hpp>
21 #include <boost/math/constants/constants.hpp>
22 #include <boost/math/tools/roots.hpp>
23
24 #if defined(HEYOKA_HAVE_REAL128)
25
26 #include <mp++/real128.hpp>
27
28 #endif
29
30 #include <heyoka/expression.hpp>
31 #include <heyoka/llvm_state.hpp>
32 #include <heyoka/math/kepE.hpp>
33 #include <heyoka/number.hpp>
34 #include <heyoka/taylor.hpp>
35
36 #include "catch.hpp"
37 #include "test_utils.hpp"
38
39 static std::mt19937 rng;
40
41 using namespace heyoka;
42 using namespace heyoka_test;
43
44 const auto fp_types = std::tuple<double
45 #if !defined(HEYOKA_ARCH_PPC)
46 ,
47 long double
48 #endif
49 #if defined(HEYOKA_HAVE_REAL128)
50 ,
51 mppp::real128
52 #endif
53 >{};
54
55 // Boost-based Kepler solver.
__anon9f23fef20102(auto ecc, auto M) 56 auto bmt_inv_kep_E = [](auto ecc, auto M) {
57 using std::sin;
58 using std::cos;
59
60 using fp_t = decltype(ecc);
61
62 // Initial guess.
63 auto ig = ecc < 0.8 ? M : static_cast<fp_t>(boost::math::constants::pi<double>());
64
65 auto func = [ecc, M](auto E) { return std::make_pair(E - ecc * sin(E) - M, 1 - ecc * cos(E)); };
66
67 boost::uintmax_t max_iter = 50;
68
69 return boost::math::tools::newton_raphson_iterate(func, ig, fp_t(0), fp_t(2 * boost::math::constants::pi<double>()),
70 std::numeric_limits<fp_t>::digits - 2, max_iter);
71 };
72
73 template <typename T, typename U>
compare_batch_scalar(std::initializer_list<U> sys,unsigned opt_level,bool high_accuracy,bool compact_mode)74 void compare_batch_scalar(std::initializer_list<U> sys, unsigned opt_level, bool high_accuracy, bool compact_mode)
75 {
76 for (auto batch_size : {2u, 4u, 8u, 5u}) {
77 llvm_state s{kw::opt_level = opt_level};
78
79 taylor_add_jet<T>(s, "jet_batch", sys, 3, batch_size, high_accuracy, compact_mode);
80 taylor_add_jet<T>(s, "jet_scalar", sys, 3, 1, high_accuracy, compact_mode);
81
82 s.compile();
83
84 auto jptr_batch = reinterpret_cast<void (*)(T *, const T *, const T *)>(s.jit_lookup("jet_batch"));
85 auto jptr_scalar = reinterpret_cast<void (*)(T *, const T *, const T *)>(s.jit_lookup("jet_scalar"));
86
87 std::vector<T> jet_batch;
88 jet_batch.resize(8 * batch_size);
89 std::uniform_real_distribution<float> dist(.1f, .9f);
90 std::generate(jet_batch.begin(), jet_batch.end(), [&dist]() { return T{dist(rng)}; });
91
92 std::vector<T> jet_scalar;
93 jet_scalar.resize(8);
94
95 jptr_batch(jet_batch.data(), nullptr, nullptr);
96
97 for (auto batch_idx = 0u; batch_idx < batch_size; ++batch_idx) {
98 // Assign the initial values of x and y.
99 for (auto i = 0u; i < 2u; ++i) {
100 jet_scalar[i] = jet_batch[i * batch_size + batch_idx];
101 }
102
103 jptr_scalar(jet_scalar.data(), nullptr, nullptr);
104
105 for (auto i = 2u; i < 8u; ++i) {
106 REQUIRE(jet_scalar[i] == approximately(jet_batch[i * batch_size + batch_idx], T(10000)));
107 }
108 }
109 }
110 }
111
112 // Issue in the decomposition when e = 0.
113 TEST_CASE("taylor kepE decompose bug 00")
114 {
115 llvm_state s;
116
117 auto M = "M"_var;
118
119 taylor_add_jet<double>(s, "jet", {kepE(0_dbl, M)}, 1, 1, false, false);
120 }
121
122 TEST_CASE("taylor kepE")
123 {
__anon9f23fef20402(auto fp_x, unsigned opt_level, bool high_accuracy, bool compact_mode) 124 auto tester = [](auto fp_x, unsigned opt_level, bool high_accuracy, bool compact_mode) {
125 using std::sin;
126 using std::cos;
127
128 using fp_t = decltype(fp_x);
129
130 using Catch::Matchers::Message;
131
132 auto x = "x"_var, y = "y"_var;
133
134 const auto a = fp_t{1} / 3;
135 const auto b = 1 + fp_t{3} / fp_t{7};
136 const auto c = fp_t{2} / 7;
137
138 // Number-number tests.
139 {
140 llvm_state s{kw::opt_level = opt_level};
141
142 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, expression{number{b}}), x + y}, 1, 1,
143 high_accuracy, compact_mode);
144
145 s.compile();
146
147 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
148
149 std::vector<fp_t> jet{fp_t{2}, fp_t{3}};
150 jet.resize(4);
151
152 jptr(jet.data(), nullptr, nullptr);
153
154 REQUIRE(jet[0] == 2);
155 REQUIRE(jet[1] == 3);
156 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(a, b)));
157 REQUIRE(jet[3] == 5);
158 }
159
160 {
161 llvm_state s{kw::opt_level = opt_level};
162
163 taylor_add_jet<fp_t>(s, "jet", {kepE(par[0], expression{number{b}}), x + y}, 1, 1, high_accuracy,
164 compact_mode);
165
166 s.compile();
167
168 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
169
170 std::vector<fp_t> jet{fp_t{2}, fp_t{3}};
171 jet.resize(4);
172
173 std::vector<fp_t> pars{a};
174
175 jptr(jet.data(), pars.data(), nullptr);
176
177 REQUIRE(jet[0] == 2);
178 REQUIRE(jet[1] == 3);
179 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(a, b)));
180 REQUIRE(jet[3] == 5);
181 }
182
183 {
184 llvm_state s{kw::opt_level = opt_level};
185
186 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, expression{number{b}}), x + y}, 1, 2,
187 high_accuracy, compact_mode);
188
189 s.compile();
190
191 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
192
193 std::vector<fp_t> jet{fp_t{2}, fp_t{-1}, fp_t{3}, fp_t{5}};
194 jet.resize(8);
195
196 jptr(jet.data(), nullptr, nullptr);
197
198 REQUIRE(jet[0] == 2);
199 REQUIRE(jet[1] == -1);
200
201 REQUIRE(jet[2] == 3);
202 REQUIRE(jet[3] == 5);
203
204 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(a, b)));
205 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(a, b)));
206
207 REQUIRE(jet[6] == 5);
208 REQUIRE(jet[7] == 4);
209 }
210
211 {
212 llvm_state s{kw::opt_level = opt_level};
213
214 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, par[1]), x + y}, 1, 2, high_accuracy,
215 compact_mode);
216
217 s.compile();
218
219 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
220
221 std::vector<fp_t> jet{fp_t{2}, fp_t{-1}, fp_t{3}, fp_t{5}};
222 jet.resize(8);
223
224 std::vector<fp_t> pars{fp_t{0}, fp_t{0}, b, b};
225
226 jptr(jet.data(), pars.data(), nullptr);
227
228 REQUIRE(jet[0] == 2);
229 REQUIRE(jet[1] == -1);
230
231 REQUIRE(jet[2] == 3);
232 REQUIRE(jet[3] == 5);
233
234 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(a, b)));
235 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(a, b)));
236
237 REQUIRE(jet[6] == 5);
238 REQUIRE(jet[7] == 4);
239 }
240
241 {
242 llvm_state s{kw::opt_level = opt_level};
243
244 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, expression{number{b}}), x + y}, 2, 1,
245 high_accuracy, compact_mode);
246
247 s.compile();
248
249 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
250
251 std::vector<fp_t> jet{fp_t{2}, fp_t{3}};
252 jet.resize(6);
253
254 jptr(jet.data(), nullptr, nullptr);
255
256 REQUIRE(jet[0] == 2);
257 REQUIRE(jet[1] == 3);
258 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(a, b)));
259 REQUIRE(jet[3] == 5);
260 REQUIRE(jet[4] == 0);
261 REQUIRE(jet[5] == fp_t{1} / 2 * (jet[2] + jet[3]));
262 }
263
264 {
265 llvm_state s{kw::opt_level = opt_level};
266
267 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, expression{number{b}}), x + y}, 2, 2,
268 high_accuracy, compact_mode);
269
270 s.compile();
271
272 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
273
274 std::vector<fp_t> jet{fp_t{2}, fp_t{-1}, fp_t{3}, fp_t{5}};
275 jet.resize(12);
276
277 jptr(jet.data(), nullptr, nullptr);
278
279 REQUIRE(jet[0] == 2);
280 REQUIRE(jet[1] == -1);
281
282 REQUIRE(jet[2] == 3);
283 REQUIRE(jet[3] == 5);
284
285 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(a, b)));
286 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(a, b)));
287
288 REQUIRE(jet[6] == 5);
289 REQUIRE(jet[7] == 4);
290
291 REQUIRE(jet[8] == 0);
292 REQUIRE(jet[9] == 0);
293
294 REQUIRE(jet[10] == fp_t{1} / 2 * (jet[4] + jet[6]));
295 REQUIRE(jet[11] == fp_t{1} / 2 * (jet[5] + jet[7]));
296 }
297
298 {
299 llvm_state s{kw::opt_level = opt_level};
300
301 taylor_add_jet<fp_t>(s, "jet", {kepE(par[0], par[1]), x + y}, 3, 3, high_accuracy, compact_mode);
302
303 s.compile();
304
305 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
306
307 std::vector<fp_t> jet{fp_t{2}, fp_t{-1}, fp_t{-4}, fp_t{3}, fp_t{5}, fp_t{6}};
308 jet.resize(24);
309
310 std::vector<fp_t> pars{a, a, a, b, b, b};
311
312 jptr(jet.data(), pars.data(), nullptr);
313
314 REQUIRE(jet[0] == 2);
315 REQUIRE(jet[1] == -1);
316 REQUIRE(jet[2] == -4);
317
318 REQUIRE(jet[3] == 3);
319 REQUIRE(jet[4] == 5);
320 REQUIRE(jet[5] == 6);
321
322 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(a, b)));
323 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(a, b)));
324 REQUIRE(jet[8] == approximately(bmt_inv_kep_E(a, b)));
325
326 REQUIRE(jet[9] == 5);
327 REQUIRE(jet[10] == 4);
328 REQUIRE(jet[11] == 2);
329
330 REQUIRE(jet[12] == 0);
331 REQUIRE(jet[13] == 0);
332 REQUIRE(jet[14] == 0);
333
334 REQUIRE(jet[15] == fp_t{1} / 2 * (jet[6] + jet[9]));
335 REQUIRE(jet[16] == fp_t{1} / 2 * (jet[7] + jet[10]));
336 REQUIRE(jet[17] == fp_t{1} / 2 * (jet[8] + jet[11]));
337
338 REQUIRE(jet[18] == 0);
339 REQUIRE(jet[19] == 0);
340 REQUIRE(jet[20] == 0);
341
342 REQUIRE(jet[21] == approximately(fp_t{1} / 6 * (2 * jet[12] + 2 * jet[15])));
343 REQUIRE(jet[22] == approximately(fp_t{1} / 6 * (2 * jet[13] + 2 * jet[16])));
344 REQUIRE(jet[23] == approximately(fp_t{1} / 6 * (2 * jet[14] + 2 * jet[17])));
345 }
346
347 // Do the batch/scalar comparison.
348 compare_batch_scalar<fp_t>({kepE(expression{number{a}}, expression{number{b}}), x + y}, opt_level,
349 high_accuracy, compact_mode);
350
351 // Variable-number tests.
352 {
353 llvm_state s{kw::opt_level = opt_level};
354
355 taylor_add_jet<fp_t>(s, "jet", {kepE(y, expression{number{a}}), kepE(x, expression{number{b}})}, 1, 1,
356 high_accuracy, compact_mode);
357
358 s.compile();
359
360 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
361
362 std::vector<fp_t> jet{fp_t{.2}, fp_t{.3}};
363 jet.resize(4);
364
365 jptr(jet.data(), nullptr, nullptr);
366
367 REQUIRE(jet[0] == fp_t{.2});
368 REQUIRE(jet[1] == fp_t{.3});
369 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(jet[1], a)));
370 REQUIRE(jet[3] == approximately(bmt_inv_kep_E(jet[0], b)));
371 }
372
373 {
374 llvm_state s{kw::opt_level = opt_level};
375
376 taylor_add_jet<fp_t>(s, "jet", {kepE(y, par[0]), kepE(x, expression{number{b}})}, 1, 1, high_accuracy,
377 compact_mode);
378
379 s.compile();
380
381 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
382
383 std::vector<fp_t> jet{fp_t{.2}, fp_t{.3}};
384 jet.resize(4);
385
386 std::vector<fp_t> pars{a};
387
388 jptr(jet.data(), pars.data(), nullptr);
389
390 REQUIRE(jet[0] == fp_t{.2});
391 REQUIRE(jet[1] == fp_t{.3});
392 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(jet[1], a)));
393 REQUIRE(jet[3] == approximately(bmt_inv_kep_E(jet[0], b)));
394 }
395
396 {
397 llvm_state s{kw::opt_level = opt_level};
398
399 taylor_add_jet<fp_t>(s, "jet", {kepE(y, expression{number{b}}), kepE(x, expression{number{b}})}, 1, 2,
400 high_accuracy, compact_mode);
401
402 s.compile();
403
404 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
405
406 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.3}, fp_t{.4}};
407 jet.resize(8);
408
409 jptr(jet.data(), nullptr, nullptr);
410
411 REQUIRE(jet[0] == .2);
412 REQUIRE(jet[1] == .5);
413
414 REQUIRE(jet[2] == .3);
415 REQUIRE(jet[3] == .4);
416
417 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(jet[2], b)));
418 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(jet[3], b)));
419
420 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(jet[0], b)));
421 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(jet[1], b)));
422 }
423
424 {
425 llvm_state s{kw::opt_level = opt_level};
426
427 taylor_add_jet<fp_t>(s, "jet", {kepE(y, expression{number{b}}), kepE(x, par[1])}, 1, 2, high_accuracy,
428 compact_mode);
429
430 s.compile();
431
432 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
433
434 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.3}, fp_t{.4}};
435 jet.resize(8);
436
437 std::vector<fp_t> pars{fp_t{0}, fp_t{0}, b, b};
438
439 jptr(jet.data(), pars.data(), nullptr);
440
441 REQUIRE(jet[0] == .2);
442 REQUIRE(jet[1] == .5);
443
444 REQUIRE(jet[2] == .3);
445 REQUIRE(jet[3] == .4);
446
447 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(jet[2], b)));
448 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(jet[3], b)));
449
450 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(jet[0], b)));
451 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(jet[1], b)));
452 }
453
454 {
455 llvm_state s{kw::opt_level = opt_level};
456
457 taylor_add_jet<fp_t>(s, "jet", {kepE(y, expression{number{b}}), kepE(x, expression{number{b}})}, 2, 1,
458 high_accuracy, compact_mode);
459
460 s.compile();
461
462 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
463
464 std::vector<fp_t> jet{fp_t{.2}, fp_t{.3}};
465 jet.resize(6);
466
467 jptr(jet.data(), nullptr, nullptr);
468
469 REQUIRE(jet[0] == .2);
470 REQUIRE(jet[1] == .3);
471 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(jet[1], b)));
472 REQUIRE(jet[3] == approximately(bmt_inv_kep_E(jet[0], b)));
473 REQUIRE(jet[4] == approximately(fp_t{1} / 2 * jet[3] * sin(jet[2]) / (1 - jet[1] * cos(jet[2]))));
474 REQUIRE(jet[5] == approximately(fp_t{1} / 2 * jet[2] * sin(jet[3]) / (1 - jet[0] * cos(jet[3]))));
475 }
476
477 {
478 llvm_state s{kw::opt_level = opt_level};
479
480 taylor_add_jet<fp_t>(s, "jet", {kepE(y, expression{number{b}}), kepE(x, expression{number{b}})}, 2, 2,
481 high_accuracy, compact_mode);
482
483 s.compile();
484
485 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
486
487 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.3}, fp_t{.4}};
488 jet.resize(12);
489
490 jptr(jet.data(), nullptr, nullptr);
491
492 REQUIRE(jet[0] == .2);
493 REQUIRE(jet[1] == .5);
494
495 REQUIRE(jet[2] == .3);
496 REQUIRE(jet[3] == .4);
497
498 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(jet[2], b)));
499 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(jet[3], b)));
500
501 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(jet[0], b)));
502 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(jet[1], b)));
503
504 REQUIRE(jet[8] == approximately(fp_t{1} / 2 * jet[6] * sin(jet[4]) / (1 - jet[2] * cos(jet[4]))));
505 REQUIRE(jet[9] == approximately(fp_t{1} / 2 * jet[7] * sin(jet[5]) / (1 - jet[3] * cos(jet[5]))));
506
507 REQUIRE(jet[10] == approximately(fp_t{1} / 2 * jet[4] * sin(jet[6]) / (1 - jet[0] * cos(jet[6]))));
508 REQUIRE(jet[11] == approximately(fp_t{1} / 2 * jet[5] * sin(jet[7]) / (1 - jet[1] * cos(jet[7]))));
509 }
510
511 {
512 llvm_state s{kw::opt_level = opt_level};
513
514 taylor_add_jet<fp_t>(s, "jet", {kepE(y, expression{number{b}}), kepE(x, expression{number{b}})}, 3, 3,
515 high_accuracy, compact_mode);
516
517 s.compile();
518
519 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
520
521 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.1}, fp_t{.3}, fp_t{.4}, fp_t{.6}};
522 jet.resize(24);
523
524 jptr(jet.data(), nullptr, nullptr);
525
526 REQUIRE(jet[0] == .2);
527 REQUIRE(jet[1] == .5);
528 REQUIRE(jet[2] == .1);
529
530 REQUIRE(jet[3] == .3);
531 REQUIRE(jet[4] == .4);
532 REQUIRE(jet[5] == .6);
533
534 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(jet[3], b)));
535 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(jet[4], b)));
536 REQUIRE(jet[8] == approximately(bmt_inv_kep_E(jet[5], b)));
537
538 REQUIRE(jet[9] == approximately(bmt_inv_kep_E(jet[0], b)));
539 REQUIRE(jet[10] == approximately(bmt_inv_kep_E(jet[1], b)));
540 REQUIRE(jet[11] == approximately(bmt_inv_kep_E(jet[2], b)));
541
542 REQUIRE(jet[12] == approximately(fp_t{1} / 2 * jet[9] * sin(jet[6]) / (1 - jet[3] * cos(jet[6]))));
543 REQUIRE(jet[13] == approximately(fp_t{1} / 2 * jet[10] * sin(jet[7]) / (1 - jet[4] * cos(jet[7]))));
544 REQUIRE(jet[14] == approximately(fp_t{1} / 2 * jet[11] * sin(jet[8]) / (1 - jet[5] * cos(jet[8]))));
545
546 REQUIRE(jet[15] == approximately(fp_t{1} / 2 * jet[6] * sin(jet[9]) / (1 - jet[0] * cos(jet[9]))));
547 REQUIRE(jet[16] == approximately(fp_t{1} / 2 * jet[7] * sin(jet[10]) / (1 - jet[1] * cos(jet[10]))));
548 REQUIRE(jet[17] == approximately(fp_t{1} / 2 * jet[8] * sin(jet[11]) / (1 - jet[2] * cos(jet[11]))));
549
550 REQUIRE(jet[18]
551 == approximately(
552 fp_t{1} / 6
553 * (2 * jet[15] * sin(jet[6]) / (1 - jet[3] * cos(jet[6]))
554 + jet[9]
555 * (2 * jet[12] * cos(jet[6]) * (1 - jet[3] * cos(jet[6]))
556 + sin(jet[6]) * (jet[9] * cos(jet[6]) - jet[3] * 2 * jet[12] * sin(jet[6])))
557 / ((1 - jet[3] * cos(jet[6])) * (1 - jet[3] * cos(jet[6]))))));
558 REQUIRE(jet[19]
559 == approximately(
560 fp_t{1} / 6
561 * (2 * jet[16] * sin(jet[7]) / (1 - jet[4] * cos(jet[7]))
562 + jet[10]
563 * (2 * jet[13] * cos(jet[7]) * (1 - jet[4] * cos(jet[7]))
564 + sin(jet[7]) * (jet[10] * cos(jet[7]) - jet[4] * 2 * jet[13] * sin(jet[7])))
565 / ((1 - jet[4] * cos(jet[7])) * (1 - jet[4] * cos(jet[7]))))));
566 REQUIRE(jet[20]
567 == approximately(
568 fp_t{1} / 6
569 * (2 * jet[17] * sin(jet[8]) / (1 - jet[5] * cos(jet[8]))
570 + jet[11]
571 * (2 * jet[14] * cos(jet[8]) * (1 - jet[5] * cos(jet[8]))
572 + sin(jet[8]) * (jet[11] * cos(jet[8]) - jet[5] * 2 * jet[14] * sin(jet[8])))
573 / ((1 - jet[5] * cos(jet[8])) * (1 - jet[5] * cos(jet[8]))))));
574
575 REQUIRE(jet[21]
576 == approximately(
577 fp_t{1} / 6
578 * (2 * jet[12] * sin(jet[9]) / (1 - jet[0] * cos(jet[9]))
579 + jet[6]
580 * (2 * jet[15] * cos(jet[9]) * (1 - jet[0] * cos(jet[9]))
581 + sin(jet[9]) * (jet[6] * cos(jet[9]) - jet[0] * 2 * jet[15] * sin(jet[9])))
582 / ((1 - jet[0] * cos(jet[9])) * (1 - jet[0] * cos(jet[9]))))));
583 REQUIRE(jet[22]
584 == approximately(
585 fp_t{1} / 6
586 * (2 * jet[13] * sin(jet[10]) / (1 - jet[1] * cos(jet[10]))
587 + jet[7]
588 * (2 * jet[16] * cos(jet[10]) * (1 - jet[1] * cos(jet[10]))
589 + sin(jet[10]) * (jet[7] * cos(jet[10]) - jet[1] * 2 * jet[16] * sin(jet[10])))
590 / ((1 - jet[1] * cos(jet[10])) * (1 - jet[1] * cos(jet[10]))))));
591 REQUIRE(jet[23]
592 == approximately(
593 fp_t{1} / 6
594 * (2 * jet[14] * sin(jet[11]) / (1 - jet[2] * cos(jet[11]))
595 + jet[8]
596 * (2 * jet[17] * cos(jet[11]) * (1 - jet[2] * cos(jet[11]))
597 + sin(jet[11]) * (jet[8] * cos(jet[11]) - jet[2] * 2 * jet[17] * sin(jet[11])))
598 / ((1 - jet[2] * cos(jet[11])) * (1 - jet[2] * cos(jet[11]))))));
599 }
600
601 // Do the batch/scalar comparison.
602 compare_batch_scalar<fp_t>({kepE(y, expression{number{b}}), kepE(x, expression{number{b}})}, opt_level,
603 high_accuracy, compact_mode);
604
605 // Number-variable tests.
606 {
607 llvm_state s{kw::opt_level = opt_level};
608
609 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, y), kepE(expression{number{c}}, x)}, 1, 1,
610 high_accuracy, compact_mode);
611
612 s.compile();
613
614 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
615
616 std::vector<fp_t> jet{fp_t{.2}, fp_t{.3}};
617 jet.resize(4);
618
619 jptr(jet.data(), nullptr, nullptr);
620
621 REQUIRE(jet[0] == fp_t{.2});
622 REQUIRE(jet[1] == fp_t{.3});
623 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(a, jet[1])));
624 REQUIRE(jet[3] == approximately(bmt_inv_kep_E(c, jet[0])));
625 }
626
627 {
628 llvm_state s{kw::opt_level = opt_level};
629
630 taylor_add_jet<fp_t>(s, "jet", {kepE(par[0], y), kepE(expression{number{c}}, x)}, 1, 1, high_accuracy,
631 compact_mode);
632
633 s.compile();
634
635 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
636
637 std::vector<fp_t> jet{fp_t{.2}, fp_t{.3}};
638 jet.resize(4);
639
640 std::vector<fp_t> pars{a};
641
642 jptr(jet.data(), pars.data(), nullptr);
643
644 REQUIRE(jet[0] == fp_t{.2});
645 REQUIRE(jet[1] == fp_t{.3});
646 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(a, jet[1])));
647 REQUIRE(jet[3] == approximately(bmt_inv_kep_E(c, jet[0])));
648 }
649
650 {
651 llvm_state s{kw::opt_level = opt_level};
652
653 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, y), kepE(expression{number{c}}, x)}, 1, 2,
654 high_accuracy, compact_mode);
655
656 s.compile();
657
658 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
659
660 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.3}, fp_t{.4}};
661 jet.resize(8);
662
663 jptr(jet.data(), nullptr, nullptr);
664
665 REQUIRE(jet[0] == .2);
666 REQUIRE(jet[1] == .5);
667
668 REQUIRE(jet[2] == .3);
669 REQUIRE(jet[3] == .4);
670
671 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(a, jet[2])));
672 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(a, jet[3])));
673
674 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(c, jet[0])));
675 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(c, jet[1])));
676 }
677
678 {
679 llvm_state s{kw::opt_level = opt_level};
680
681 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, y), kepE(par[1], x)}, 1, 2, high_accuracy,
682 compact_mode);
683
684 s.compile();
685
686 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
687
688 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.3}, fp_t{.4}};
689 jet.resize(8);
690
691 std::vector<fp_t> pars{fp_t{0}, fp_t{0}, c, c};
692
693 jptr(jet.data(), pars.data(), nullptr);
694
695 REQUIRE(jet[0] == .2);
696 REQUIRE(jet[1] == .5);
697
698 REQUIRE(jet[2] == .3);
699 REQUIRE(jet[3] == .4);
700
701 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(a, jet[2])));
702 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(a, jet[3])));
703
704 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(c, jet[0])));
705 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(c, jet[1])));
706 }
707
708 {
709 llvm_state s{kw::opt_level = opt_level};
710
711 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, y), kepE(expression{number{c}}, x)}, 2, 1,
712 high_accuracy, compact_mode);
713
714 s.compile();
715
716 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
717
718 std::vector<fp_t> jet{fp_t{.2}, fp_t{.3}};
719 jet.resize(6);
720
721 jptr(jet.data(), nullptr, nullptr);
722
723 REQUIRE(jet[0] == .2);
724 REQUIRE(jet[1] == .3);
725 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(a, jet[1])));
726 REQUIRE(jet[3] == approximately(bmt_inv_kep_E(c, jet[0])));
727 REQUIRE(jet[4] == approximately(fp_t{1} / 2 * jet[3] / (1 - a * cos(jet[2]))));
728 REQUIRE(jet[5] == approximately(fp_t{1} / 2 * jet[2] / (1 - c * cos(jet[3]))));
729 }
730
731 {
732 llvm_state s{kw::opt_level = opt_level};
733
734 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, y), kepE(expression{number{c}}, x)}, 2, 2,
735 high_accuracy, compact_mode);
736
737 s.compile();
738
739 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
740
741 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.3}, fp_t{.4}};
742 jet.resize(12);
743
744 jptr(jet.data(), nullptr, nullptr);
745
746 REQUIRE(jet[0] == .2);
747 REQUIRE(jet[1] == .5);
748
749 REQUIRE(jet[2] == .3);
750 REQUIRE(jet[3] == .4);
751
752 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(a, jet[2])));
753 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(a, jet[3])));
754
755 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(c, jet[0])));
756 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(c, jet[1])));
757
758 REQUIRE(jet[8] == approximately(fp_t{1} / 2 * jet[6] / (1 - a * cos(jet[4]))));
759 REQUIRE(jet[9] == approximately(fp_t{1} / 2 * jet[7] / (1 - a * cos(jet[5]))));
760
761 REQUIRE(jet[10] == approximately(fp_t{1} / 2 * jet[4] / (1 - c * cos(jet[6]))));
762 REQUIRE(jet[11] == approximately(fp_t{1} / 2 * jet[5] / (1 - c * cos(jet[7]))));
763 }
764
765 {
766 llvm_state s{kw::opt_level = opt_level};
767
768 taylor_add_jet<fp_t>(s, "jet", {kepE(expression{number{a}}, y), kepE(expression{number{c}}, x)}, 3, 3,
769 high_accuracy, compact_mode);
770
771 s.compile();
772
773 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
774
775 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.1}, fp_t{.3}, fp_t{.4}, fp_t{.6}};
776 jet.resize(24);
777
778 jptr(jet.data(), nullptr, nullptr);
779
780 REQUIRE(jet[0] == .2);
781 REQUIRE(jet[1] == .5);
782 REQUIRE(jet[2] == .1);
783
784 REQUIRE(jet[3] == .3);
785 REQUIRE(jet[4] == .4);
786 REQUIRE(jet[5] == .6);
787
788 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(a, jet[3])));
789 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(a, jet[4])));
790 REQUIRE(jet[8] == approximately(bmt_inv_kep_E(a, jet[5])));
791
792 REQUIRE(jet[9] == approximately(bmt_inv_kep_E(c, jet[0])));
793 REQUIRE(jet[10] == approximately(bmt_inv_kep_E(c, jet[1])));
794 REQUIRE(jet[11] == approximately(bmt_inv_kep_E(c, jet[2])));
795
796 REQUIRE(jet[12] == approximately(fp_t{1} / 2 * jet[9] / (1 - a * cos(jet[6]))));
797 REQUIRE(jet[13] == approximately(fp_t{1} / 2 * jet[10] / (1 - a * cos(jet[7]))));
798 REQUIRE(jet[14] == approximately(fp_t{1} / 2 * jet[11] / (1 - a * cos(jet[8]))));
799
800 REQUIRE(jet[15] == approximately(fp_t{1} / 2 * jet[6] / (1 - c * cos(jet[9]))));
801 REQUIRE(jet[16] == approximately(fp_t{1} / 2 * jet[7] / (1 - c * cos(jet[10]))));
802 REQUIRE(jet[17] == approximately(fp_t{1} / 2 * jet[8] / (1 - c * cos(jet[11]))));
803
804 REQUIRE(jet[18]
805 == approximately(fp_t{1} / 6
806 * (2 * jet[15] * (1 - a * cos(jet[6])) - a * jet[9] * 2 * jet[12] * sin(jet[6]))
807 / ((1 - a * cos(jet[6])) * (1 - a * cos(jet[6])))));
808 REQUIRE(jet[19]
809 == approximately(fp_t{1} / 6
810 * (2 * jet[16] * (1 - a * cos(jet[7])) - a * jet[10] * 2 * jet[13] * sin(jet[7]))
811 / ((1 - a * cos(jet[7])) * (1 - a * cos(jet[7])))));
812 REQUIRE(jet[20]
813 == approximately(fp_t{1} / 6
814 * (2 * jet[17] * (1 - a * cos(jet[8])) - a * jet[11] * 2 * jet[14] * sin(jet[8]))
815 / ((1 - a * cos(jet[8])) * (1 - a * cos(jet[8])))));
816
817 REQUIRE(jet[21]
818 == approximately(fp_t{1} / 6
819 * (2 * jet[12] * (1 - c * cos(jet[9])) - c * jet[6] * 2 * jet[15] * sin(jet[9]))
820 / ((1 - c * cos(jet[9])) * (1 - c * cos(jet[9])))));
821 REQUIRE(jet[22]
822 == approximately(fp_t{1} / 6
823 * (2 * jet[13] * (1 - c * cos(jet[10])) - c * jet[7] * 2 * jet[16] * sin(jet[10]))
824 / ((1 - c * cos(jet[10])) * (1 - c * cos(jet[10])))));
825 REQUIRE(jet[23]
826 == approximately(fp_t{1} / 6
827 * (2 * jet[14] * (1 - c * cos(jet[11])) - c * jet[8] * 2 * jet[17] * sin(jet[11]))
828 / ((1 - c * cos(jet[11])) * (1 - c * cos(jet[11])))));
829 }
830
831 // Do the batch/scalar comparison.
832 compare_batch_scalar<fp_t>({kepE(expression{number{a}}, y), kepE(expression{number{c}}, x)}, opt_level,
833 high_accuracy, compact_mode);
834
835 // Variable-variable tests.
836 {
837 llvm_state s{kw::opt_level = opt_level};
838
839 taylor_add_jet<fp_t>(s, "jet", {kepE(x, y), kepE(y, x)}, 1, 1, high_accuracy, compact_mode);
840
841 s.compile();
842
843 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
844
845 std::vector<fp_t> jet{fp_t{.2}, fp_t{.3}};
846 jet.resize(4);
847
848 jptr(jet.data(), nullptr, nullptr);
849
850 REQUIRE(jet[0] == fp_t{.2});
851 REQUIRE(jet[1] == fp_t{.3});
852 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(jet[0], jet[1])));
853 REQUIRE(jet[3] == approximately(bmt_inv_kep_E(jet[1], jet[0])));
854 }
855
856 {
857 llvm_state s{kw::opt_level = opt_level};
858
859 taylor_add_jet<fp_t>(s, "jet", {kepE(x, y), kepE(y, x)}, 1, 2, high_accuracy, compact_mode);
860
861 s.compile();
862
863 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
864
865 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.3}, fp_t{.4}};
866 jet.resize(8);
867
868 jptr(jet.data(), nullptr, nullptr);
869
870 REQUIRE(jet[0] == .2);
871 REQUIRE(jet[1] == .5);
872
873 REQUIRE(jet[2] == .3);
874 REQUIRE(jet[3] == .4);
875
876 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(jet[0], jet[2])));
877 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(jet[1], jet[3])));
878
879 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(jet[2], jet[0])));
880 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(jet[3], jet[1])));
881 }
882
883 {
884 llvm_state s{kw::opt_level = opt_level};
885
886 taylor_add_jet<fp_t>(s, "jet", {kepE(x, y), kepE(y, x)}, 2, 1, high_accuracy, compact_mode);
887
888 s.compile();
889
890 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
891
892 std::vector<fp_t> jet{fp_t{.2}, fp_t{.3}};
893 jet.resize(6);
894
895 jptr(jet.data(), nullptr, nullptr);
896
897 REQUIRE(jet[0] == .2);
898 REQUIRE(jet[1] == .3);
899 REQUIRE(jet[2] == approximately(bmt_inv_kep_E(jet[0], jet[1])));
900 REQUIRE(jet[3] == approximately(bmt_inv_kep_E(jet[1], jet[0])));
901 REQUIRE(jet[4]
902 == approximately(fp_t{1} / 2 * (jet[2] * sin(jet[2]) + jet[3]) / (1 - jet[0] * cos(jet[2]))));
903 REQUIRE(jet[5]
904 == approximately(fp_t{1} / 2 * (jet[3] * sin(jet[3]) + jet[2]) / (1 - jet[1] * cos(jet[3]))));
905 }
906
907 {
908 llvm_state s{kw::opt_level = opt_level};
909
910 taylor_add_jet<fp_t>(s, "jet", {kepE(x, y), kepE(y, x)}, 2, 2, high_accuracy, compact_mode);
911
912 s.compile();
913
914 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
915
916 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.3}, fp_t{.4}};
917 jet.resize(12);
918
919 jptr(jet.data(), nullptr, nullptr);
920
921 REQUIRE(jet[0] == .2);
922 REQUIRE(jet[1] == .5);
923
924 REQUIRE(jet[2] == .3);
925 REQUIRE(jet[3] == .4);
926
927 REQUIRE(jet[4] == approximately(bmt_inv_kep_E(jet[0], jet[2])));
928 REQUIRE(jet[5] == approximately(bmt_inv_kep_E(jet[1], jet[3])));
929
930 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(jet[2], jet[0])));
931 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(jet[3], jet[1])));
932
933 REQUIRE(jet[8]
934 == approximately(fp_t{1} / 2 * (jet[4] * sin(jet[4]) + jet[6]) / (1 - jet[0] * cos(jet[4]))));
935 REQUIRE(jet[9]
936 == approximately(fp_t{1} / 2 * (jet[5] * sin(jet[5]) + jet[7]) / (1 - jet[1] * cos(jet[5]))));
937
938 REQUIRE(jet[10]
939 == approximately(fp_t{1} / 2 * (jet[6] * sin(jet[6]) + jet[4]) / (1 - jet[2] * cos(jet[6]))));
940 REQUIRE(jet[11]
941 == approximately(fp_t{1} / 2 * (jet[7] * sin(jet[7]) + jet[5]) / (1 - jet[3] * cos(jet[7]))));
942 }
943
944 {
945 llvm_state s{kw::opt_level = opt_level};
946
947 taylor_add_jet<fp_t>(s, "jet", {kepE(x, y), kepE(y, x)}, 3, 3, high_accuracy, compact_mode);
948
949 s.compile();
950
951 auto jptr = reinterpret_cast<void (*)(fp_t *, const fp_t *, const fp_t *)>(s.jit_lookup("jet"));
952
953 std::vector<fp_t> jet{fp_t{.2}, fp_t{.5}, fp_t{.1}, fp_t{.3}, fp_t{.4}, fp_t{.6}};
954 jet.resize(24);
955
956 jptr(jet.data(), nullptr, nullptr);
957
958 REQUIRE(jet[0] == .2);
959 REQUIRE(jet[1] == .5);
960 REQUIRE(jet[2] == .1);
961
962 REQUIRE(jet[3] == .3);
963 REQUIRE(jet[4] == .4);
964 REQUIRE(jet[5] == .6);
965
966 REQUIRE(jet[6] == approximately(bmt_inv_kep_E(jet[0], jet[3])));
967 REQUIRE(jet[7] == approximately(bmt_inv_kep_E(jet[1], jet[4])));
968 REQUIRE(jet[8] == approximately(bmt_inv_kep_E(jet[2], jet[5])));
969
970 REQUIRE(jet[9] == approximately(bmt_inv_kep_E(jet[3], jet[0])));
971 REQUIRE(jet[10] == approximately(bmt_inv_kep_E(jet[4], jet[1])));
972 REQUIRE(jet[11] == approximately(bmt_inv_kep_E(jet[5], jet[2])));
973
974 REQUIRE(jet[12]
975 == approximately(fp_t{1} / 2 * (jet[6] * sin(jet[6]) + jet[9]) / (1 - jet[0] * cos(jet[6]))));
976 REQUIRE(jet[13]
977 == approximately(fp_t{1} / 2 * (jet[7] * sin(jet[7]) + jet[10]) / (1 - jet[1] * cos(jet[7]))));
978 REQUIRE(jet[14]
979 == approximately(fp_t{1} / 2 * (jet[8] * sin(jet[8]) + jet[11]) / (1 - jet[2] * cos(jet[8]))));
980
981 REQUIRE(jet[15]
982 == approximately(fp_t{1} / 2 * (jet[9] * sin(jet[9]) + jet[6]) / (1 - jet[3] * cos(jet[9]))));
983 REQUIRE(jet[16]
984 == approximately(fp_t{1} / 2 * (jet[10] * sin(jet[10]) + jet[7]) / (1 - jet[4] * cos(jet[10]))));
985 REQUIRE(jet[17]
986 == approximately(fp_t{1} / 2 * (jet[11] * sin(jet[11]) + jet[8]) / (1 - jet[5] * cos(jet[11]))));
987
988 REQUIRE(jet[18]
989 == approximately(fp_t{1} / 6
990 * ((2 * jet[12] * sin(jet[6]) + jet[6] * 2 * jet[12] * cos(jet[6]) + 2 * jet[15])
991 * (1 - jet[0] * cos(jet[6]))
992 - (jet[6] * sin(jet[6]) + jet[9])
993 * (jet[0] * 2 * jet[12] * sin(jet[6]) - jet[6] * cos(jet[6])))
994 / ((1 - jet[0] * cos(jet[6])) * (1 - jet[0] * cos(jet[6])))));
995 REQUIRE(jet[19]
996 == approximately(fp_t{1} / 6
997 * ((2 * jet[13] * sin(jet[7]) + jet[7] * 2 * jet[13] * cos(jet[7]) + 2 * jet[16])
998 * (1 - jet[1] * cos(jet[7]))
999 - (jet[7] * sin(jet[7]) + jet[10])
1000 * (jet[1] * 2 * jet[13] * sin(jet[7]) - jet[7] * cos(jet[7])))
1001 / ((1 - jet[1] * cos(jet[7])) * (1 - jet[1] * cos(jet[7])))));
1002 REQUIRE(jet[20]
1003 == approximately(fp_t{1} / 6
1004 * ((2 * jet[14] * sin(jet[8]) + jet[8] * 2 * jet[14] * cos(jet[8]) + 2 * jet[17])
1005 * (1 - jet[2] * cos(jet[8]))
1006 - (jet[8] * sin(jet[8]) + jet[11])
1007 * (jet[2] * 2 * jet[14] * sin(jet[8]) - jet[8] * cos(jet[8])))
1008 / ((1 - jet[2] * cos(jet[8])) * (1 - jet[2] * cos(jet[8])))));
1009
1010 REQUIRE(jet[21]
1011 == approximately(fp_t{1} / 6
1012 * ((2 * jet[15] * sin(jet[9]) + jet[9] * 2 * jet[15] * cos(jet[9]) + 2 * jet[12])
1013 * (1 - jet[3] * cos(jet[9]))
1014 - (jet[9] * sin(jet[9]) + jet[6])
1015 * (jet[3] * 2 * jet[15] * sin(jet[9]) - jet[9] * cos(jet[9])))
1016 / ((1 - jet[3] * cos(jet[9])) * (1 - jet[3] * cos(jet[9])))));
1017 REQUIRE(
1018 jet[22]
1019 == approximately(fp_t{1} / 6
1020 * ((2 * jet[16] * sin(jet[10]) + jet[10] * 2 * jet[16] * cos(jet[10]) + 2 * jet[13])
1021 * (1 - jet[4] * cos(jet[10]))
1022 - (jet[10] * sin(jet[10]) + jet[7])
1023 * (jet[4] * 2 * jet[16] * sin(jet[10]) - jet[10] * cos(jet[10])))
1024 / ((1 - jet[4] * cos(jet[10])) * (1 - jet[4] * cos(jet[10])))));
1025 REQUIRE(
1026 jet[23]
1027 == approximately(fp_t{1} / 6
1028 * ((2 * jet[17] * sin(jet[11]) + jet[11] * 2 * jet[17] * cos(jet[11]) + 2 * jet[14])
1029 * (1 - jet[5] * cos(jet[11]))
1030 - (jet[11] * sin(jet[11]) + jet[8])
1031 * (jet[5] * 2 * jet[17] * sin(jet[11]) - jet[11] * cos(jet[11])))
1032 / ((1 - jet[5] * cos(jet[11])) * (1 - jet[5] * cos(jet[11])))));
1033 }
1034 };
1035
1036 for (auto cm : {true, false}) {
1037 for (auto f : {false, true}) {
__anon9f23fef20502(auto x) 1038 tuple_for_each(fp_types, [&tester, f, cm](auto x) { tester(x, 0, f, cm); });
__anon9f23fef20602(auto x) 1039 tuple_for_each(fp_types, [&tester, f, cm](auto x) { tester(x, 1, f, cm); });
__anon9f23fef20702(auto x) 1040 tuple_for_each(fp_types, [&tester, f, cm](auto x) { tester(x, 2, f, cm); });
__anon9f23fef20802(auto x) 1041 tuple_for_each(fp_types, [&tester, f, cm](auto x) { tester(x, 3, f, cm); });
1042 }
1043 }
1044 }
1045