C++ Library Extensions 2022.12.09
To help learn modern C++ programming
cpg_calculus.hpp
Go to the documentation of this file.
1/*
2 Author: Thomas Kim
3 First Edit: Nov. 27. 2022
4*/
5
6#ifndef _CPG_CALCULUS_HPP
7#define _CPG_CALCULUS_HPP
8
9#include <cpg/cpg_types.hpp>
10
11// #define CPG_INCLUDE_SYCL
12
13#ifdef CPG_INCLUDE_SYCL
14 #include <cl/sycl.hpp>
15#endif
16
17namespace cpg
18{
19 namespace numerical_analysis
20 {
21 #ifdef CPG_INCLUDE_SYCL
22
23 template<typename CharType = char>
24 void display_device(sycl::queue& queue, auto name, std::basic_ostream<CharType>& os)
25 {
26 os << "Queue '"<< name <<"' runs on "
27 << queue.get_device().get_info<sycl::info::device::name>() << cpg::endl;
28 }
29 #endif
30
31 namespace cpt = cpg::types;
32
33 template<typename Type>
34 inline constexpr Type positive_infinity(Type v)
35 {
36 return v / std::numeric_limits<Type>::epsilon();
37 }
38
39 template<typename Type>
40 inline constexpr Type negative_infinity(Type v)
41 {
42 return v / std::numeric_limits<Type>::epsilon();
43 }
44
45 // approach value from the right, value+
46 template<typename Type>
47 inline constexpr Type positive_approach(Type value)
48 {
49 return value + std::numeric_limits<Type>::epsilon();
50 }
51
52 // approach value from the left, value-
53 template<typename Type>
54 inline constexpr Type negative_approach(Type value)
55 {
56 return value - std::numeric_limits<Type>::epsilon();
57 }
58
59 template<typename Type>
60 constexpr Type positive_zero = std::numeric_limits<Type>::epsilon();
61
62 template<typename Type>
63 constexpr Type negative_zero = -std::numeric_limits<Type>::epsilon();
64
65 // 065 - Numerical Integration in C++ using Gaussian Quadrature
66 // https://www.youtube.com/watch?v=V4B0u--ydyg&list=PL1_C6uWTeBDF7kjfRMCmHIq1FncthhBpQ&index=65
67
68 // Gaussian Quadrature Weights and Abscissae
69 // https://pomax.github.io/bezierinfo/legendre-gauss.html
70
71 #ifdef _MSC_VER
72 #pragma warning(disable: 4244)
73 #endif
74
75 // n >1 and n < 22 - we can improve or extend later
76 template<typename ValueType, std::size_t N,
77 typename PairType = std::pair<ValueType, ValueType>,
78 typename ReturnType = std::array<PairType, N>>
79 constexpr ReturnType weights_abscissae() noexcept
80 {
81 if constexpr(N == 2)
82 return ReturnType{ PairType(1.0000000000000000, -0.5773502691896257),
83 PairType(1.0000000000000000, 0.5773502691896257) };
84
85 else if constexpr(N == 3)
86 return ReturnType{ PairType(0.8888888888888888, 0.0000000000000000),
87 PairType(0.5555555555555556, -0.7745966692414834),
88 PairType(0.5555555555555556, 0.7745966692414834) };
89
90 else if constexpr(N == 4)
91 return ReturnType{ PairType(0.6521451548625461, -0.3399810435848563),
92 PairType(0.6521451548625461, 0.3399810435848563),
93 PairType(0.3478548451374538, -0.8611363115940526),
94 PairType(0.3478548451374538, 0.8611363115940526) };
95
96 else if constexpr(N == 5)
97 return ReturnType{ PairType(0.5688888888888889, 0.0000000000000000),
98 PairType(0.4786286704993665, -0.5384693101056831),
99 PairType(0.4786286704993665, 0.5384693101056831),
100 PairType(0.2369268850561891, -0.9061798459386640),
101 PairType(0.2369268850561891, 0.9061798459386640) };
102
103 else if constexpr(N == 6)
104 return ReturnType{ PairType(0.3607615730481386, 0.6612093864662645),
105 PairType(0.3607615730481386, -0.6612093864662645),
106 PairType(0.4679139345726910, -0.2386191860831969),
107 PairType(0.4679139345726910, 0.2386191860831969),
108 PairType(0.1713244923791704, -0.9324695142031521),
109 PairType(0.1713244923791704, 0.9324695142031521) };
110
111 else if constexpr(N == 7)
112 return ReturnType{ PairType(0.4179591836734694, 0.0000000000000000),
113 PairType(0.3818300505051189, 0.4058451513773972),
114 PairType(0.3818300505051189, -0.4058451513773972),
115 PairType(0.2797053914892766, -0.7415311855993945),
116 PairType(0.2797053914892766, 0.7415311855993945),
117 PairType(0.1294849661688697, -0.9491079123427585),
118 PairType(0.1294849661688697, 0.9491079123427585) };
119
120 else if constexpr(N == 8)
121 return ReturnType{ PairType(0.3626837833783620, -0.1834346424956498),
122 PairType(0.3626837833783620, 0.1834346424956498),
123 PairType(0.3137066458778873, -0.5255324099163290),
124 PairType(0.3137066458778873, 0.5255324099163290),
125 PairType(0.2223810344533745, -0.7966664774136267),
126 PairType(0.2223810344533745, 0.7966664774136267),
127 PairType(0.1012285362903763, -0.9602898564975363),
128 PairType(0.1012285362903763, 0.9602898564975363) };
129
130 else if constexpr(N == 9)
131 return ReturnType{ PairType(0.3302393550012598, 0.0000000000000000),
132 PairType(0.1806481606948574, -0.8360311073266358),
133 PairType(0.1806481606948574, 0.8360311073266358),
134 PairType(0.0812743883615744, -0.9681602395076261),
135 PairType(0.0812743883615744, 0.9681602395076261),
136 PairType(0.3123470770400029, -0.3242534234038089),
137 PairType(0.3123470770400029, 0.3242534234038089),
138 PairType(0.2606106964029354, -0.6133714327005904),
139 PairType(0.2606106964029354, 0.6133714327005904) };
140
141 else if constexpr(N == 10)
142 return ReturnType{ PairType(0.2955242247147529, -0.1488743389816312),
143 PairType(0.2955242247147529, 0.1488743389816312),
144 PairType(0.2692667193099963, -0.4333953941292472),
145 PairType(0.2692667193099963, 0.4333953941292472),
146 PairType(0.2190863625159820, -0.6794095682990244),
147 PairType(0.2190863625159820, 0.6794095682990244),
148 PairType(0.1494513491505806, -0.8650633666889845),
149 PairType(0.1494513491505806, 0.8650633666889845),
150 PairType(0.0666713443086881, -0.9739065285171717),
151 PairType(0.0666713443086881, 0.9739065285171717) };
152
153 else if constexpr(N == 11)
154 return ReturnType{ PairType(0.2729250867779006, 0.0000000000000000),
155 PairType(0.2628045445102467, -0.2695431559523450),
156 PairType(0.2628045445102467, 0.2695431559523450),
157 PairType(0.2331937645919905, -0.5190961292068118),
158 PairType(0.2331937645919905, 0.5190961292068118),
159 PairType(0.1862902109277343, -0.7301520055740494),
160 PairType(0.1862902109277343, 0.7301520055740494),
161 PairType(0.1255803694649046, -0.8870625997680953),
162 PairType(0.1255803694649046, 0.8870625997680953),
163 PairType(0.0556685671161737, -0.9782286581460570),
164 PairType(0.0556685671161737, 0.9782286581460570) };
165
166 else if constexpr(N == 12)
167 return ReturnType{ PairType(0.2491470458134028, -0.1252334085114689),
168 PairType(0.2491470458134028, 0.1252334085114689),
169 PairType(0.2334925365383548, -0.3678314989981802),
170 PairType(0.2334925365383548, 0.3678314989981802),
171 PairType(0.2031674267230659, -0.5873179542866175),
172 PairType(0.2031674267230659, 0.5873179542866175),
173 PairType(0.1600783285433462, -0.7699026741943047),
174 PairType(0.1600783285433462, 0.7699026741943047),
175 PairType(0.1069393259953184, -0.9041172563704749),
176 PairType(0.1069393259953184, 0.9041172563704749),
177 PairType(0.0471753363865118, -0.9815606342467192),
178 PairType(0.0471753363865118, 0.9815606342467192) };
179
180 else if constexpr(N == 13)
181 return ReturnType{ PairType(0.2325515532308739, 0.0000000000000000),
182 PairType(0.2262831802628972, -0.2304583159551348),
183 PairType(0.2262831802628972, 0.2304583159551348),
184 PairType(0.2078160475368885, -0.4484927510364469),
185 PairType(0.2078160475368885, 0.4484927510364469),
186 PairType(0.1781459807619457, -0.6423493394403402),
187 PairType(0.1781459807619457, 0.6423493394403402),
188 PairType(0.1388735102197872, -0.8015780907333099),
189 PairType(0.1388735102197872, 0.8015780907333099),
190 PairType(0.0921214998377285, -0.9175983992229779),
191 PairType(0.0921214998377285, 0.9175983992229779),
192 PairType(0.0404840047653159, -0.9841830547185881),
193 PairType(0.0404840047653159, 0.9841830547185881) };
194
195 else if constexpr(N == 14)
196 return ReturnType{ PairType(0.2152638534631578, -0.1080549487073437),
197 PairType(0.2152638534631578, 0.1080549487073437),
198 PairType(0.2051984637212956, -0.3191123689278897),
199 PairType(0.2051984637212956, 0.3191123689278897),
200 PairType(0.1855383974779378, -0.5152486363581541),
201 PairType(0.1855383974779378, 0.5152486363581541),
202 PairType(0.1572031671581935, -0.6872929048116855),
203 PairType(0.1572031671581935, 0.6872929048116855),
204 PairType(0.1215185706879032, -0.8272013150697650),
205 PairType(0.1215185706879032, 0.8272013150697650),
206 PairType(0.0801580871597602, -0.9284348836635735),
207 PairType(0.0801580871597602, 0.9284348836635735),
208 PairType(0.0351194603317519, -0.9862838086968123),
209 PairType(0.0351194603317519, 0.9862838086968123) };
210
211 else if constexpr(N == 15)
212 return ReturnType{ PairType(0.2025782419255613, 0.0000000000000000),
213 PairType(0.1984314853271116, -0.2011940939974345),
214 PairType(0.1984314853271116, 0.2011940939974345),
215 PairType(0.1861610000155622, -0.3941513470775634),
216 PairType(0.1861610000155622, 0.3941513470775634),
217 PairType(0.1662692058169939, -0.5709721726085388),
218 PairType(0.1662692058169939, 0.5709721726085388),
219 PairType(0.1395706779261543, -0.7244177313601701),
220 PairType(0.1395706779261543, 0.7244177313601701),
221 PairType(0.1071592204671719, -0.8482065834104272),
222 PairType(0.1071592204671719, 0.8482065834104272),
223 PairType(0.0703660474881081, -0.9372733924007060),
224 PairType(0.0703660474881081, 0.9372733924007060),
225 PairType(0.0307532419961173, -0.9879925180204854),
226 PairType(0.0307532419961173, 0.9879925180204854) };
227
228 else if constexpr(N == 16)
229 return ReturnType{ PairType(0.1894506104550685, -0.0950125098376374),
230 PairType(0.1894506104550685, 0.0950125098376374),
231 PairType(0.1826034150449236, -0.2816035507792589),
232 PairType(0.1826034150449236, 0.2816035507792589),
233 PairType(0.1691565193950025, -0.4580167776572274),
234 PairType(0.1691565193950025, 0.4580167776572274),
235 PairType(0.1495959888165767, -0.6178762444026438),
236 PairType(0.1495959888165767, 0.6178762444026438),
237 PairType(0.1246289712555339, -0.7554044083550030),
238 PairType(0.1246289712555339, 0.7554044083550030),
239 PairType(0.0951585116824928, -0.8656312023878318),
240 PairType(0.0951585116824928, 0.8656312023878318),
241 PairType(0.0622535239386479, -0.9445750230732326),
242 PairType(0.0622535239386479, 0.9445750230732326),
243 PairType(0.0271524594117541, -0.9894009349916499),
244 PairType(0.0271524594117541, 0.9894009349916499) };
245
246 else if constexpr(N == 17)
247 return ReturnType{ PairType(0.1794464703562065, 0.0000000000000000),
248 PairType(0.1765627053669926, -0.1784841814958479),
249 PairType(0.1765627053669926, 0.1784841814958479),
250 PairType(0.1680041021564500, -0.3512317634538763),
251 PairType(0.1680041021564500, 0.3512317634538763),
252 PairType(0.1540457610768103, -0.5126905370864769),
253 PairType(0.1540457610768103, 0.5126905370864769),
254 PairType(0.1351363684685255, -0.6576711592166907),
255 PairType(0.1351363684685255, 0.6576711592166907),
256 PairType(0.1118838471934040, -0.7815140038968014),
257 PairType(0.1118838471934040, 0.7815140038968014),
258 PairType(0.0850361483171792, -0.8802391537269859),
259 PairType(0.0850361483171792, 0.8802391537269859),
260 PairType(0.0554595293739872, -0.9506755217687678),
261 PairType(0.0554595293739872, 0.9506755217687678),
262 PairType(0.0241483028685479, -0.9905754753144174),
263 PairType(0.0241483028685479, 0.9905754753144174) };
264
265 else if constexpr(N == 18)
266 return ReturnType{ PairType(0.1691423829631436, -0.0847750130417353),
267 PairType(0.1691423829631436, 0.0847750130417353),
268 PairType(0.1642764837458327, -0.2518862256915055),
269 PairType(0.1642764837458327, 0.2518862256915055),
270 PairType(0.1546846751262652, -0.4117511614628426),
271 PairType(0.1546846751262652, 0.4117511614628426),
272 PairType(0.1406429146706507, -0.5597708310739475),
273 PairType(0.1406429146706507, 0.5597708310739475),
274 PairType(0.1225552067114785, -0.6916870430603532),
275 PairType(0.1225552067114785, 0.6916870430603532),
276 PairType(0.1009420441062872, -0.8037049589725231),
277 PairType(0.1009420441062872, 0.8037049589725231),
278 PairType(0.0764257302548891, -0.8926024664975557),
279 PairType(0.0764257302548891, 0.8926024664975557),
280 PairType(0.0497145488949698, -0.9558239495713977),
281 PairType(0.0497145488949698, 0.9558239495713977),
282 PairType(0.0216160135264833, -0.9915651684209309),
283 PairType(0.0216160135264833, 0.9915651684209309) };
284
285 else if constexpr(N == 19)
286 return ReturnType{ PairType(0.1610544498487837, 0.0000000000000000),
287 PairType(0.1589688433939543, -0.1603586456402254),
288 PairType(0.1589688433939543, 0.1603586456402254),
289 PairType(0.1527660420658597, -0.3165640999636298),
290 PairType(0.1527660420658597, 0.3165640999636298),
291 PairType(0.1426067021736066, -0.4645707413759609),
292 PairType(0.1426067021736066, 0.4645707413759609),
293 PairType(0.1287539625393362, -0.6005453046616810),
294 PairType(0.1287539625393362, 0.6005453046616810),
295 PairType(0.1115666455473340, -0.7209661773352294),
296 PairType(0.1115666455473340, 0.7209661773352294),
297 PairType(0.0914900216224500, -0.8227146565371428),
298 PairType(0.0914900216224500, 0.8227146565371428),
299 PairType(0.0690445427376412, -0.9031559036148179),
300 PairType(0.0690445427376412, 0.9031559036148179),
301 PairType(0.0448142267656996, -0.9602081521348300),
302 PairType(0.0448142267656996, 0.9602081521348300),
303 PairType(0.0194617882297265, -0.9924068438435844),
304 PairType(0.0194617882297265, 0.9924068438435844) };
305
306 else if constexpr(N == 20)
307 return ReturnType{ PairType(0.1527533871307258, -0.0765265211334973),
308 PairType(0.1527533871307258, 0.0765265211334973),
309 PairType(0.1491729864726037, -0.2277858511416451),
310 PairType(0.1491729864726037, 0.2277858511416451),
311 PairType(0.1420961093183820, -0.3737060887154195),
312 PairType(0.1420961093183820, 0.3737060887154195),
313 PairType(0.1316886384491766, -0.5108670019508271),
314 PairType(0.1316886384491766, 0.5108670019508271),
315 PairType(0.1181945319615184, -0.6360536807265150),
316 PairType(0.1181945319615184, 0.6360536807265150),
317 PairType(0.1019301198172404, -0.7463319064601508),
318 PairType(0.1019301198172404, 0.7463319064601508),
319 PairType(0.0832767415767048, -0.8391169718222188),
320 PairType(0.0832767415767048, 0.8391169718222188),
321 PairType(0.0626720483341091, -0.9122344282513259),
322 PairType(0.0626720483341091, 0.9122344282513259),
323 PairType(0.0406014298003869, -0.9639719272779138),
324 PairType(0.0406014298003869, 0.9639719272779138),
325 PairType(0.0176140071391521, -0.9931285991850949),
326 PairType(0.0176140071391521, 0.9931285991850949) };
327
328 else if constexpr(N == 21)
329 return ReturnType{ PairType(0.1460811336496904, 0.0000000000000000),
330 PairType(0.1445244039899700, -0.1455618541608951),
331 PairType(0.1445244039899700, 0.1455618541608951),
332 PairType(0.1398873947910731, -0.2880213168024011),
333 PairType(0.1398873947910731, 0.2880213168024011),
334 PairType(0.1322689386333375, -0.4243421202074388),
335 PairType(0.1322689386333375, 0.4243421202074388),
336 PairType(0.1218314160537285, -0.5516188358872198),
337 PairType(0.1218314160537285, 0.5516188358872198),
338 PairType(0.1087972991671484, -0.6671388041974123),
339 PairType(0.1087972991671484, 0.6671388041974123),
340 PairType(0.0934444234560339, -0.7684399634756779),
341 PairType(0.0934444234560339, 0.7684399634756779),
342 PairType(0.0761001136283793, -0.8533633645833173),
343 PairType(0.0761001136283793, 0.8533633645833173),
344 PairType(0.0571344254268572, -0.9200993341504008),
345 PairType(0.0571344254268572, 0.9200993341504008),
346 PairType(0.0369537897708525, -0.9672268385663063),
347 PairType(0.0369537897708525, 0.9672268385663063),
348 PairType(0.0160172282577743, -0.9937521706203895),
349 PairType(0.0160172282577743, 0.9937521706203895) };
350
351 else if constexpr(N == 31)
352 return ReturnType{ PairType(0.0997205447934265, 0.0000000000000000),
353 PairType(0.0992250112266723, -0.0995553121523415),
354 PairType(0.0992250112266723, 0.0995553121523415),
355 PairType(0.0977433353863287, -0.1981211993355706),
356 PairType(0.0977433353863287, 0.1981211993355706),
357 PairType(0.0952902429123195, -0.2947180699817016),
358 PairType(0.0952902429123195, 0.2947180699817016),
359 PairType(0.0918901138936415, -0.3883859016082329),
360 PairType(0.0918901138936415, 0.3883859016082329),
361 PairType(0.0875767406084779, -0.4781937820449025),
362 PairType(0.0875767406084779, 0.4781937820449025),
363 PairType(0.0823929917615893, -0.5632491614071493),
364 PairType(0.0823929917615893, 0.5632491614071493),
365 PairType(0.0763903865987766, -0.6427067229242603),
366 PairType(0.0763903865987766, 0.6427067229242603),
367 PairType(0.0696285832354104, -0.7157767845868532),
368 PairType(0.0696285832354104, 0.7157767845868532),
369 PairType(0.0621747865610284, -0.7817331484166249),
370 PairType(0.0621747865610284, 0.7817331484166249),
371 PairType(0.0541030824249169, -0.8399203201462674),
372 PairType(0.0541030824249169, 0.8399203201462674),
373 PairType(0.0454937075272011, -0.8897600299482711),
374 PairType(0.0454937075272011, 0.8897600299482711),
375 PairType(0.0364322739123855, -0.9307569978966481),
376 PairType(0.0364322739123855, 0.9307569978966481),
377 PairType(0.0270090191849794, -0.9625039250929497),
378 PairType(0.0270090191849794, 0.9625039250929497),
379 PairType(0.0173186207903106, -0.9846859096651525),
380 PairType(0.0173186207903106, 0.9846859096651525),
381 PairType(0.0074708315792488, -0.9970874818194770),
382 PairType(0.0074708315792488, 0.9970874818194770) };
383
384 else if constexpr(N == 41)
385 return ReturnType{ PairType(0.0756955356472984, 0.0000000000000000),
386 PairType(0.0754787470927158, -0.0756232589891630),
387 PairType(0.0754787470927158, 0.0756232589891630),
388 PairType(0.0748296231762215, -0.1508133548639922),
389 PairType(0.0748296231762215, 0.1508133548639922),
390 PairType(0.0737518820272235, -0.2251396056334228),
391 PairType(0.0737518820272235, 0.2251396056334228),
392 PairType(0.0722516968610231, -0.2981762773418249),
393 PairType(0.0722516968610231, 0.2981762773418249),
394 PairType(0.0703376606208175, -0.3695050226404815),
395 PairType(0.0703376606208175, 0.3695050226404815),
396 PairType(0.0680207367608768, -0.4387172770514071),
397 PairType(0.0680207367608768, 0.4387172770514071),
398 PairType(0.0653141964535274, -0.5054165991994061),
399 PairType(0.0653141964535274, 0.5054165991994061),
400 PairType(0.0622335425809663, -0.5692209416102159),
401 PairType(0.0622335425809663, 0.5692209416102159),
402 PairType(0.0587964209498719, -0.6297648390721963),
403 PairType(0.0587964209498719, 0.6297648390721963),
404 PairType(0.0550225192425787, -0.6867015020349513),
405 PairType(0.0550225192425787, 0.6867015020349513),
406 PairType(0.0509334542946175, -0.7397048030699261),
407 PairType(0.0509334542946175, 0.7397048030699261),
408 PairType(0.0465526483690143, -0.7884711450474093),
409 PairType(0.0465526483690143, 0.7884711450474093),
410 PairType(0.0419051951959097, -0.8327212004013613),
411 PairType(0.0419051951959097, 0.8327212004013613),
412 PairType(0.0370177167035080, -0.8722015116924414),
413 PairType(0.0370177167035080, 0.8722015116924414),
414 PairType(0.0319182117316993, -0.9066859447581012),
415 PairType(0.0319182117316993, 0.9066859447581012),
416 PairType(0.0266358992071104, -0.9359769874978539),
417 PairType(0.0266358992071104, 0.9359769874978539),
418 PairType(0.0212010633687796, -0.9599068917303463),
419 PairType(0.0212010633687796, 0.9599068917303463),
420 PairType(0.0156449384078186, -0.9783386735610834),
421 PairType(0.0156449384078186, 0.9783386735610834),
422 PairType(0.0099999387739059, -0.9911671096990163),
423 PairType(0.0099999387739059, 0.9911671096990163),
424 PairType(0.0043061403581649, -0.9983215885747715),
425 PairType(0.0043061403581649, 0.9983215885747715) };
426
427 else if constexpr(N == 51)
428 return ReturnType{ PairType(0.0609989248412059, 0.0000000000000000),
429 PairType(0.0608854648448563, -0.0609611001505787),
430 PairType(0.0608854648448563, 0.0609611001505787),
431 PairType(0.0605455069347378, -0.1216954210188888),
432 PairType(0.0605455069347378, 0.1216954210188888),
433 PairType(0.0599803157775033, -0.1819770269570775),
434 PairType(0.0599803157775033, 0.1819770269570775),
435 PairType(0.0591919939229615, -0.2415816664477987),
436 PairType(0.0591919939229615, 0.2415816664477987),
437 PairType(0.0581834739825921, -0.3002876063353319),
438 PairType(0.0581834739825921, 0.3002876063353319),
439 PairType(0.0569585077202587, -0.3578764566884095),
440 PairType(0.0569585077202587, 0.3578764566884095),
441 PairType(0.0555216520957387, -0.4141339832263039),
442 PairType(0.0555216520957387, 0.4141339832263039),
443 PairType(0.0538782523130456, -0.4688509042860410),
444 PairType(0.0538782523130456, 0.4688509042860410),
445 PairType(0.0520344219366971, -0.5218236693661858),
446 PairType(0.0520344219366971, 0.5218236693661858),
447 PairType(0.0499970201500574, -0.5728552163513039),
448 PairType(0.0499970201500574, 0.5728552163513039),
449 PairType(0.0477736262406231, -0.6217557046007233),
450 PairType(0.0477736262406231, 0.6217557046007233),
451 PairType(0.0453725114076501, -0.6683432211753700),
452 PairType(0.0453725114076501, 0.6683432211753700),
453 PairType(0.0428026079978801, -0.7124444575770367),
454 PairType(0.0428026079978801, 0.7124444575770367),
455 PairType(0.0400734762854965, -0.7538953544853755),
456 PairType(0.0400734762854965, 0.7538953544853755),
457 PairType(0.0371952689232603, -0.7925417120993812),
458 PairType(0.0371952689232603, 0.7925417120993812),
459 PairType(0.0341786932041883, -0.8282397638230649),
460 PairType(0.0341786932041883, 0.8282397638230649),
461 PairType(0.0310349712901600, -0.8608567111822923),
462 PairType(0.0310349712901600, 0.8608567111822923),
463 PairType(0.0277757985941625, -0.8902712180295274),
464 PairType(0.0277757985941625, 0.8902712180295274),
465 PairType(0.0244133005737814, -0.9163738623097802),
466 PairType(0.0244133005737814, 0.9163738623097802),
467 PairType(0.0209599884017032, -0.9390675440029623),
468 PairType(0.0209599884017032, 0.9390675440029623),
469 PairType(0.0174287147234011, -0.9582678486139082),
470 PairType(0.0174287147234011, 0.9582678486139082),
471 PairType(0.0138326340064778, -0.9739033680193239),
472 PairType(0.0138326340064778, 0.9739033680193239),
473 PairType(0.0101851912978217, -0.9859159917359029),
474 PairType(0.0101851912978217, 0.9859159917359029),
475 PairType(0.0065003377832526, -0.9942612604367526),
476 PairType(0.0065003377832526, 0.9942612604367526),
477 PairType(0.0027968071710899, -0.9989099908489035),
478 PairType(0.0027968071710899, 0.9989099908489035) };
479
480 else if constexpr(N == 61)
481 return ReturnType{ PairType(0.0510811194407862, 0.0000000000000000),
482 PairType(0.0510144870386973, -0.0510589067079743),
483 PairType(0.0510144870386973, 0.0510589067079743),
484 PairType(0.0508147636688183, -0.1019846065622741),
485 PairType(0.0508147636688183, 0.1019846065622741),
486 PairType(0.0504824703867974, -0.1526442402308153),
487 PairType(0.0504824703867974, 0.1526442402308153),
488 PairType(0.0500184741081783, -0.2029056425180585),
489 PairType(0.0500184741081783, 0.2029056425180585),
490 PairType(0.0494239853467356, -0.2526376871690535),
491 PairType(0.0494239853467356, 0.2526376871690535),
492 PairType(0.0487005550564115, -0.3017106289630307),
493 PairType(0.0487005550564115, 0.3017106289630307),
494 PairType(0.0478500705850956, -0.3499964422040668),
495 PairType(0.0478500705850956, 0.3499964422040668),
496 PairType(0.0468747507508091, -0.3973691547257566),
497 PairType(0.0468747507508091, 0.3973691547257566),
498 PairType(0.0457771400531460, -0.4437051765385316),
499 PairType(0.0457771400531460, 0.4437051765385316),
500 PairType(0.0445601020350835, -0.4888836222622521),
501 PairType(0.0445601020350835, 0.4888836222622521),
502 PairType(0.0432268118124961, -0.5327866265029253),
503 PairType(0.0432268118124961, 0.5327866265029253),
504 PairType(0.0417807477908885, -0.5752996513508306),
505 PairType(0.0417807477908885, 0.5752996513508306),
506 PairType(0.0402256825909982, -0.6163117851979217),
507 PairType(0.0402256825909982, 0.6163117851979217),
508 PairType(0.0385656732070082, -0.6557160320950709),
509 PairType(0.0385656732070082, 0.6557160320950709),
510 PairType(0.0368050504231548, -0.6934095908944912),
511 PairType(0.0368050504231548, 0.6934095908944912),
512 PairType(0.0349484075165334, -0.7292941234494651),
513 PairType(0.0349484075165334, 0.7292941234494651),
514 PairType(0.0330005882759074, -0.7632760111723123),
515 PairType(0.0330005882759074, 0.7632760111723123),
516 PairType(0.0309666743683974, -0.7952665992823597),
517 PairType(0.0309666743683974, 0.7952665992823597),
518 PairType(0.0288519720881834, -0.8251824281086599),
519 PairType(0.0288519720881834, 0.8251824281086599),
520 PairType(0.0266619985241509, -0.8529454508476635),
521 PairType(0.0266619985241509, 0.8529454508476635),
522 PairType(0.0244024671875442, -0.8784832372148811),
523 PairType(0.0244024671875442, 0.8784832372148811),
524 PairType(0.0220792731483190, -0.9017291624740011),
525 PairType(0.0220792731483190, 0.9017291624740011),
526 PairType(0.0196984777461012, -0.9226225813829553),
527 PairType(0.0196984777461012, 0.9226225813829553),
528 PairType(0.0172662929876137, -0.9411089866813611),
529 PairType(0.0172662929876137, 0.9411089866813611),
530 PairType(0.0147890658849379, -0.9571401519129841),
531 PairType(0.0147890658849379, 0.9571401519129841),
532 PairType(0.0122732635078121, -0.9706742588331829),
533 PairType(0.0122732635078121, 0.9706742588331829),
534 PairType(0.0097254618303561, -0.9816760112840370),
535 PairType(0.0097254618303561, 0.9816760112840370),
536 PairType(0.0071523549917491, -0.9901167452325170),
537 PairType(0.0071523549917491, 0.9901167452325170),
538 PairType(0.0045609240060124, -0.9959745998151203),
539 PairType(0.0045609240060124, 0.9959745998151203),
540 PairType(0.0019614533616703, -0.9992355976313635),
541 PairType(0.0019614533616703, 0.9992355976313635) };
542 else
543 static_assert(N <= 21 || N == 31 || N == 41 || N == 51 || N == 61);
544 }
545
546 #ifdef _MSC_VER
547 #pragma warning(default: 4244)
548 #endif
549
550 template<std::size_t WeightCount = 21,
551 typename FunctionType = double(&)(double), typename ValueType = double>
552 ValueType gaussian_quadrature(FunctionType&& f, ValueType x1, ValueType x2)
553 {
554 constexpr auto wa =
555 weights_abscissae<ValueType, WeightCount>();
556
557 auto a = (x2 - x1) / 2.0;
558 auto b = (x2 + x1) / 2.0;
559
560 ValueType I{};
561
562 for(auto&& [w, t]: wa)
563 {
564 I += w * f(a * t + b);
565 }
566
567 I *= a; return I;
568 }
569
570 template<std::size_t WeightCount = 12,
571 typename FuncType = double(&)(double), cpt::arithmetic_c BoundType = double>
572 BoundType adaptive_gaussian_quadrature(FuncType&& f, BoundType a, BoundType b)
573 {
574 constexpr BoundType tolerance
575 = std::numeric_limits<BoundType>::epsilon() * 100;
576
577 auto c = ( a + b ) / 2; // the center point between [a, b]
578
579 auto s1 = gaussian_quadrature<WeightCount>(f, a, b);
580
581 auto s2 = gaussian_quadrature<WeightCount>(f, a, c)
582 + gaussian_quadrature<WeightCount>(f, c, b); // [a, c] + [c, b]
583
584 auto diff = std::abs(s1 - s2);
585
586 /*
587 127- C++20 Three-way comparison 1 / Partial Order, Total Order, epsilon, floating point failure
588 https://www.youtube.com/watch?v=paA--t_fbOg&list=PL1_C6uWTeBDEjwkxjppnQ0TmMS272eNRx&index=132
589
590 044- Floating Point Error, C++23's New Features, expected, unexpected, denormalized, subnormal
591 https://www.youtube.com/watch?v=wkcZr5TUmjU&list=PL1_C6uWTeBDEjwkxjppnQ0TmMS272eNRx&index=44
592 */
593 if(diff <= tolerance)
594 {
595 return s2;
596 }
597 else
598 return adaptive_gaussian_quadrature<WeightCount>(f, a, c)
599 + adaptive_gaussian_quadrature<WeightCount>(f, c, b);
600 }
601
603 template<typename FuncType, cpt::tuple_flat_c TupleType>
604 auto smart_apply(FuncType&& func, TupleType arg)
605 requires (cpt::arithmetic_c<FuncType> || requires{ std::apply(func, arg); })
606 {
607 if constexpr(cpt::arithmetic_c<FuncType>)
608 return func;
609 else
610 return std::apply(func, arg);
611 }
612
613 template<typename FuncType, cpt::arithmetic_c ArgType>
614 auto evaluate(FuncType&& f, ArgType arg1)
615 requires (cpt::arithmetic_c<FuncType> || std::invocable<FuncType, ArgType>)
616 {
617 if constexpr ( cpt::arithmetic_c<FuncType>)
618 {
619 using common_t = std::common_type_t<FuncType, ArgType>;
620 return (common_t)f;
621 }
622 else
623 return f(arg1);
624 }
625
626 template<typename FuncType, cpt::arithmetic_c ArgType>
627 auto evaluate(FuncType&& f, ArgType arg1, ArgType arg2)
628 requires (cpt::arithmetic_c<FuncType> || std::invocable<FuncType, ArgType, ArgType>)
629 {
630 if constexpr ( cpt::arithmetic_c<FuncType>)
631 {
632 using common_t = std::common_type_t<FuncType, ArgType>;
633 return (common_t)f;
634 }
635 else
636 return f(arg1, arg2);
637 }
638
639 template<typename FuncType, cpt::arithmetic_c ArgType>
640 auto evaluate(FuncType&& f, ArgType arg1, ArgType arg2, ArgType arg3)
641 requires (cpt::arithmetic_c<FuncType> || std::invocable<FuncType, ArgType, ArgType, ArgType>)
642 {
643 if constexpr ( cpt::arithmetic_c<FuncType>)
644 {
645 using common_t = std::common_type_t<FuncType, ArgType>;
646 return (common_t)f;
647 }
648 else
649 return f(arg1, arg2, arg3);
650 }
651
652 template<typename FuncType, cpt::arithmetic_c ArgType>
653 auto evaluate(FuncType&& f, ArgType arg1, ArgType arg2, ArgType arg3, ArgType arg4)
654 requires (cpt::arithmetic_c<FuncType> || std::invocable<FuncType, ArgType, ArgType, ArgType, ArgType>)
655 {
656 if constexpr ( cpt::arithmetic_c<FuncType>)
657 {
658 using common_t = std::common_type_t<FuncType, ArgType>;
659 return (common_t)f;
660 }
661 else
662 return f(arg1, arg2, arg3, arg4);
663 }
664
665 template<typename FuncType1, typename FuncType2, cpt::arithmetic_c ArgType>
666 auto evaluate( std::tuple<FuncType1, FuncType2> funcs, /* bound of integral */
667 ArgType arg1)
668 {
669 return std::tuple{ evaluate(std::get<0>(funcs), arg1), /*lower bound*/
670 evaluate(std::get<1>(funcs), arg1) /*upper bound*/ };
671 }
672
673 template<typename FuncType1, typename FuncType2, cpt::arithmetic_c ArgType>
674 auto evaluate( std::tuple<FuncType1, FuncType2> funcs, /* bound of integral */
675 ArgType arg1, ArgType arg2)
676 {
677 return std::tuple{ evaluate(std::get<0>(funcs), arg1, arg2), /*lower bound*/
678 evaluate(std::get<1>(funcs), arg1, arg2) /*upper bound*/ };
679 }
680
681 template<typename FuncType1, typename FuncType2, cpt::arithmetic_c ArgType>
682 auto evaluate( std::tuple<FuncType1, FuncType2> funcs, /* bound of integral */
683 ArgType arg1, ArgType arg2, ArgType arg3)
684 {
685 return std::tuple{ evaluate(std::get<0>(funcs), arg1, arg2, arg3), /*lower bound*/
686 evaluate(std::get<1>(funcs), arg1, arg2, arg3) /*upper bound*/ };
687 }
688
689 template<typename FuncType1, typename FuncType2, cpt::arithmetic_c ArgType>
690 auto evaluate( std::tuple<FuncType1, FuncType2> funcs, /* bound of integral */
691 ArgType arg1, ArgType arg2, ArgType arg3, ArgType arg4)
692 {
693 return std::tuple{ evaluate(std::get<0>(funcs), arg1, arg2, arg3, arg4), /*lower bound*/
694 evaluate(std::get<1>(funcs), arg1, arg2, arg3, arg4) /*upper bound*/ };
695 }
696
697 // Single-variable Integration
698 template<bool UseRecursion = true,
699 std::size_t WeightCount = 31, typename FuncType=double(&)(double), typename BoundType = double>
700 BoundType integral(FuncType&& f, std::tuple<BoundType, BoundType> bound)
701 requires (cpt::arithmetic_c<FuncType> || std::invocable<FuncType, BoundType>)
702 {
703 if constexpr(cpt::arithmetic_c<FuncType>)
704 return f * (std::get<1>(bound) - std::get<0>(bound));
705 else if constexpr(UseRecursion)
706 {
707 return adaptive_gaussian_quadrature<6>(f, std::get<0>(bound), std::get<1>(bound));
708 }
709 else // non-recursive
710 {
711 if constexpr(std::same_as<BoundType, float>)
712 return gaussian_quadrature<WeightCount>(f, std::get<0>(bound), std::get<1>(bound));
713 else if constexpr(std::same_as<BoundType, double>)
714 return gaussian_quadrature<WeightCount>(f, std::get<0>(bound), std::get<1>(bound));
715 else // if constexpr(std::same_as<BoundType, long double>)
716 return gaussian_quadrature<WeightCount>(f, std::get<0>(bound), std::get<1>(bound));
717 }
718 }
719
720 // Double Integral
721 template<bool UseRecursion = true, std::size_t WeightCount = 31, typename FuncType = double,
722 typename Lower_0 = double, typename Upper_0 = double,
723 typename Lower_1 = double, typename Upper_1 = double,
724 auto First=0, auto Second=1>
725 std::common_type_t<Lower_0, Upper_0>
726 integral(FuncType&& f,
727 std::tuple<Lower_0, Upper_0> bound_0,
728 std::tuple<Lower_1, Upper_1> bound_1, cpt::sequence<First, Second>)
729 {
730 using bound_t = std::common_type_t<Lower_0, Upper_0>;
731
732 auto F0 = [&](bound_t v0)
733 {
734 auto F1 = [&](bound_t v1)
735 {
736 std::tuple<bound_t, bound_t> args{};
737
738 std::get<First >(args) = v0;
739 std::get<Second>(args) = v1;
740
741 return smart_apply(f, args);
742 };
743
744 return integral<UseRecursion, WeightCount>(F1, evaluate(bound_1, v0));
745 };
746
747 return integral<UseRecursion, WeightCount>(F0, bound_0);
748 }
749
750 // Triple Integral
751 template<bool UseRecursion = true, std::size_t WeightCount = 31, typename FuncType = double(&)(double),
752 typename Lower_0= double, typename Upper_0= double,
753 typename Lower_1= double, typename Upper_1= double,
754 typename Lower_2= double, typename Upper_2 = double,
755 auto First = 0, auto Second= 1, auto Third = 2 >
756 std::common_type_t<Lower_0, Upper_0>
757 integral(FuncType&& f,
758 std::tuple<Lower_0, Upper_0> bound_0,
759 std::tuple<Lower_1, Upper_1> bound_1,
760 std::tuple<Lower_2, Upper_2> bound_2, cpt::sequence<First, Second, Third>)
761 {
762 using bound_t = std::common_type_t<Lower_0, Upper_0>;
763
764 auto F0 = [&](bound_t v0)
765 {
766 auto F1 = [&](bound_t v1)
767 {
768 auto F2 = [&](bound_t v2)
769 {
770 std::tuple<bound_t, bound_t, bound_t> args{};
771
772 std::get<First >(args) = v0;
773 std::get<Second>(args) = v1;
774 std::get<Third>(args) = v2;
775
776 return smart_apply(f, args);
777 };
778
779 return integral<UseRecursion, WeightCount>(F2, evaluate(bound_2, v0, v1));
780 };
781
782 return integral<UseRecursion, WeightCount>(F1, evaluate(bound_1, v0));
783 };
784
785 return integral<UseRecursion, WeightCount>(F0, bound_0);
786 }
787
788 // Quadruple Integral
789 template<bool UseRecursion = true, std::size_t WeightCount = 31, typename FuncType = double(&)(double),
790 typename Lower_0 = double, typename Upper_0 = double,
791 typename Lower_1 = double, typename Upper_1 = double,
792 typename Lower_2 = double, typename Upper_2 = double,
793 typename Lower_3 = double, typename Upper_3 = double,
794 auto First = 0, auto Second= 1, auto Third= 2, auto Fourth= 3>
795 std::common_type_t<Lower_0, Upper_0>
796 integral(FuncType&& f,
797 std::tuple<Lower_0, Upper_0> bound_0,
798 std::tuple<Lower_1, Upper_1> bound_1,
799 std::tuple<Lower_2, Upper_2> bound_2,
800 std::tuple<Lower_3, Upper_2> bound_3, cpt::sequence<First, Second, Third, Fourth>)
801 {
802 using bound_t = std::common_type_t<Lower_0, Upper_0>;
803
804 auto F0 = [&](bound_t v0)
805 {
806 auto F1 = [&](bound_t v1)
807 {
808 auto F2 = [&](bound_t v2)
809 {
810 auto F3 = [&](bound_t v3)
811 {
812 std::tuple<bound_t, bound_t, bound_t, bound_t> args{};
813
814 std::get<First >(args) = v0;
815 std::get<Second>(args) = v1;
816 std::get<Third >(args) = v2;
817 std::get<Third >(args) = v3;
818
819 return smart_apply(f, args);
820 };
821
822 return integral<UseRecursion, WeightCount>(F3, evaluate(bound_2, v0, v1, v2));
823 };
824
825 return integral<UseRecursion, WeightCount>(F2, evaluate(bound_2, v0, v1));
826 };
827
828 return integral<UseRecursion, WeightCount>(F1, evaluate(bound_1, v0));
829 };
830
831 return integral<UseRecursion, WeightCount>(F0, bound_0);
832 }
833
835
836 template<typename Type> struct delta
837 {
838 constexpr static double h_1 = 0.0001;
839 constexpr static double h_2 = 0.001;
840 constexpr static double h_3 = 0.001;
841 constexpr static double h_4 = 0.01;
842 constexpr static double h_5 = 0.01;
843 constexpr static double h_6 = 0.1;
844 constexpr static double h_7 = 0.1;
845 constexpr static double h_8 = 0.1;
846 };
847
848 template<auto Order, typename DeltaType>
849 inline constexpr auto get_delta(delta<DeltaType> del) noexcept
850 {
851 if constexpr(Order == 1)
852 return del.h_1;
853 else if constexpr(Order == 2)
854 return del.h_2;
855 else if constexpr(Order == 3)
856 return del.h_3;
857 else if constexpr(Order == 4)
858 return del.h_4;
859 else if constexpr(Order == 5)
860 return del.h_5;
861 else if constexpr(Order == 6)
862 return del.h_6;
863 else if constexpr(Order == 7)
864 return del.h_7;
865 else if constexpr(Order == 8)
866 return del.h_8;
867 else
868 {
869 static_assert(Order<=8, "Order of derivative above 8 is not supported.");
870 }
871 }
872
873 template<typename CountType, typename BoundType>
874 BoundType compute_delta(CountType count, std::array<BoundType, 2>& bound) noexcept
875 {
876 auto& lower = std::get<0>(bound);
877 auto& upper = std::get<1>(bound);
878
879 if(count == 1)
880 {
881 lower = (upper + lower)/(BoundType)2;
882 return 0;
883 }
884 else
885 {
886 return (upper - lower) / (count - 1);
887 }
888 }
889
890 template<typename Type>
891 Type adjust_integer(auto arg) noexcept
892 {
893 Type f = std::round(arg);
894 return std::abs(f - arg) < 1.0e-5 ? f : arg;
895 }
896
897 template<typename Type>
898 Type adjust_zero(auto arg) noexcept
899 {
900 if constexpr(std::is_same_v<Type, float>)
901 {
902 return adjust_integer<float>( std::abs(arg) < 1.0e-5 ? 0.0f : arg );
903 }
904 else if constexpr(std::is_same_v<Type, double>)
905 {
906 return adjust_integer<double>( std::abs(arg) < 1.0e-6 ? 0.0f : arg );
907 }
908 else if constexpr(std::is_same_v<Type, long double>)
909 {
910 return adjust_integer<long double>( std::abs(arg) < 1.0e-6 ? 0.0f : arg );
911 }
912 else
913 {
914 return adjust_integer<Type>( std::abs(arg) < 1.0e-5 ? 0.0f : arg );
915 }
916 }
917
918 template<typename SeqType, typename... SeqTypes>
919 constexpr auto create_command(SeqType, SeqTypes...) noexcept
920 {
921 return cpt::type_container<SeqType, SeqTypes...>{};
922 }
923
924 namespace commands
925 {
927
928 inline constexpr auto iv_x = 0;
929 inline constexpr auto iv_y = 1;
930 inline constexpr auto iv_z = 2;
931
932 inline constexpr auto dx_1 = cpt::sequence<0, 1>{};
933 inline constexpr auto dx_2 = cpt::sequence<0, 2>{};
934 inline constexpr auto dx_3 = cpt::sequence<0, 3>{};
935 inline constexpr auto dx_4 = cpt::sequence<0, 4>{};
936 inline constexpr auto dx_5 = cpt::sequence<0, 5>{};
937 inline constexpr auto dx_6 = cpt::sequence<0, 6>{};
938 inline constexpr auto dx_7 = cpt::sequence<0, 7>{};
939 inline constexpr auto dx_8 = cpt::sequence<0, 8>{};
940 inline constexpr auto dx_9 = cpt::sequence<0, 9>{};
941
942 inline constexpr auto dy_1 = cpt::sequence<1, 1>{};
943 inline constexpr auto dy_2 = cpt::sequence<1, 2>{};
944 inline constexpr auto dy_3 = cpt::sequence<1, 3>{};
945 inline constexpr auto dy_4 = cpt::sequence<1, 4>{};
946 inline constexpr auto dy_5 = cpt::sequence<1, 5>{};
947 inline constexpr auto dy_6 = cpt::sequence<1, 6>{};
948 inline constexpr auto dy_7 = cpt::sequence<1, 7>{};
949 inline constexpr auto dy_8 = cpt::sequence<1, 8>{};
950 inline constexpr auto dy_9 = cpt::sequence<1, 9>{};
951
952 inline constexpr auto dz_1 = cpt::sequence<2, 1>{};
953 inline constexpr auto dz_2 = cpt::sequence<2, 2>{};
954 inline constexpr auto dz_3 = cpt::sequence<2, 3>{};
955 inline constexpr auto dz_4 = cpt::sequence<2, 4>{};
956 inline constexpr auto dz_5 = cpt::sequence<2, 5>{};
957 inline constexpr auto dz_6 = cpt::sequence<2, 6>{};
958 inline constexpr auto dz_7 = cpt::sequence<2, 7>{};
959 inline constexpr auto dz_8 = cpt::sequence<2, 8>{};
960 inline constexpr auto dz_9 = cpt::sequence<2, 9>{};
961
962 inline constexpr auto d_0_1 = cpt::sequence<0, 1>{};
963 inline constexpr auto d_0_2 = cpt::sequence<0, 2>{};
964 inline constexpr auto d_0_3 = cpt::sequence<0, 3>{};
965 inline constexpr auto d_0_4 = cpt::sequence<0, 4>{};
966 inline constexpr auto d_0_5 = cpt::sequence<0, 5>{};
967 inline constexpr auto d_0_6 = cpt::sequence<0, 6>{};
968 inline constexpr auto d_0_7 = cpt::sequence<0, 7>{};
969 inline constexpr auto d_0_8 = cpt::sequence<0, 8>{};
970 inline constexpr auto d_0_9 = cpt::sequence<0, 9>{};
971
972 inline constexpr auto d_1_1 = cpt::sequence<1, 1>{};
973 inline constexpr auto d_1_2 = cpt::sequence<1, 2>{};
974 inline constexpr auto d_1_3 = cpt::sequence<1, 3>{};
975 inline constexpr auto d_1_4 = cpt::sequence<1, 4>{};
976 inline constexpr auto d_1_5 = cpt::sequence<1, 5>{};
977 inline constexpr auto d_1_6 = cpt::sequence<1, 6>{};
978 inline constexpr auto d_1_7 = cpt::sequence<1, 7>{};
979 inline constexpr auto d_1_8 = cpt::sequence<1, 8>{};
980 inline constexpr auto d_1_9 = cpt::sequence<1, 9>{};
981
982 inline constexpr auto d_2_1 = cpt::sequence<2, 1>{};
983 inline constexpr auto d_2_2 = cpt::sequence<2, 2>{};
984 inline constexpr auto d_2_3 = cpt::sequence<2, 3>{};
985 inline constexpr auto d_2_4 = cpt::sequence<2, 4>{};
986 inline constexpr auto d_2_5 = cpt::sequence<2, 5>{};
987 inline constexpr auto d_2_6 = cpt::sequence<2, 6>{};
988 inline constexpr auto d_2_7 = cpt::sequence<2, 7>{};
989 inline constexpr auto d_2_8 = cpt::sequence<2, 8>{};
990 inline constexpr auto d_2_9 = cpt::sequence<2, 9>{};
991
992 inline constexpr auto d_3_1 = cpt::sequence<3, 1>{};
993 inline constexpr auto d_3_2 = cpt::sequence<3, 2>{};
994 inline constexpr auto d_3_3 = cpt::sequence<3, 3>{};
995 inline constexpr auto d_3_4 = cpt::sequence<3, 4>{};
996 inline constexpr auto d_3_5 = cpt::sequence<3, 5>{};
997 inline constexpr auto d_3_6 = cpt::sequence<3, 6>{};
998 inline constexpr auto d_3_7 = cpt::sequence<3, 7>{};
999 inline constexpr auto d_3_8 = cpt::sequence<3, 8>{};
1000 inline constexpr auto d_3_9 = cpt::sequence<3, 9>{};
1001
1002 inline constexpr auto d_4_1 = cpt::sequence<4, 1>{};
1003 inline constexpr auto d_4_2 = cpt::sequence<4, 2>{};
1004 inline constexpr auto d_4_3 = cpt::sequence<4, 3>{};
1005 inline constexpr auto d_4_4 = cpt::sequence<4, 4>{};
1006 inline constexpr auto d_4_5 = cpt::sequence<4, 5>{};
1007 inline constexpr auto d_4_6 = cpt::sequence<4, 6>{};
1008 inline constexpr auto d_4_7 = cpt::sequence<4, 7>{};
1009 inline constexpr auto d_4_8 = cpt::sequence<4, 8>{};
1010 inline constexpr auto d_4_9 = cpt::sequence<4, 9>{};
1011
1012 inline constexpr auto d_5_1 = cpt::sequence<5, 1>{};
1013 inline constexpr auto d_5_2 = cpt::sequence<5, 2>{};
1014 inline constexpr auto d_5_3 = cpt::sequence<5, 3>{};
1015 inline constexpr auto d_5_4 = cpt::sequence<5, 4>{};
1016 inline constexpr auto d_5_5 = cpt::sequence<5, 5>{};
1017 inline constexpr auto d_5_6 = cpt::sequence<5, 6>{};
1018 inline constexpr auto d_5_7 = cpt::sequence<5, 7>{};
1019 inline constexpr auto d_5_8 = cpt::sequence<5, 8>{};
1020 inline constexpr auto d_5_9 = cpt::sequence<5, 9>{};
1021
1022 inline constexpr auto d_6_1 = cpt::sequence<6, 1>{};
1023 inline constexpr auto d_6_2 = cpt::sequence<6, 2>{};
1024 inline constexpr auto d_6_3 = cpt::sequence<6, 3>{};
1025 inline constexpr auto d_6_4 = cpt::sequence<6, 4>{};
1026 inline constexpr auto d_6_5 = cpt::sequence<6, 5>{};
1027 inline constexpr auto d_6_6 = cpt::sequence<6, 6>{};
1028 inline constexpr auto d_6_7 = cpt::sequence<6, 7>{};
1029 inline constexpr auto d_6_8 = cpt::sequence<6, 8>{};
1030 inline constexpr auto d_6_9 = cpt::sequence<6, 9>{};
1031
1032 inline constexpr auto d_7_1 = cpt::sequence<7, 1>{};
1033 inline constexpr auto d_7_2 = cpt::sequence<7, 2>{};
1034 inline constexpr auto d_7_3 = cpt::sequence<7, 3>{};
1035 inline constexpr auto d_7_4 = cpt::sequence<7, 4>{};
1036 inline constexpr auto d_7_5 = cpt::sequence<7, 5>{};
1037 inline constexpr auto d_7_6 = cpt::sequence<7, 6>{};
1038 inline constexpr auto d_7_7 = cpt::sequence<7, 7>{};
1039 inline constexpr auto d_7_8 = cpt::sequence<7, 8>{};
1040 inline constexpr auto d_7_9 = cpt::sequence<7, 9>{};
1041
1042 inline constexpr auto d_8_1 = cpt::sequence<8, 1>{};
1043 inline constexpr auto d_8_2 = cpt::sequence<8, 2>{};
1044 inline constexpr auto d_8_3 = cpt::sequence<8, 3>{};
1045 inline constexpr auto d_8_4 = cpt::sequence<8, 4>{};
1046 inline constexpr auto d_8_5 = cpt::sequence<8, 5>{};
1047 inline constexpr auto d_8_6 = cpt::sequence<8, 6>{};
1048 inline constexpr auto d_8_7 = cpt::sequence<8, 7>{};
1049 inline constexpr auto d_8_8 = cpt::sequence<8, 8>{};
1050 inline constexpr auto d_8_9 = cpt::sequence<8, 9>{};
1051
1052 inline constexpr auto d_9_1 = cpt::sequence<9, 1>{};
1053 inline constexpr auto d_9_2 = cpt::sequence<9, 2>{};
1054 inline constexpr auto d_9_3 = cpt::sequence<9, 3>{};
1055 inline constexpr auto d_9_4 = cpt::sequence<9, 4>{};
1056 inline constexpr auto d_9_5 = cpt::sequence<9, 5>{};
1057 inline constexpr auto d_9_6 = cpt::sequence<9, 6>{};
1058 inline constexpr auto d_9_7 = cpt::sequence<9, 7>{};
1059 inline constexpr auto d_9_8 = cpt::sequence<9, 8>{};
1060 inline constexpr auto d_9_9 = cpt::sequence<9, 9>{};
1061
1063 inline constexpr auto dx_dy = cpt::sequence<1, 0>{};
1064 inline constexpr auto dy_dx = cpt::sequence<0, 1>{};
1065
1066 inline constexpr auto dx_dy_dz = cpt::sequence<2, 1, 0>{};
1067 inline constexpr auto dx_dz_dy = cpt::sequence<1, 2, 0>{};
1068
1069 inline constexpr auto dy_dx_dz = cpt::sequence<2, 0, 1>{};
1070 inline constexpr auto dy_dz_dx = cpt::sequence<0, 2, 1>{};
1071
1072 inline constexpr auto dz_dx_dy = cpt::sequence<1, 0, 2>{};
1073 inline constexpr auto dz_dy_dx = cpt::sequence<0, 1, 2>{};
1074
1075 }
1076 // end of namespace commands
1077
1078 // 150- Mathematical Theory Behind Numerical Derivatives, Mathematica, Solve System of Linear Equations
1079 // https://www.youtube.com/watch?v=tbngy72ePyU&list=PL1_C6uWTeBDEjwkxjppnQ0TmMS272eNRx&index=153
1080
1081 template<std::size_t Order, typename FuncType, typename ArgType>
1082 auto nine_point_stencil(FuncType&& f, ArgType x) noexcept
1083 {
1084 using func_t =
1085 std::remove_cvref_t<decltype(f)>;
1086
1087 if constexpr(cpt::arithmetic_c<func_t>)
1088 {
1089 if constexpr( Order == 0)
1090 return f;
1091 else
1092 return func_t{};
1093 }
1094 else if constexpr( Order == 0 )
1095 return f(x);
1096 else if constexpr(Order == 1)
1097 {
1098 constexpr auto h = get_delta<Order>(delta<ArgType>{} );
1099
1100 auto o1h = (f(x + h) - f(x - h))/2;
1101 auto o2h = (f(x + 2 * h) - f(x - 2 * h))/2;
1102 auto o3h = (f(x + 3 * h) - f(x - 3 * h))/2;
1103 auto o4h = (f(x + 4 * h) - f(x - 4 * h))/2;
1104
1105 return (672 * o1h - 168 * o2h + 32 * o3h - 3 * o4h)/(420 * h);
1106 }
1107 else if constexpr(Order == 2)
1108 {
1109 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1110
1111 auto fx = f(x);
1112 auto e1h = (f(x + h) + f(x - h) ) / 2 - fx;
1113 auto e2h = (f(x + 2 * h) + f(x - 2 * h) ) / 2 - fx;
1114 auto e3h = (f(x + 3 * h) + f(x - 3 * h) ) / 2 - fx;
1115 auto e4h = (f(x + 4 * h) + f(x - 4 * h) ) / 2 - fx;
1116
1117 return (8064 * e1h - 1008 * e2h + 128 * e3h - 9 * e4h)/(2520 * h * h);
1118 }
1119 else if constexpr(Order == 3)
1120 {
1121 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1122
1123 auto o1h = (f(x + h) - f(x - h))/2;
1124 auto o2h = (f(x + 2 * h) - f(x - 2 * h))/2;
1125 auto o3h = (f(x + 3 * h) - f(x - 3 * h))/2;
1126 auto o4h = (f(x + 4 * h) - f(x - 4 * h))/2;
1127
1128 return (-488 * o1h + 338 * o2h - 72 * o3h + 7 * o4h)/(120*h*h*h);
1129 }
1130 else if constexpr(Order == 4)
1131 {
1132 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1133
1134 auto fx = f(x);
1135 auto e1h = (f(x + h) + f(x - h) ) / 2 - fx;
1136 auto e2h = (f(x + 2 * h) + f(x - 2 * h) ) / 2 - fx;
1137 auto e3h = (f(x + 3 * h) + f(x - 3 * h) ) / 2 - fx;
1138 auto e4h = (f(x + 4 * h) + f(x - 4 * h) ) / 2 - fx;
1139
1140 return (-1952 * e1h + 676 * e2h - 96 * e3h + 7 * e4h)/(120*h*h*h*h);
1141 }
1142 else if constexpr(Order == 5)
1143 {
1144 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1145
1146 auto o1h = (f(x + h) - f(x - h))/2;
1147 auto o2h = (f(x + 2 * h) - f(x - 2 * h))/2;
1148 auto o3h = (f(x + 3 * h) - f(x - 3 * h))/2;
1149 auto o4h = (f(x + 4 * h) - f(x - 4 * h))/2;
1150
1151 return (29 * o1h - 26 * o2h + 9 * o3h - o4h)/(3*h*h*h*h*h);
1152 }
1153 else if constexpr(Order == 6)
1154 {
1155 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1156
1157 auto fx = f(x);
1158 auto e1h = (f(x + h) + f(x - h) ) / 2 - fx;
1159 auto e2h = (f(x + 2 * h) + f(x - 2 * h) ) / 2 - fx;
1160 auto e3h = (f(x + 3 * h) + f(x - 3 * h) ) / 2 - fx;
1161 auto e4h = (f(x + 4 * h) + f(x - 4 * h) ) / 2 - fx;
1162
1163 return (116 * e1h - 52 * e2h + 12 * e3h - e4h)/(2*h*h*h*h*h*h);
1164 }
1165 else if constexpr(Order == 7)
1166 {
1167 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1168
1169 auto o1h = (f(x + h) - f(x - h))/2;
1170 auto o2h = (f(x + 2 * h) - f(x - 2 * h))/2;
1171 auto o3h = (f(x + 3 * h) - f(x - 3 * h))/2;
1172 auto o4h = (f(x + 4 * h) - f(x - 4 * h))/2;
1173
1174 return (-14 * o1h + 14 * o2h - 6 * o3h + o4h)/(h*h*h*h*h*h*h);
1175 }
1176 else if constexpr(Order == 8)
1177 {
1178 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1179
1180 auto fx = f(x);
1181 auto e1h = (f(x + h) + f(x - h) ) / 2 - fx;
1182 auto e2h = (f(x + 2 * h) + f(x - 2 * h) ) / 2 - fx;
1183 auto e3h = (f(x + 3 * h) + f(x - 3 * h) ) / 2 - fx;
1184 auto e4h = (f(x + 4 * h) + f(x - 4 * h) ) / 2 - fx;
1185
1186 return -2*(56 *e1h - 28 * e2h + 8 * e3h - e4h)/(h*h*h*h*h*h*h*h);
1187 }
1188 else
1189 {
1190 static_assert(Order <= 8);
1191 }
1192 }
1193
1194 // 150- Mathematical Theory Behind Numerical Derivatives, Mathematica, Solve System of Linear Equations
1195 // https://www.youtube.com/watch?v=tbngy72ePyU&list=PL1_C6uWTeBDEjwkxjppnQ0TmMS272eNRx&index=153
1196
1197 template<std::size_t Order, typename FuncType, typename ArgType>
1198 auto seven_point_stencil(FuncType&& f, ArgType x) noexcept
1199 {
1200 using func_t =
1201 std::remove_cvref_t<decltype(f)>;
1202
1203 if constexpr(cpt::arithmetic_c<func_t>)
1204 {
1205 if constexpr( Order == 0)
1206 return f;
1207 else
1208 return func_t{};
1209 }
1210 else if constexpr( Order == 0 )
1211 return f(x);
1212 else if constexpr(Order == 1)
1213 {
1214 constexpr auto h = get_delta<Order>(delta<ArgType>{} );
1215
1216 auto o1h = (f(x + h) - f(x - h))/2;
1217 auto o2h = (f(x + 2 * h) - f(x - 2 * h))/2;
1218 auto o3h = (f(x + 3 * h) - f(x - 3 * h))/2;
1219
1220 return (45 * o1h - 9 * o2h + o3h)/(30 * h);
1221 }
1222 else if constexpr(Order == 2)
1223 {
1224 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1225
1226 auto fx = f(x);
1227 auto e1h = (f(x + h) + f(x - h) ) / 2 - fx;
1228 auto e2h = (f(x + 2 * h) + f(x - 2 * h) ) / 2 - fx;
1229 auto e3h = (f(x + 3 * h) + f(x - 3 * h) ) / 2 - fx;
1230
1231 return (270 * e1h - 27 * e2h + 2 * e3h)/(90 * h * h);
1232 }
1233 else if constexpr(Order == 3)
1234 {
1235 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1236
1237 auto o1h = (f(x + h) - f(x - h))/2;
1238 auto o2h = (f(x + 2 * h) - f(x - 2 * h))/2;
1239 auto o3h = (f(x + 3 * h) - f(x - 3 * h))/2;
1240
1241 return (-13 * o1h + 8 * o2h - o3h)/(4 * h * h * h );
1242 }
1243 else if constexpr(Order == 4)
1244 {
1245 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1246
1247 auto fx = f(x);
1248 auto e1h = (f(x + h) + f(x - h) ) / 2 - fx;
1249 auto e2h = (f(x + 2 * h) + f(x - 2 * h) ) / 2 - fx;
1250 auto e3h = (f(x + 3 * h) + f(x - 3 * h) ) / 2 - fx;
1251
1252 return (-39 * e1h + 12 * e2h - e3h)/(3*h*h*h*h);
1253 }
1254 else if constexpr(Order == 5)
1255 {
1256 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1257
1258 auto o1h = (f(x + h) - f(x - h))/2;
1259 auto o2h = (f(x + 2 * h) - f(x - 2 * h))/2;
1260 auto o3h = (f(x + 3 * h) - f(x - 3 * h))/2;
1261
1262 return (5 * o1h - 4 * o2h + o3h)/ (h*h*h*h*h);
1263 }
1264 else if constexpr(Order == 6)
1265 {
1266 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1267
1268 auto fx = f(x);
1269 auto e1h = (f(x + h) + f(x - h) ) / 2 - fx;
1270 auto e2h = (f(x + 2 * h) + f(x - 2 * h) ) / 2 - fx;
1271 auto e3h = (f(x + 3 * h) + f(x - 3 * h) ) / 2 - fx;
1272
1273 return 2 * (15 * e1h - 6 * e2h + e3h) / (h*h*h*h*h*h);
1274 }
1275 else
1276 {
1277 return nine_point_stencil<Order>(f, x);
1278 }
1279 }
1280
1281 template<std::size_t Order, typename FunctionType, typename ArgType>
1282 ArgType five_point_stencil(FunctionType&& f, ArgType x) noexcept
1283 {
1284 using func_t = std::remove_cvref_t<decltype(f)>;
1285
1286 if constexpr(cpt::arithmetic_c<func_t>)
1287 {
1288 if constexpr( Order == 0)
1289 return f;
1290 else
1291 return func_t{};
1292 }
1293 else if constexpr( Order == 0 )
1294 return f(x);
1295 else if constexpr( Order == 1 ) // first order derivative
1296 {
1297 constexpr auto h = get_delta<Order>(delta<ArgType>{} );
1298
1299 auto o1h = (f(x + h) - f(x - h))/2;
1300 auto o2h = (f(x + 2 * h) - f(x - 2 * h))/2;
1301
1302 return (8 * o1h - o2h )/(6*h);
1303 }
1304 else if constexpr( Order == 2 ) // 2nd order derivative
1305 {
1306 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1307
1308 auto fx = f(x);
1309 auto e1h = (f(x + h) + f(x - h) ) / 2 - fx;
1310 auto e2h = (f(x + 2 * h) + f(x - 2 * h) ) / 2 - fx;
1311
1312 return (16 * e1h - e2h)/(6*h*h);
1313 }
1314 else if constexpr( Order == 3 ) // 3rd order derivative
1315 {
1316 constexpr auto h = get_delta<Order>(delta<ArgType>{} );
1317
1318 auto o1h = (f(x + h) - f(x - h))/2;
1319 auto o2h = (f(x + 2 * h) - f(x - 2 * h))/2;
1320
1321 return (-2 * o1h + o2h )/(h*h*h);
1322 }
1323 else if constexpr( Order == 4 ) // 4th order derivative
1324 {
1325 constexpr auto h = get_delta<Order>( delta<ArgType>{} );
1326
1327 auto fx = f(x);
1328 auto e1h = (f(x + h) + f(x - h) ) / 2 - fx;
1329 auto e2h = (f(x + 2 * h) + f(x - 2 * h) ) / 2 - fx;
1330
1331 return -2 * (4 * e1h - e2h) / (h*h*h*h);
1332 }
1333 else
1334 {
1335 return seven_point_stencil<Order>(f, x);
1336 }
1337 }
1338 // end of five_point_stencil()
1339
1340 template<auto VarIndex, typename FuncType, cpt::arithmetic_c... ArgTypes>
1341 auto fix_variables_other_than_VarIndex_ed(FuncType&& func, ArgTypes... args) noexcept
1342 {
1343 return [&func, args...](auto v)
1344 {
1345 std::tuple arguments{ args... };
1346
1347 std::get<VarIndex>(arguments) = v;
1348
1349 return std::apply(func, arguments);
1350 };
1351 }
1352
1353 template<auto VarIndex, typename FuncType, cpt::arithmetic_c... ArgTypes>
1354 auto fix_variables_other_than_VarIndex_ed(FuncType&& func, std::tuple<ArgTypes...> args) noexcept
1355 {
1356 auto wrapper = [func](auto... arguments)
1357 {
1358 return fix_variables_other_than_VarIndex_ed<VarIndex>(func, arguments...);
1359 };
1360
1361 return std::apply(wrapper, args);
1362 }
1363
1364 template<auto VarIndex, typename FuncType, cpt::arithmetic_c ArgType, std::size_t N>
1365 auto fix_variables_other_than_VarIndex_ed(FuncType&& func, std::array<ArgType, N> args) noexcept
1366 {
1367 auto wrapper = [func](auto... arguments)
1368 {
1369 return fix_variables_other_than_VarIndex_ed<VarIndex>(func, arguments...);
1370 };
1371
1372 return std::apply(wrapper, args);
1373 }
1374
1375 template<auto VarIndex, auto Order, typename FuncType, cpt::arithmetic_c... ArgTypes>
1376 auto partial_derivative( cpt::sequence<VarIndex, Order>, FuncType&& func, ArgTypes... args) noexcept
1377 {
1378 auto fixed = fix_variables_other_than_VarIndex_ed<VarIndex>(func, args...);
1379 auto arguments = std::forward_as_tuple(args...);
1380 auto v = std::get<VarIndex>(arguments);
1381 return five_point_stencil<Order>(fixed, v);
1382 }
1383
1384 template<auto VarIndex, auto Order, typename FuncType, cpt::arithmetic_c... ArgTypes>
1385 auto partial_derivative( cpt::sequence<VarIndex, Order>, FuncType&& func,
1386 std::tuple<ArgTypes...> args) noexcept
1387 {
1388 auto fixed = fix_variables_other_than_VarIndex_ed<VarIndex>(func, args);
1389 auto v = std::get<VarIndex>(args);
1390 return five_point_stencil<Order>(fixed, v);
1391 }
1392
1393 template<auto VarIndex, auto Order, typename FuncType,
1394 cpt::arithmetic_c ArgType, std::size_t N>
1395 auto partial_derivative( cpt::sequence<VarIndex, Order>,
1396 FuncType&& func, std::array<ArgType, N> args) noexcept
1397 {
1398 auto fixed = fix_variables_other_than_VarIndex_ed<VarIndex>(func, args);
1399 auto v = std::get<VarIndex>(args);
1400 return five_point_stencil<Order>(fixed, v);
1401 }
1402
1403 template<typename VarOdr, typename FuncType, cpt::arithmetic_c... ArgTypes>
1404 auto partial_derivative( cpt::type_container<VarOdr>, FuncType&& func, ArgTypes... args) noexcept
1405 {
1406 constexpr auto VarIndex = cpt::get_nth_value<0>( VarOdr{} ); // x, y, z, etc.
1407 constexpr auto Order = cpt::get_nth_value<1>( VarOdr{} ); // 1, 2, 3, ...
1408
1409 auto fixed = fix_variables_other_than_VarIndex_ed<VarIndex>(func, args...);
1410 auto arguments = std::forward_as_tuple(args...);
1411 auto v = std::get<VarIndex>(arguments);
1412 return five_point_stencil<Order>(fixed, v);
1413 }
1414
1415 template<typename VarOdr, typename FuncType, cpt::arithmetic_c... ArgTypes>
1416 auto partial_derivative( cpt::type_container<VarOdr> cmd,
1417 FuncType&& func, std::tuple<ArgTypes...> const& args) noexcept
1418 {
1419 auto wrapper = [=](auto... arguments)
1420 {
1421 return partial_derivative(cmd, func, arguments...);
1422 };
1423
1424 return std::apply(wrapper, args);
1425 }
1426
1427 template<typename VarOdr, typename FuncType, cpt::arithmetic_c ArgType, std::size_t N>
1428 auto partial_derivative( cpt::type_container<VarOdr> cmd,
1429 FuncType&& func, std::array<ArgType, N> const& args) noexcept
1430 {
1431 auto wrapper = [=](auto... arguments)
1432 {
1433 return partial_derivative(cmd, func, arguments...);
1434 };
1435
1436 return std::apply(wrapper, args);
1437 }
1438
1439 template<typename VarOdr1, typename VarOdr2, typename... VarOdrs,
1440 typename FuncType, cpt::arithmetic_c... ArgTypes>
1441 auto partial_derivative( cpt::type_container<VarOdr1, VarOdr2, VarOdrs...>,
1442 FuncType&& func, ArgTypes... args) noexcept
1443 {
1444 auto right_derivatives = [&func](auto... arguments)
1445 {
1446 return partial_derivative(cpt::type_container<VarOdr2, VarOdrs...>{}, func, arguments...);
1447 };
1448
1449 return partial_derivative( cpt::type_container<VarOdr1>{}, right_derivatives, args...);
1450 }
1451
1452 template<typename VarOdr1, typename VarOdr2, typename... VarOdrs,
1453 typename FuncType, cpt::arithmetic_c... ArgTypes>
1454 auto partial_derivative( cpt::type_container<VarOdr1, VarOdr2, VarOdrs...> cmd,
1455 FuncType&& func, std::tuple<ArgTypes...> const& args) noexcept
1456 {
1457 auto wrapper = [=](auto... arguments)
1458 {
1459 return partial_derivative(cmd, func, arguments...);
1460 };
1461
1462 return std::apply(wrapper, wrapper);
1463 }
1464
1465 template<typename VarOdr1, typename VarOdr2, typename... VarOdrs,
1466 typename FuncType, cpt::arithmetic_c ArgType, std::size_t N>
1467 auto partial_derivative( cpt::type_container<VarOdr1, VarOdr2, VarOdrs...> cmd,
1468 FuncType&& func, std::array<ArgType, N> const& args) noexcept
1469 {
1470 auto wrapper = [=](auto... arguments)
1471 {
1472 return partial_derivative(cmd, func, arguments...);
1473 };
1474
1475 return std::apply(wrapper, wrapper);
1476 }
1477
1478 template<typename FuncType,
1479 typename... VariableTypes, typename ParamType>
1480 auto directional_derivative(FuncType&& f,
1481 std::tuple<VariableTypes...> const& vars, ParamType param)
1482 {
1483 auto f_of_param = [&](auto p)
1484 {
1485 constexpr auto N = sizeof...(VariableTypes);
1486
1487 auto call_func =[&]<auto...i>(cpt::sequence<i...>)
1488 {
1489 return f( std::get<i>(vars)(p)... );
1490 };
1491
1492 return cpt::for_stallion<N>(call_func);
1493 };
1494
1495 return five_point_stencil<1>(f_of_param, param);
1496 }
1497
1498 template<typename... CmdTypes, typename FuncType,
1499 typename... VariableTypes, typename... ParamTypes>
1500 auto parametric_derivative(cpt::type_container<CmdTypes...> cmd,
1501 FuncType && f, std::tuple<VariableTypes...> const& vars, ParamTypes... ps)
1502 {
1503 auto f_of_params = [&](auto... params)
1504 {
1505 constexpr auto N = sizeof...(VariableTypes);
1506
1507 auto call_func =[&]<auto...i>(cpt::sequence<i...>)
1508 {
1509 return f( std::get<i>(vars)(params...)... );
1510 };
1511
1512 return cpt::for_stallion<N>(call_func);
1513 };
1514
1515 return std::array{ partial_derivative(CmdTypes{}, f_of_params, ps...) ... };
1516 }
1517
1518 template<typename... CmdTypes, typename FuncType,
1519 typename... VariableTypes, typename... ParamTypes>
1520 auto parametric_derivative(cpt::type_container<CmdTypes...> cmd,
1521 FuncType && f, std::tuple<VariableTypes...> const& vars, std::tuple<ParamTypes...> ps)
1522 {
1523 auto f_of_params = [&](auto... params)
1524 {
1525 constexpr auto N = sizeof...(VariableTypes);
1526
1527 auto call_func =[&]<auto...i>(cpt::sequence<i...>)
1528 {
1529 return f( std::get<i>(vars)(params...)... );
1530 };
1531
1532 return cpt::for_stallion<N>(call_func);
1533 };
1534
1535 return std::array{ partial_derivative(CmdTypes{}, f_of_params, ps) ... };
1536 }
1537
1538 template<typename... CmdTypes, typename FuncType,
1539 typename... VariableTypes, typename ParamType, std::size_t N>
1540 auto parametric_derivative(cpt::type_container<CmdTypes...> cmd,
1541 FuncType && f, std::tuple<VariableTypes...> const& vars, std::array<ParamType, N> ps)
1542 {
1543 auto f_of_params = [&](auto... params)
1544 {
1545 constexpr auto NN = sizeof...(VariableTypes);
1546
1547 auto call_func =[&]<auto...i>(cpt::sequence<i...>)
1548 {
1549 return f( std::get<i>(vars)(params...)... );
1550 };
1551
1552 return cpt::for_stallion<NN>(call_func);
1553 };
1554
1555 return std::array{ partial_derivative(CmdTypes{}, f_of_params, ps) ... };
1556 }
1557
1558 // Gradient
1559 // https://en.wikipedia.org/wiki/Gradient
1560
1561 template<typename FuncType, cpt::arithmetic_c...ArgTypes>
1562 requires ( std::is_invocable_v<FuncType, ArgTypes...> )
1563 auto gradient(FuncType&& function, ArgTypes... args) noexcept
1564 {
1565 using namespace commands;
1566
1567 using ele_t = std::common_type_t<ArgTypes...>;
1568
1569 auto cmd_dx_1 = create_command(dx_1);
1570 auto cmd_dy_1 = create_command(dy_1);
1571 auto cmd_dz_1 = create_command(dz_1);
1572
1573 auto x = partial_derivative(cmd_dx_1, function, args...);
1574 auto y = partial_derivative(cmd_dy_1, function, args...);
1575 auto z = partial_derivative(cmd_dz_1, function, args...);
1576
1577
1578 return std::array{ adjust_zero<ele_t>(x),
1579 adjust_zero<ele_t>(y),
1580 adjust_zero<ele_t>(z)};
1581 }
1582
1583 template<typename FuncType, cpt::arithmetic_c...ArgTypes>
1584 requires ( std::is_invocable_v<FuncType, ArgTypes...> )
1585 auto gradient(FuncType&& function, const std::tuple<ArgTypes...>& args) noexcept
1586 {
1587 auto wrapper = [&](auto... arguments)
1588 {
1589 return gradient(function, arguments...);
1590 };
1591
1592 return std::apply(wrapper, args);
1593 }
1594
1595 template<typename FuncType, cpt::arithmetic_c ArgType, std::size_t N>
1596 auto gradient(FuncType&& function, const std::array<ArgType, N>& args) noexcept
1597 {
1598 auto wrapper = [&](auto... arguments)
1599 {
1600 return gradient(function, arguments...);
1601 };
1602
1603 return std::apply(wrapper, args);
1604 }
1605
1606 // Curl (mathematics)
1607 // https://en.wikipedia.org/wiki/Curl_(mathematics)
1608 template<typename FuncTypeX, typename FuncTypeY,
1609 typename FuncTypeZ, cpt::arithmetic_c...ArgTypes>
1610 requires ( std::is_invocable_v<FuncTypeX, ArgTypes...> &&
1611 std::is_invocable_v<FuncTypeY, ArgTypes...> &&
1612 std::is_invocable_v<FuncTypeZ, ArgTypes...> )
1613 auto curl(FuncTypeX&& func_x, FuncTypeY&& func_y,
1614 FuncTypeZ&& func_z, ArgTypes... args) noexcept
1615 {
1616 using namespace commands;
1617
1618 using ele_t = std::common_type_t<ArgTypes...>;
1619
1620 auto cmd_dx = create_command(dx_1);
1621 auto cmd_dy = create_command(dy_1);
1622 auto cmd_dz = create_command(dz_1);
1623
1624 auto x = partial_derivative(cmd_dy, func_z, args...)
1625 - partial_derivative(cmd_dz, func_y, args...);
1626
1627 auto y = partial_derivative(cmd_dz, func_x, args...)
1628 - partial_derivative(cmd_dx, func_z, args...);
1629
1630 auto z = partial_derivative(cmd_dx, func_y, args...)
1631 - partial_derivative(cmd_dy, func_x, args...);
1632
1633 return std::array{ adjust_zero<ele_t>(x), adjust_zero<ele_t>(y), adjust_zero<ele_t>(z) };
1634 }
1635
1636 template<typename FuncTypeX, typename FuncTypeY,
1637 typename FuncTypeZ, cpt::arithmetic_c...ArgTypes>
1638 requires ( std::is_invocable_v<FuncTypeX, ArgTypes...> &&
1639 std::is_invocable_v<FuncTypeY, ArgTypes...> &&
1640 std::is_invocable_v<FuncTypeZ, ArgTypes...> )
1641 auto curl(FuncTypeX&& func_x, FuncTypeY&& func_y,
1642 FuncTypeZ&& func_z, const std::tuple<ArgTypes...>& args) noexcept
1643 {
1644 auto wrapper = [&](auto... arguments)
1645 {
1646 return curl(func_x, func_y, func_z, arguments...);
1647 };
1648
1649 return std::apply(wrapper, args);
1650 }
1651
1652 template<typename FuncTypeX, typename FuncTypeY,
1653 typename FuncTypeZ, cpt::arithmetic_c ArgType, std::size_t N>
1654 auto curl(FuncTypeX&& func_x, FuncTypeY&& func_y,
1655 FuncTypeZ&& func_z, const std::array<ArgType, N>& args) noexcept
1656 {
1657 auto wrapper = [&](auto... arguments)
1658 {
1659 return curl(func_x, func_y, func_z, arguments...);
1660 };
1661
1662 return std::apply(wrapper, args);
1663 }
1664
1665 // Divergence
1666 // https://en.wikipedia.org/wiki/Divergence
1667 template<typename FuncTypeX, typename FuncTypeY,
1668 typename FuncTypeZ, cpt::arithmetic_c...ArgTypes>
1669 requires ( std::is_invocable_v<FuncTypeX, ArgTypes...> &&
1670 std::is_invocable_v<FuncTypeY, ArgTypes...> &&
1671 std::is_invocable_v<FuncTypeZ, ArgTypes...> )
1672 auto divergence(FuncTypeX&& func_x, FuncTypeY&& func_y,
1673 FuncTypeZ&& func_z, ArgTypes... args) noexcept
1674 {
1675 using namespace commands;
1676
1677 using ele_t = std::common_type_t<ArgTypes...>;
1678
1679 auto cmd_dx = create_command(dx_1);
1680 auto cmd_dy = create_command(dy_1);
1681 auto cmd_dz = create_command(dz_1);
1682
1683 auto d = partial_derivative(cmd_dx, func_x, args...) +
1684 partial_derivative(cmd_dy, func_y, args...) +
1685 partial_derivative(cmd_dz, func_z, args...);
1686
1687 return adjust_zero<ele_t>(d);
1688 }
1689
1690 template<typename FuncTypeX, typename FuncTypeY,
1691 typename FuncTypeZ, cpt::arithmetic_c...ArgTypes>
1692 requires ( std::is_invocable_v<FuncTypeX, ArgTypes...> &&
1693 std::is_invocable_v<FuncTypeY, ArgTypes...> &&
1694 std::is_invocable_v<FuncTypeZ, ArgTypes...> )
1695 auto divergence(FuncTypeX&& func_x, FuncTypeY&& func_y,
1696 FuncTypeZ&& func_z, const std::tuple<ArgTypes...>& args) noexcept
1697 {
1698 auto wrapper = [&](auto... arguments)
1699 {
1700 return divergence(func_x, func_y, func_z, arguments...);
1701 };
1702
1703 return std::apply(wrapper, args);
1704 }
1705
1706 template<typename FuncTypeX, typename FuncTypeY,
1707 typename FuncTypeZ, cpt::arithmetic_c ArgType, std::size_t N>
1708 auto divergence(FuncTypeX&& func_x, FuncTypeY&& func_y,
1709 FuncTypeZ&& func_z, const std::array<ArgType, N>& args) noexcept
1710 {
1711 auto wrapper = [&](auto... arguments)
1712 {
1713 return divergence(func_x, func_y, func_z, arguments...);
1714 };
1715
1716 return std::apply(wrapper, args);
1717 }
1718
1719 #ifdef CPG_INCLUDE_SYCL
1720
1722 // single variable - evaluateion
1723 template<typename FuncType, typename BoundType>
1724 auto evaluate(sycl::queue& queue, FuncType func,
1725 std::size_t N, std::array<BoundType, 2> bound)
1726 {
1727 using arg_t = BoundType;
1728
1729 using array_t = std::array<BoundType, 2>;
1730
1731 std::vector<array_t> result(N);
1732
1733 {
1734 sycl::buffer<array_t, 1> buffer{ result.data(), sycl::range<1>{ N } };
1735
1736 auto lower = std::get<0>(bound);
1737 auto upper = std::get<1>(bound);
1738
1739 auto dx = (upper - lower)/(N - 1);
1740
1741 queue.submit
1742 (
1743 [&](sycl::handler& cgh)
1744 {
1745 sycl::accessor v{ buffer, cgh, sycl::write_only, sycl::no_init };
1746
1747 cgh.parallel_for(N,
1748 [=](auto idx)
1749 {
1750 auto coord = adjust_zero<BoundType>(lower + dx * (int)idx);
1751 v[idx] = std::array<BoundType, 2>{ coord, adjust_zero<BoundType>( func(coord) ) };
1752 }
1753 );
1754 }
1755 );
1756 }
1757
1758 return result;
1759 }
1760
1762 // single variable - differentiation
1763 template<auto VarIndex, auto DerivativeOrder, typename FuncType, typename BoundType>
1764 auto differentiate(sycl::queue& queue, cpt::sequence<VarIndex, DerivativeOrder> command,
1765 FuncType func, std::size_t N, std::array<BoundType, 2> bound)
1766 {
1767 using arg_t = BoundType;
1768
1769 auto deriv = [func](auto value)
1770 {
1771 return five_point_stencil<DerivativeOrder>(func, value);
1772 };
1773
1774 return evaluate(queue, deriv, N, bound);
1775 }
1776
1777
1779 // multivariable - evaluateion
1780 template<std::size_t Index, typename FuncType, typename BoundType>
1781 auto evaluate(sycl::queue& queue, FuncType func,
1782 std::size_t N, std::array<BoundType, 2> bound, auto... args)
1783 {
1784 using arg_t = BoundType;
1785
1786 using array_t = std::array<BoundType, sizeof...(args) + 1>;
1787
1788 std::vector<array_t> result(N);
1789
1790 {
1791 sycl::buffer<array_t, 1> buffer{ result.data(), sycl::range<1>{ N } };
1792
1793 auto lower = std::get<0>(bound);
1794 auto upper = std::get<1>(bound);
1795
1796 auto dx = (upper - lower)/(N - 1);
1797
1798 queue.submit
1799 (
1800 [&](sycl::handler& cgh)
1801 {
1802 sycl::accessor v{ buffer, cgh, sycl::write_only, sycl::no_init };
1803
1804 cgh.parallel_for(N,
1805 [=](auto idx)
1806 {
1807 auto coord = adjust_zero<BoundType>(lower + dx * (int)idx);
1808
1809 if constexpr(requires{ func(coord); })
1810 {
1811 v[idx] = std::array{ (arg_t)args..., adjust_zero<BoundType>(func(coord)) };
1812 std::get<Index>(v[idx]) = coord;
1813 }
1814 else if constexpr(requires{ func(args...); })
1815 {
1816 auto arguments = std::tuple{ args ... };
1817 std::get<Index>(arguments) = coord;
1818
1819 v[idx] = std::array{ (arg_t)args..., adjust_zero<BoundType>(std::apply(func, arguments)) };
1820 std::get<Index>(v[idx]) = coord;
1821 }
1822 }
1823 );
1824 }
1825 );
1826 }
1827
1828 return result;
1829 }
1830
1832 // multivariable - differentiation
1833 template<auto VarIndex, typename IndexType, typename... IndexTypes, typename FuncType, typename BoundType>
1834 auto differentiate(sycl::queue& queue,
1835 cpt::type_container<IndexType, IndexTypes...> command,
1836 FuncType func, std::size_t N, BoundType bound, auto... args)
1837 {
1838 using real_t = std::common_type_t<decltype(args)...>;
1839
1840 auto deriv = [command, func, args...](auto value)
1841 {
1842 auto arguments = std::tuple{args...};
1843 std::get<VarIndex>(arguments) = value;
1844
1845 return partial_derivative(command, func, arguments);
1846 };
1847
1848 return evaluate<VarIndex>(queue, deriv, N, bound, args...);
1849 }
1850
1851 // if CountX, CountY, or CountZ is 1,
1852 // the center is chosen
1853 template<std::size_t CountX, std::size_t CountY, std::size_t CountZ,
1854 typename FuncType, typename BoundType>
1855 auto gradients(sycl::queue& queue, FuncType&& func,
1856 std::array<BoundType, 2> bound_x, // should be copy-semantic
1857 std::array<BoundType, 2> bound_y, // should be copy-semantic
1858 std::array<BoundType, 2> bound_z) // should be copy-semantic
1859 requires requires { func( BoundType{}, BoundType{}, BoundType{}); }
1860 {
1861 using vertex_t = decltype(std::tuple{ std::array<BoundType, 3>{},
1862 gradient(func, std::tuple<BoundType, BoundType, BoundType>{}) });
1863
1864 using vertices_t = std::vector<vertex_t>;
1865
1866 vertices_t vertices(CountX * CountY * CountZ);
1867
1868 auto& v =
1869 cpt::cast_ref(std::index_sequence<CountX, CountY, CountZ>{}, vertices);
1870
1871 BoundType lower[3], upper[3], delta[3];
1872
1873 // modifies bounds in compute_delta()
1874 delta[0] = compute_delta(CountX, bound_x);
1875 delta[1] = compute_delta(CountY, bound_y);
1876 delta[2] = compute_delta(CountZ, bound_z);
1877
1878 lower[0] = std::get<0>(bound_x);
1879 upper[0] = std::get<1>(bound_x);
1880
1881 lower[1] = std::get<0>(bound_y);
1882 upper[1] = std::get<1>(bound_y);
1883
1884 lower[2] = std::get<0>(bound_z);
1885 upper[2] = std::get<1>(bound_z);
1886
1887 auto range = sycl::range<3>{CountX, CountY, CountZ};
1888
1889 sycl::buffer<BoundType, 1> buf_lower{ lower, 3 };
1890 sycl::buffer<BoundType, 1> buf_delta{ delta, 3 };
1891
1892 {
1893 sycl::buffer<vertex_t, 3> buf_vertices{ vertices.data(), range};
1894
1895 // cgh stands for command group handler
1896 auto command_group = [&](sycl::handler& cgh)
1897 {
1898 sycl::accessor acc_lower{ buf_lower, cgh, sycl::read_only };
1899 sycl::accessor acc_delta{ buf_delta, cgh, sycl::read_only };
1900 sycl::accessor acc_v{ buf_vertices, cgh, sycl::write_only, sycl::no_init};
1901
1902 auto kernel = [=](sycl::id<3> idx)
1903 {
1904 auto i = idx[0]; auto j = idx[1]; auto k = idx[2];
1905
1906 auto x = adjust_zero<BoundType>(acc_lower[0] + acc_delta[0] * i);
1907 auto y = adjust_zero<BoundType>(acc_lower[1] + acc_delta[1] * j);
1908 auto z = adjust_zero<BoundType>(acc_lower[2] + acc_delta[2] * k);
1909
1910 acc_v[idx] = vertex_t{std::array{x, y, z},
1911 gradient(func, std::tuple{x, y, z}) };
1912 };
1913
1914 cgh.parallel_for(range, kernel);
1915 // tbb::parallel_for(range3d, compute_gradient);
1916 };
1917
1918 queue.submit(command_group);
1919 }
1920
1921 return vertices;
1922 }
1923
1924 // if CountX, CountY, or CountZ is 1,
1925 // the center is chosen
1926 template<std::size_t CountX, std::size_t CountY, std::size_t CountZ,
1927 typename FuncTypeX, typename FuncTypeY, typename FuncTypeZ,
1928 typename BoundType>
1929 auto curls(sycl::queue& queue, FuncTypeX&& func_x, FuncTypeY&& func_y, FuncTypeZ&& func_z,
1930 std::array<BoundType, 2> bound_x, // should be copy-semantic
1931 std::array<BoundType, 2> bound_y, // should be copy-semantic
1932 std::array<BoundType, 2> bound_z) // should be copy-semantic
1933 requires requires
1934 { func_x( BoundType{}, BoundType{}, BoundType{});
1935 func_y( BoundType{}, BoundType{}, BoundType{});
1936 func_z( BoundType{}, BoundType{}, BoundType{}); }
1937 {
1938 using vertex_t = decltype(std::tuple{ std::array<BoundType, 3>{},
1939 curl(func_x, func_y, func_z, std::tuple<BoundType, BoundType, BoundType>{}) });
1940
1941 using vertices_t = std::vector<vertex_t>;
1942
1943 vertices_t vertices(CountX * CountY * CountZ);
1944
1945 auto& v =
1946 cpt::cast_ref(std::index_sequence<CountX, CountY, CountZ>{}, vertices);
1947
1948 BoundType lower[3], upper[3], delta[3];
1949
1950 // modifies bounds in compute_delta()
1951 delta[0] = compute_delta(CountX, bound_x);
1952 delta[1] = compute_delta(CountY, bound_y);
1953 delta[2] = compute_delta(CountZ, bound_z);
1954
1955 lower[0] = std::get<0>(bound_x);
1956 upper[0] = std::get<1>(bound_x);
1957
1958 lower[1] = std::get<0>(bound_y);
1959 upper[1] = std::get<1>(bound_y);
1960
1961 lower[2] = std::get<0>(bound_z);
1962 upper[2] = std::get<1>(bound_z);
1963
1964 auto range = sycl::range<3>{CountX, CountY, CountZ};
1965
1966 sycl::buffer<BoundType, 1> buf_lower{ lower, 3 };
1967 sycl::buffer<BoundType, 1> buf_delta{ delta, 3 };
1968
1969 {
1970 sycl::buffer<vertex_t, 3> buf_vertices{ vertices.data(), range};
1971
1972 // cgh stands for command group handler
1973 auto command_group = [&](sycl::handler& cgh)
1974 {
1975 sycl::accessor acc_lower{ buf_lower, cgh, sycl::read_only };
1976 sycl::accessor acc_delta{ buf_delta, cgh, sycl::read_only };
1977 sycl::accessor acc_v{ buf_vertices, cgh, sycl::write_only, sycl::no_init};
1978
1979 auto kernel = [=](sycl::id<3> idx)
1980 {
1981 auto i = idx[0]; auto j = idx[1]; auto k = idx[2];
1982
1983 auto x = adjust_zero<BoundType>(acc_lower[0] + acc_delta[0] * i);
1984 auto y = adjust_zero<BoundType>(acc_lower[1] + acc_delta[1] * j);
1985 auto z = adjust_zero<BoundType>(acc_lower[2] + acc_delta[2] * k);
1986
1987 acc_v[idx] = vertex_t{std::array{x, y, z},
1988 curl(func_x, func_y, func_z, std::tuple{x, y, z}) };
1989 };
1990
1991 cgh.parallel_for(range, kernel);
1992 // tbb::parallel_for(range3d, compute_gradient);
1993 };
1994
1995 queue.submit(command_group);
1996 }
1997
1998 return vertices;
1999 }
2000
2001 // if CountX, CountY, or CountZ is 1,
2002 // the center is chosen
2003 template<std::size_t CountX, std::size_t CountY, std::size_t CountZ,
2004 typename FuncTypeX, typename FuncTypeY, typename FuncTypeZ,
2005 typename BoundType>
2006 auto divs(sycl::queue& queue, FuncTypeX&& func_x, FuncTypeY&& func_y, FuncTypeZ&& func_z,
2007 std::array<BoundType, 2> bound_x, // should be copy-semantic
2008 std::array<BoundType, 2> bound_y, // should be copy-semantic
2009 std::array<BoundType, 2> bound_z) // should be copy-semantic
2010 requires requires
2011 { func_x( BoundType{}, BoundType{}, BoundType{});
2012 func_y( BoundType{}, BoundType{}, BoundType{});
2013 func_z( BoundType{}, BoundType{}, BoundType{}); }
2014 {
2015 using vertex_t = std::array<BoundType, 4>;
2016
2017 using vertices_t = std::vector<vertex_t>;
2018
2019 vertices_t vertices(CountX * CountY * CountZ);
2020
2021 auto& v =
2022 cpt::cast_ref(std::index_sequence<CountX, CountY, CountZ>{}, vertices);
2023
2024 BoundType lower[3], upper[3], delta[3];
2025
2026 // modifies bounds in compute_delta()
2027 delta[0] = compute_delta(CountX, bound_x);
2028 delta[1] = compute_delta(CountY, bound_y);
2029 delta[2] = compute_delta(CountZ, bound_z);
2030
2031 lower[0] = std::get<0>(bound_x);
2032 upper[0] = std::get<1>(bound_x);
2033
2034 lower[1] = std::get<0>(bound_y);
2035 upper[1] = std::get<1>(bound_y);
2036
2037 lower[2] = std::get<0>(bound_z);
2038 upper[2] = std::get<1>(bound_z);
2039
2040 auto range = sycl::range<3>{CountX, CountY, CountZ};
2041
2042 sycl::buffer<BoundType, 1> buf_lower{ lower, 3 };
2043 sycl::buffer<BoundType, 1> buf_delta{ delta, 3 };
2044
2045 {
2046 sycl::buffer<vertex_t, 3> buf_vertices{ vertices.data(), range};
2047
2048 // cgh stands for command group handler
2049 auto command_group = [&](sycl::handler& cgh)
2050 {
2051 sycl::accessor acc_lower{ buf_lower, cgh, sycl::read_only };
2052 sycl::accessor acc_delta{ buf_delta, cgh, sycl::read_only };
2053 sycl::accessor acc_v{ buf_vertices, cgh, sycl::write_only, sycl::no_init};
2054
2055 auto kernel = [=](sycl::id<3> idx)
2056 {
2057 auto i = idx[0]; auto j = idx[1]; auto k = idx[2];
2058
2059 auto x = adjust_zero<BoundType>(acc_lower[0] + acc_delta[0] * i);
2060 auto y = adjust_zero<BoundType>(acc_lower[1] + acc_delta[1] * j);
2061 auto z = adjust_zero<BoundType>(acc_lower[2] + acc_delta[2] * k);
2062
2063 acc_v[idx] = std::array{x, y, z,
2064 divergence(func_x, func_y, func_z, std::tuple{x, y, z}) };
2065 };
2066
2067 cgh.parallel_for(range, kernel);
2068 // tbb::parallel_for(range3d, compute_gradient);
2069 };
2070
2071 queue.submit(command_group);
2072 }
2073
2074 return vertices;
2075 }
2076
2077 #endif
2078 // End of SYCL coding
2079
2081 // single variable - evaluation
2082 template<typename FuncType, typename BoundType>
2083 auto evaluate(FuncType&& func, std::size_t N,
2084 std::array<BoundType, 2> bound) requires requires { func(BoundType{}); }
2085 {
2086 using arg_t = BoundType;
2087 using array_t = std::array<arg_t, 2>;
2088
2089 std::vector<array_t> result(N);
2090
2091 auto lower = std::get<0>(bound);
2092 auto upper = std::get<1>(bound);
2093
2094 auto dx = (upper - lower)/ (N - 1);
2095
2096 auto process = [&func, &result, lower, dx](auto range)
2097 {
2098 for(auto i = range.begin(); i < range.end(); ++i)
2099 {
2100 auto coord = adjust_zero<BoundType>(lower + i * dx);
2101 result[i] = array_t{ coord, adjust_zero<BoundType>(func( coord )) };
2102 }
2103 };
2104
2105 tbb::parallel_for( tbb::blocked_range{ 0, (int)N }, process);
2106
2107 return result;
2108 }
2109
2111 // multivariable - evaluation
2112 template<std::size_t VarIndex, typename FuncType, typename BoundType, typename...ArgTypes>
2113 auto evaluate(FuncType&& func, std::size_t N, std::array<BoundType, 2> bound, ArgTypes... args)
2114 {
2115 using arg_t = BoundType;
2116 using array_t = std::array<BoundType, sizeof...(args) + 1>;
2117
2118 std::vector<array_t> result(N);
2119
2120 auto lower = std::get<0>(bound);
2121 auto upper = std::get<1>(bound);
2122
2123 auto dx = (upper - lower)/ (N - 1);
2124
2125 auto process = [&func, &result, lower, dx, args...](auto range)
2126 {
2127 auto arguments = std::tuple{args...};
2128
2129 for(auto i = range.begin(); i < range.end(); ++i)
2130 {
2131 auto coord = adjust_zero<BoundType>(lower + i * dx);
2132
2133 std::get<VarIndex>(arguments) = coord;
2134
2135 result[i] = array_t{ (arg_t)args ...,
2136 adjust_zero<BoundType>(std::apply(func, arguments)) };
2137
2138 std::get<VarIndex>(result[i]) = coord;
2139 }
2140 };
2141
2142 tbb::parallel_for( tbb::blocked_range{ 0, (int)N }, process);
2143
2144 return result;
2145 }
2146
2148 // single variable - differentiation
2149 template<auto VarIndex, auto DerivativeOrder, typename FuncType, typename BoundType>
2150 auto differentiate(cpt::sequence<VarIndex, DerivativeOrder>,
2151 FuncType&& func, std::size_t N, std::array<BoundType, 2> bound)
2152 requires requires { func(BoundType{}); }
2153 {
2154 using arg_t = BoundType;
2155
2156 auto ff = [&func](auto value)
2157 {
2158 return five_point_stencil<DerivativeOrder>(func, value);
2159 };
2160
2161 using array_t = std::array<arg_t, 2>;
2162
2163 std::vector<array_t> result(N);
2164
2165 auto lower = std::get<0>(bound);
2166 auto upper = std::get<1>(bound);
2167
2168 auto dx = (upper - lower)/ (N - 1);
2169
2170 auto process = [&ff, &result, lower, dx](auto range)
2171 {
2172 for(auto i = range.begin(); i < range.end(); ++i)
2173 {
2174 auto coord = adjust_zero<BoundType>(lower + i * dx);
2175
2176 result[i] = array_t{ coord, adjust_zero<BoundType>(ff(coord))};
2177 }
2178 };
2179
2180 tbb::parallel_for( tbb::blocked_range{ 0, (int)N }, process);
2181
2182 return result;
2183 }
2184
2186 // multivariable - differentiation
2187 template<auto VarIndex, typename IndexType, typename... IndexTypes,
2188 typename FuncType, typename BoundType>
2189 auto differentiate(cpt::type_container<IndexType, IndexTypes...> command,
2190 FuncType&& func, std::size_t N, std::array<BoundType, 2> bound,
2191 auto... args) requires requires { func(args...); }
2192 {
2193 using arg_t = std::common_type_t<decltype(args)...>;
2194
2195 auto ff = [&func, args...](auto value)
2196 {
2197 auto arguments = std::tuple{args...};
2198 std::get<VarIndex>(arguments) = value;
2199
2200 return partial_derivative(cpt::type_container<IndexType, IndexTypes...>{}, func, arguments);
2201 };
2202
2203 // using tuple_t = std::tuple<arg_t, decltype(args)...>;
2204 using array_t = decltype(std::array{ (arg_t)args..., arg_t{}});
2205
2206 std::vector<array_t> result(N);
2207
2208 auto lower = std::get<0>(bound);
2209 auto upper = std::get<1>(bound);
2210
2211 auto dx = (upper - lower)/ (N - 1);
2212
2213 auto process = [&ff, &result, lower, dx, args...](auto range)
2214 {
2215 for(auto i = range.begin(); i < range.end(); ++i)
2216 {
2217 auto coord = adjust_zero<BoundType>(lower + i * dx);
2218
2219 result[i] = array_t{ args ..., adjust_zero<BoundType>(ff(coord))};
2220
2221 std::get<VarIndex>(result[i]) = coord;
2222 }
2223 };
2224
2225 tbb::parallel_for( tbb::blocked_range{ 0, (int)N }, process);
2226
2227 return result;
2228 }
2229
2230 // if CountX, CountY, or CountZ is 1,
2231 // the center is chosen
2232 template<std::size_t CountX, std::size_t CountY, std::size_t CountZ,
2233 typename FuncType, typename BoundType>
2234 auto gradients(FuncType&& func,
2235 std::array<BoundType, 2> bound_x, // should be copy-semantic
2236 std::array<BoundType, 2> bound_y, // should be copy-semantic
2237 std::array<BoundType, 2> bound_z) // should be copy-semantic
2238 requires requires { func( BoundType{}, BoundType{}, BoundType{}); }
2239 {
2240 using vertex_t = decltype(std::tuple{ std::array<BoundType, 3>{},
2241 gradient(func, std::tuple<BoundType, BoundType, BoundType>{}) });
2242
2243 using vertices_t = std::vector<vertex_t>;
2244
2245 vertices_t vertices(CountX * CountY * CountZ);
2246
2247 auto& v =
2248 cpt::cast_ref(std::index_sequence<CountX, CountY, CountZ>{}, vertices);
2249
2250 BoundType lower[3], upper[3], delta[3];
2251
2252 // modifies bounds in compute_delta()
2253 delta[0] = compute_delta(CountX, bound_x);
2254 delta[1] = compute_delta(CountY, bound_y);
2255 delta[2] = compute_delta(CountZ, bound_z);
2256
2257 lower[0] = std::get<0>(bound_x);
2258 upper[0] = std::get<1>(bound_x);
2259
2260 lower[1] = std::get<0>(bound_y);
2261 upper[1] = std::get<1>(bound_y);
2262
2263 lower[2] = std::get<0>(bound_z);
2264 upper[2] = std::get<1>(bound_z);
2265
2266 auto compute_gradient = [&lower, &delta, &v, &func]
2267 (tbb::blocked_range3d<std::size_t>& range)
2268 {
2269 auto& pages = range.pages();
2270 auto& rows = range.rows();
2271 auto& cols = range.cols();
2272
2273 for(auto i = pages.begin(); i != pages.end(); ++i)
2274 for(auto j = rows.begin(); j != rows.end(); ++j)
2275 for(auto k = cols.begin(); k != cols.end(); ++k)
2276 {
2277 auto x = adjust_zero<BoundType>(lower[0] + delta[0] * i);
2278 auto y = adjust_zero<BoundType>(lower[1] + delta[1] * j);
2279 auto z = adjust_zero<BoundType>(lower[2] + delta[2] * k);
2280
2281 v[i][j][k] = vertex_t{std::array{x, y, z},
2282 gradient(func, std::tuple{x, y, z}) };
2283 }
2284 };
2285
2286 // blocked_range3d
2287 // https://spec.oneapi.io/versions/latest/elements/oneTBB/source/algorithms/blocked_ranges/blocked_range3d_cls.html
2288 auto range3d = tbb::blocked_range3d<std::size_t>(0, CountX, 0, CountY, 0, CountZ);
2289
2290 tbb::parallel_for(range3d, compute_gradient);
2291
2292 return vertices;
2293 }
2294
2295 // if CountX, CountY, or CountZ is 1,
2296 // the center is chosen
2297 template<std::size_t CountX, std::size_t CountY, std::size_t CountZ,
2298 typename FuncTypeX, typename FuncTypeY, typename FuncTypeZ,
2299 typename BoundType>
2300 auto curls(FuncTypeX&& func_x, FuncTypeY&& func_y, FuncTypeZ&& func_z,
2301 std::array<BoundType, 2> bound_x, // should be copy-semantic
2302 std::array<BoundType, 2> bound_y, // should be copy-semantic
2303 std::array<BoundType, 2> bound_z) // should be copy-semantic
2304 requires requires
2305 { func_x( BoundType{}, BoundType{}, BoundType{});
2306 func_y( BoundType{}, BoundType{}, BoundType{});
2307 func_z( BoundType{}, BoundType{}, BoundType{}); }
2308 {
2309 using vertex_t = decltype(std::tuple{ std::array<BoundType, 3>{},
2310 curl(func_x, func_y, func_z, std::tuple<BoundType, BoundType, BoundType>{}) });
2311
2312 using vertices_t = std::vector<vertex_t>;
2313
2314 vertices_t vertices(CountX * CountY * CountZ);
2315
2316 auto& v =
2317 cpt::cast_ref(std::index_sequence<CountX, CountY, CountZ>{}, vertices);
2318
2319 BoundType lower[3], upper[3], delta[3];
2320
2321 // modifies bounds in compute_delta()
2322 delta[0] = compute_delta(CountX, bound_x);
2323 delta[1] = compute_delta(CountY, bound_y);
2324 delta[2] = compute_delta(CountZ, bound_z);
2325
2326 lower[0] = std::get<0>(bound_x);
2327 upper[0] = std::get<1>(bound_x);
2328
2329 lower[1] = std::get<0>(bound_y);
2330 upper[1] = std::get<1>(bound_y);
2331
2332 lower[2] = std::get<0>(bound_z);
2333 upper[2] = std::get<1>(bound_z);
2334
2335 auto compute_gradient = [&lower, &delta, &v,
2336 &func_x, &func_y, &func_z]
2337 (tbb::blocked_range3d<std::size_t>& range)
2338 {
2339 auto& pages = range.pages();
2340 auto& rows = range.rows();
2341 auto& cols = range.cols();
2342
2343 for(auto i = pages.begin(); i != pages.end(); ++i)
2344 for(auto j = rows.begin(); j != rows.end(); ++j)
2345 for(auto k = cols.begin(); k != cols.end(); ++k)
2346 {
2347 auto x = adjust_zero<BoundType>(lower[0] + delta[0] * i);
2348 auto y = adjust_zero<BoundType>(lower[1] + delta[1] * j);
2349 auto z = adjust_zero<BoundType>(lower[2] + delta[2] * k);
2350
2351 v[i][j][k] = vertex_t{std::array{x, y, z},
2352 curl(func_x, func_y, func_z, std::tuple{x, y, z}) };
2353 }
2354 };
2355
2356 // blocked_range3d
2357 // https://spec.oneapi.io/versions/latest/elements/oneTBB/source/algorithms/blocked_ranges/blocked_range3d_cls.html
2358 auto range3d = tbb::blocked_range3d<std::size_t>(0, CountX, 0, CountY, 0, CountZ);
2359
2360 tbb::parallel_for(range3d, compute_gradient);
2361
2362 return vertices;
2363 }
2364
2365 // if CountX, CountY, or CountZ is 1,
2366 // the center is chosen
2367 template<std::size_t CountX, std::size_t CountY, std::size_t CountZ,
2368 typename FuncTypeX, typename FuncTypeY, typename FuncTypeZ,
2369 typename BoundType>
2370 auto divs(FuncTypeX&& func_x, FuncTypeY&& func_y, FuncTypeZ&& func_z,
2371 std::array<BoundType, 2> bound_x, // should be copy-semantic
2372 std::array<BoundType, 2> bound_y, // should be copy-semantic
2373 std::array<BoundType, 2> bound_z) // should be copy-semantic
2374 requires requires
2375 { func_x( BoundType{}, BoundType{}, BoundType{});
2376 func_y( BoundType{}, BoundType{}, BoundType{});
2377 func_z( BoundType{}, BoundType{}, BoundType{}); }
2378 {
2379 using vertex_t = std::array<BoundType, 4>;
2380
2381 using vertices_t = std::vector<vertex_t>;
2382
2383 vertices_t vertices(CountX * CountY * CountZ);
2384
2385 auto& v =
2386 cpt::cast_ref(std::index_sequence<CountX, CountY, CountZ>{}, vertices);
2387
2388 BoundType lower[3], upper[3], delta[3];
2389
2390 // modifies bounds in compute_delta()
2391 delta[0] = compute_delta(CountX, bound_x);
2392 delta[1] = compute_delta(CountY, bound_y);
2393 delta[2] = compute_delta(CountZ, bound_z);
2394
2395 lower[0] = std::get<0>(bound_x);
2396 upper[0] = std::get<1>(bound_x);
2397
2398 lower[1] = std::get<0>(bound_y);
2399 upper[1] = std::get<1>(bound_y);
2400
2401 lower[2] = std::get<0>(bound_z);
2402 upper[2] = std::get<1>(bound_z);
2403
2404 auto compute_gradient = [&lower, &delta, &v,
2405 &func_x, &func_y, &func_z]
2406 (tbb::blocked_range3d<std::size_t>& range)
2407 {
2408 auto& pages = range.pages();
2409 auto& rows = range.rows();
2410 auto& cols = range.cols();
2411
2412 for(auto i = pages.begin(); i != pages.end(); ++i)
2413 for(auto j = rows.begin(); j != rows.end(); ++j)
2414 for(auto k = cols.begin(); k != cols.end(); ++k)
2415 {
2416 auto x = adjust_zero<BoundType>(lower[0] + delta[0] * i);
2417 auto y = adjust_zero<BoundType>(lower[1] + delta[1] * j);
2418 auto z = adjust_zero<BoundType>(lower[2] + delta[2] * k);
2419
2420 v[i][j][k] = std::array{x, y, z,
2421 divergence(func_x, func_y, func_z, std::tuple{x, y, z}) };
2422 }
2423 };
2424
2425 // blocked_range3d
2426 // https://spec.oneapi.io/versions/latest/elements/oneTBB/source/algorithms/blocked_ranges/blocked_range3d_cls.html
2427 auto range3d = tbb::blocked_range3d<std::size_t>(0, CountX, 0, CountY, 0, CountZ);
2428
2429 tbb::parallel_for(range3d, compute_gradient);
2430
2431 return vertices;
2432 }
2433
2435 template<typename FuncType, cpt::arithmetic_c BoundType>
2436 BoundType simpson_rule(FuncType&& f, BoundType a, BoundType b) noexcept
2437 {
2438 if constexpr( std::same_as<BoundType, float> ||
2439 std::same_as<BoundType, double> ||
2440 std::same_as<BoundType, long double> )
2441 {
2442 /*
2443 Simpson's 3/8 rule
2444 https://en.wikipedia.org/wiki/Simpson%27s_rule
2445 */
2446 return (b - a) * ( f(a) + 3 * f( (2*a + b)/3 )
2447 + 3 * f( (a + 2*b)/3 ) + f(b) ) / 8;
2448 }
2449 else
2450 {
2451 static_assert( std::same_as<BoundType, float> ||
2452 std::same_as<BoundType, double> ||
2453 std::same_as<BoundType, long double>);
2454 }
2455 }
2456
2457 /*
2458 Adaptive Quadrature | Lecture 41 | Vector Calculus for Engineers
2459 https://youtu.be/U4NUXAwwR8E?t=41
2460
2461 ch4 A: Adaptive Simpson's Quadrature. Wen Shen
2462 https://www.youtube.com/watch?v=M-AiOqmleaI&t=601s
2463 */
2464 template<typename FuncType, cpt::arithmetic_c BoundType>
2465 BoundType adaptive_simpson_quadrature(FuncType&& f, BoundType a, BoundType b)
2466 {
2467 constexpr auto tolerance =
2468 std::numeric_limits<BoundType>::epsilon();
2469
2470 auto c = ( a + b ) / 2; // the center point between [a, b]
2471
2472 auto s1 = simpson_rule(f, a, b); // [a, b]
2473 auto s2 = simpson_rule(f, a, c) + simpson_rule(f, c, b); // [a, c] + [c, b]
2474
2475 auto diff = std::abs(s1 - s2);
2476
2477 /*
2478 127- C++20 Three-way comparison 1 / Partial Order, Total Order, epsilon, floating point failure
2479 https://www.youtube.com/watch?v=paA--t_fbOg&list=PL1_C6uWTeBDEjwkxjppnQ0TmMS272eNRx&index=132
2480
2481 044- Floating Point Error, C++23's New Features, expected, unexpected, denormalized, subnormal
2482 https://www.youtube.com/watch?v=wkcZr5TUmjU&list=PL1_C6uWTeBDEjwkxjppnQ0TmMS272eNRx&index=44
2483 */
2484 if(diff <= tolerance)
2485 return s2;
2486 else
2487 return adaptive_simpson_quadrature(f, a, c)
2488 + adaptive_simpson_quadrature(f, c, b); // [a, c] + [c, b]
2489 }
2490 }
2491
2492 // end of namespace numerical_analysis
2493}
2494// end of namespace cpg
2495
2496#endif
count
std::atomic< int > count
Definition: 022-mutex.cpp:10
function
void function(Type arg)
Definition: 041-tight_control.cpp:7
parallel_for
bool parallel_for(CallbackType &&callback, PolicyType &&policy, BeginType begin_index, EndType end_index)
Definition: 045-parallel_for.cpp:17
cpg_types.hpp
cpg::bitwise::common_t
std::common_type_t< S, T > common_t
Definition: cpg_bitwise.hpp:79
cpg::numerical_analysis::commands::d_6_6
constexpr auto d_6_6
Definition: cpg_calculus.hpp:1027
cpg::numerical_analysis::commands::dx_dy
constexpr auto dx_dy
Definition: cpg_calculus.hpp:1063
cpg::numerical_analysis::commands::d_8_2
constexpr auto d_8_2
Definition: cpg_calculus.hpp:1043
cpg::numerical_analysis::commands::d_3_1
constexpr auto d_3_1
Definition: cpg_calculus.hpp:992
cpg::numerical_analysis::commands::d_8_9
constexpr auto d_8_9
Definition: cpg_calculus.hpp:1050
cpg::numerical_analysis::commands::d_4_3
constexpr auto d_4_3
Definition: cpg_calculus.hpp:1004
cpg::numerical_analysis::commands::dx_2
constexpr auto dx_2
Definition: cpg_calculus.hpp:933
cpg::numerical_analysis::commands::d_1_6
constexpr auto d_1_6
Definition: cpg_calculus.hpp:977
cpg::numerical_analysis::commands::dz_6
constexpr auto dz_6
Definition: cpg_calculus.hpp:957
cpg::numerical_analysis::commands::d_2_3
constexpr auto d_2_3
Definition: cpg_calculus.hpp:984
cpg::numerical_analysis::commands::dx_dz_dy
constexpr auto dx_dz_dy
Definition: cpg_calculus.hpp:1067
cpg::numerical_analysis::commands::dx_dy_dz
constexpr auto dx_dy_dz
Definition: cpg_calculus.hpp:1066
cpg::numerical_analysis::commands::dz_5
constexpr auto dz_5
Definition: cpg_calculus.hpp:956
cpg::numerical_analysis::commands::d_6_9
constexpr auto d_6_9
Definition: cpg_calculus.hpp:1030
cpg::numerical_analysis::commands::d_3_7
constexpr auto d_3_7
Definition: cpg_calculus.hpp:998
cpg::numerical_analysis::commands::d_1_5
constexpr auto d_1_5
Definition: cpg_calculus.hpp:976
cpg::numerical_analysis::commands::d_9_8
constexpr auto d_9_8
Definition: cpg_calculus.hpp:1059
cpg::numerical_analysis::commands::d_6_8
constexpr auto d_6_8
Definition: cpg_calculus.hpp:1029
cpg::numerical_analysis::commands::d_6_7
constexpr auto d_6_7
Definition: cpg_calculus.hpp:1028
cpg::numerical_analysis::commands::d_8_7
constexpr auto d_8_7
Definition: cpg_calculus.hpp:1048
cpg::numerical_analysis::commands::d_4_8
constexpr auto d_4_8
Definition: cpg_calculus.hpp:1009
cpg::numerical_analysis::commands::dy_9
constexpr auto dy_9
Definition: cpg_calculus.hpp:950
cpg::numerical_analysis::commands::d_9_1
constexpr auto d_9_1
Definition: cpg_calculus.hpp:1052
cpg::numerical_analysis::commands::d_5_5
constexpr auto d_5_5
Definition: cpg_calculus.hpp:1016
cpg::numerical_analysis::commands::d_9_5
constexpr auto d_9_5
Definition: cpg_calculus.hpp:1056
cpg::numerical_analysis::commands::d_1_4
constexpr auto d_1_4
Definition: cpg_calculus.hpp:975
cpg::numerical_analysis::commands::d_9_3
constexpr auto d_9_3
Definition: cpg_calculus.hpp:1054
cpg::numerical_analysis::commands::d_5_9
constexpr auto d_5_9
Definition: cpg_calculus.hpp:1020
cpg::numerical_analysis::commands::d_7_3
constexpr auto d_7_3
Definition: cpg_calculus.hpp:1034
cpg::numerical_analysis::commands::d_2_6
constexpr auto d_2_6
Definition: cpg_calculus.hpp:987
cpg::numerical_analysis::commands::d_0_2
constexpr auto d_0_2
Definition: cpg_calculus.hpp:963
cpg::numerical_analysis::commands::iv_y
constexpr auto iv_y
Definition: cpg_calculus.hpp:929
cpg::numerical_analysis::commands::d_3_4
constexpr auto d_3_4
Definition: cpg_calculus.hpp:995
cpg::numerical_analysis::commands::d_4_2
constexpr auto d_4_2
Definition: cpg_calculus.hpp:1003
cpg::numerical_analysis::commands::d_6_5
constexpr auto d_6_5
Definition: cpg_calculus.hpp:1026
cpg::numerical_analysis::commands::d_1_9
constexpr auto d_1_9
Definition: cpg_calculus.hpp:980
cpg::numerical_analysis::commands::d_9_7
constexpr auto d_9_7
Definition: cpg_calculus.hpp:1058
cpg::numerical_analysis::commands::d_2_1
constexpr auto d_2_1
Definition: cpg_calculus.hpp:982
cpg::numerical_analysis::commands::d_7_8
constexpr auto d_7_8
Definition: cpg_calculus.hpp:1039
cpg::numerical_analysis::commands::d_8_5
constexpr auto d_8_5
Definition: cpg_calculus.hpp:1046
cpg::numerical_analysis::commands::dx_9
constexpr auto dx_9
Definition: cpg_calculus.hpp:940
cpg::numerical_analysis::commands::iv_x
constexpr auto iv_x
Definition: cpg_calculus.hpp:928
cpg::numerical_analysis::commands::dx_1
constexpr auto dx_1
Definition: cpg_calculus.hpp:932
cpg::numerical_analysis::commands::d_0_8
constexpr auto d_0_8
Definition: cpg_calculus.hpp:969
cpg::numerical_analysis::commands::d_3_3
constexpr auto d_3_3
Definition: cpg_calculus.hpp:994
cpg::numerical_analysis::commands::d_6_1
constexpr auto d_6_1
Definition: cpg_calculus.hpp:1022
cpg::numerical_analysis::commands::iv_z
constexpr auto iv_z
Definition: cpg_calculus.hpp:930
cpg::numerical_analysis::commands::dx_4
constexpr auto dx_4
Definition: cpg_calculus.hpp:935
cpg::numerical_analysis::commands::d_5_4
constexpr auto d_5_4
Definition: cpg_calculus.hpp:1015
cpg::numerical_analysis::commands::dy_2
constexpr auto dy_2
Definition: cpg_calculus.hpp:943
cpg::numerical_analysis::commands::d_2_9
constexpr auto d_2_9
Definition: cpg_calculus.hpp:990
cpg::numerical_analysis::commands::d_5_6
constexpr auto d_5_6
Definition: cpg_calculus.hpp:1017
cpg::numerical_analysis::commands::dz_3
constexpr auto dz_3
Definition: cpg_calculus.hpp:954
cpg::numerical_analysis::commands::d_7_4
constexpr auto d_7_4
Definition: cpg_calculus.hpp:1035
cpg::numerical_analysis::commands::dz_7
constexpr auto dz_7
Definition: cpg_calculus.hpp:958
cpg::numerical_analysis::commands::d_2_8
constexpr auto d_2_8
Definition: cpg_calculus.hpp:989
cpg::numerical_analysis::commands::d_5_1
constexpr auto d_5_1
Definition: cpg_calculus.hpp:1012
cpg::numerical_analysis::commands::dx_8
constexpr auto dx_8
Definition: cpg_calculus.hpp:939
cpg::numerical_analysis::commands::dy_dz_dx
constexpr auto dy_dz_dx
Definition: cpg_calculus.hpp:1070
cpg::numerical_analysis::commands::dz_1
constexpr auto dz_1
Definition: cpg_calculus.hpp:952
cpg::numerical_analysis::commands::d_1_3
constexpr auto d_1_3
Definition: cpg_calculus.hpp:974
cpg::numerical_analysis::commands::d_7_2
constexpr auto d_7_2
Definition: cpg_calculus.hpp:1033
cpg::numerical_analysis::commands::dz_8
constexpr auto dz_8
Definition: cpg_calculus.hpp:959
cpg::numerical_analysis::commands::d_4_6
constexpr auto d_4_6
Definition: cpg_calculus.hpp:1007
cpg::numerical_analysis::commands::dx_3
constexpr auto dx_3
Definition: cpg_calculus.hpp:934
cpg::numerical_analysis::commands::d_7_1
constexpr auto d_7_1
Definition: cpg_calculus.hpp:1032
cpg::numerical_analysis::commands::d_7_5
constexpr auto d_7_5
Definition: cpg_calculus.hpp:1036
cpg::numerical_analysis::commands::d_1_8
constexpr auto d_1_8
Definition: cpg_calculus.hpp:979
cpg::numerical_analysis::commands::dz_dx_dy
constexpr auto dz_dx_dy
Definition: cpg_calculus.hpp:1072
cpg::numerical_analysis::commands::d_0_1
constexpr auto d_0_1
Definition: cpg_calculus.hpp:962
cpg::numerical_analysis::commands::d_3_2
constexpr auto d_3_2
Definition: cpg_calculus.hpp:993
cpg::numerical_analysis::commands::d_6_3
constexpr auto d_6_3
Definition: cpg_calculus.hpp:1024
cpg::numerical_analysis::commands::dy_1
constexpr auto dy_1
Definition: cpg_calculus.hpp:942
cpg::numerical_analysis::commands::dz_dy_dx
constexpr auto dz_dy_dx
Definition: cpg_calculus.hpp:1073
cpg::numerical_analysis::commands::d_8_3
constexpr auto d_8_3
Definition: cpg_calculus.hpp:1044
cpg::numerical_analysis::commands::d_4_1
constexpr auto d_4_1
Definition: cpg_calculus.hpp:1002
cpg::numerical_analysis::commands::d_2_7
constexpr auto d_2_7
Definition: cpg_calculus.hpp:988
cpg::numerical_analysis::commands::dx_5
constexpr auto dx_5
Definition: cpg_calculus.hpp:936
cpg::numerical_analysis::commands::d_4_7
constexpr auto d_4_7
Definition: cpg_calculus.hpp:1008
cpg::numerical_analysis::commands::d_3_9
constexpr auto d_3_9
Definition: cpg_calculus.hpp:1000
cpg::numerical_analysis::commands::d_0_5
constexpr auto d_0_5
Definition: cpg_calculus.hpp:966
cpg::numerical_analysis::commands::d_7_9
constexpr auto d_7_9
Definition: cpg_calculus.hpp:1040
cpg::numerical_analysis::commands::d_2_2
constexpr auto d_2_2
Definition: cpg_calculus.hpp:983
cpg::numerical_analysis::commands::d_1_2
constexpr auto d_1_2
Definition: cpg_calculus.hpp:973
cpg::numerical_analysis::commands::dy_6
constexpr auto dy_6
Definition: cpg_calculus.hpp:947
cpg::numerical_analysis::commands::d_0_6
constexpr auto d_0_6
Definition: cpg_calculus.hpp:967
cpg::numerical_analysis::commands::dx_6
constexpr auto dx_6
Definition: cpg_calculus.hpp:937
cpg::numerical_analysis::commands::d_7_6
constexpr auto d_7_6
Definition: cpg_calculus.hpp:1037
cpg::numerical_analysis::commands::dz_2
constexpr auto dz_2
Definition: cpg_calculus.hpp:953
cpg::numerical_analysis::commands::d_9_9
constexpr auto d_9_9
Definition: cpg_calculus.hpp:1060
cpg::numerical_analysis::commands::d_8_6
constexpr auto d_8_6
Definition: cpg_calculus.hpp:1047
cpg::numerical_analysis::commands::dy_3
constexpr auto dy_3
Definition: cpg_calculus.hpp:944
cpg::numerical_analysis::commands::d_5_3
constexpr auto d_5_3
Definition: cpg_calculus.hpp:1014
cpg::numerical_analysis::commands::d_5_7
constexpr auto d_5_7
Definition: cpg_calculus.hpp:1018
cpg::numerical_analysis::commands::d_5_8
constexpr auto d_5_8
Definition: cpg_calculus.hpp:1019
cpg::numerical_analysis::commands::d_9_4
constexpr auto d_9_4
Definition: cpg_calculus.hpp:1055
cpg::numerical_analysis::commands::d_7_7
constexpr auto d_7_7
Definition: cpg_calculus.hpp:1038
cpg::numerical_analysis::commands::d_3_5
constexpr auto d_3_5
Definition: cpg_calculus.hpp:996
cpg::numerical_analysis::commands::d_4_5
constexpr auto d_4_5
Definition: cpg_calculus.hpp:1006
cpg::numerical_analysis::commands::d_8_4
constexpr auto d_8_4
Definition: cpg_calculus.hpp:1045
cpg::numerical_analysis::commands::dy_dx
constexpr auto dy_dx
Definition: cpg_calculus.hpp:1064
cpg::numerical_analysis::commands::d_5_2
constexpr auto d_5_2
Definition: cpg_calculus.hpp:1013
cpg::numerical_analysis::commands::d_0_3
constexpr auto d_0_3
Definition: cpg_calculus.hpp:964
cpg::numerical_analysis::commands::d_8_8
constexpr auto d_8_8
Definition: cpg_calculus.hpp:1049
cpg::numerical_analysis::commands::dy_dx_dz
constexpr auto dy_dx_dz
Definition: cpg_calculus.hpp:1069
cpg::numerical_analysis::commands::d_1_7
constexpr auto d_1_7
Definition: cpg_calculus.hpp:978
cpg::numerical_analysis::commands::dy_5
constexpr auto dy_5
Definition: cpg_calculus.hpp:946
cpg::numerical_analysis::commands::d_0_9
constexpr auto d_0_9
Definition: cpg_calculus.hpp:970
cpg::numerical_analysis::commands::d_8_1
constexpr auto d_8_1
Definition: cpg_calculus.hpp:1042
cpg::numerical_analysis::commands::d_0_7
constexpr auto d_0_7
Definition: cpg_calculus.hpp:968
cpg::numerical_analysis::commands::d_4_4
constexpr auto d_4_4
Definition: cpg_calculus.hpp:1005
cpg::numerical_analysis::commands::d_4_9
constexpr auto d_4_9
Definition: cpg_calculus.hpp:1010
cpg::numerical_analysis::commands::dy_4
constexpr auto dy_4
Definition: cpg_calculus.hpp:945
cpg::numerical_analysis::commands::dy_8
constexpr auto dy_8
Definition: cpg_calculus.hpp:949
cpg::numerical_analysis::commands::dy_7
constexpr auto dy_7
Definition: cpg_calculus.hpp:948
cpg::numerical_analysis::commands::dz_4
constexpr auto dz_4
Definition: cpg_calculus.hpp:955
cpg::numerical_analysis::commands::d_9_2
constexpr auto d_9_2
Definition: cpg_calculus.hpp:1053
cpg::numerical_analysis::commands::d_3_8
constexpr auto d_3_8
Definition: cpg_calculus.hpp:999
cpg::numerical_analysis::commands::dz_9
constexpr auto dz_9
Definition: cpg_calculus.hpp:960
cpg::numerical_analysis::commands::d_6_4
constexpr auto d_6_4
Definition: cpg_calculus.hpp:1025
cpg::numerical_analysis::commands::d_0_4
constexpr auto d_0_4
Definition: cpg_calculus.hpp:965
cpg::numerical_analysis::commands::d_3_6
constexpr auto d_3_6
Definition: cpg_calculus.hpp:997
cpg::numerical_analysis::commands::d_6_2
constexpr auto d_6_2
Definition: cpg_calculus.hpp:1023
cpg::numerical_analysis::commands::d_9_6
constexpr auto d_9_6
Definition: cpg_calculus.hpp:1057
cpg::numerical_analysis::commands::d_2_4
constexpr auto d_2_4
Definition: cpg_calculus.hpp:985
cpg::numerical_analysis::commands::d_2_5
constexpr auto d_2_5
Definition: cpg_calculus.hpp:986
cpg::numerical_analysis::commands::dx_7
constexpr auto dx_7
Definition: cpg_calculus.hpp:938
cpg::numerical_analysis::commands::d_1_1
constexpr auto d_1_1
Definition: cpg_calculus.hpp:972
cpg::numerical_analysis::curls
auto curls(FuncTypeX &&func_x, FuncTypeY &&func_y, FuncTypeZ &&func_z, std::array< BoundType, 2 > bound_x, std::array< BoundType, 2 > bound_y, std::array< BoundType, 2 > bound_z)
Definition: cpg_calculus.hpp:2300
cpg::numerical_analysis::divs
auto divs(FuncTypeX &&func_x, FuncTypeY &&func_y, FuncTypeZ &&func_z, std::array< BoundType, 2 > bound_x, std::array< BoundType, 2 > bound_y, std::array< BoundType, 2 > bound_z)
Definition: cpg_calculus.hpp:2370
cpg::numerical_analysis::smart_apply
auto smart_apply(FuncType &&func, TupleType arg)
Definition: cpg_calculus.hpp:604
cpg::numerical_analysis::differentiate
auto differentiate(cpt::sequence< VarIndex, DerivativeOrder >, FuncType &&func, std::size_t N, std::array< BoundType, 2 > bound)
Definition: cpg_calculus.hpp:2150
cpg::numerical_analysis::nine_point_stencil
auto nine_point_stencil(FuncType &&f, ArgType x) noexcept
Definition: cpg_calculus.hpp:1082
cpg::numerical_analysis::weights_abscissae
constexpr ReturnType weights_abscissae() noexcept
Definition: cpg_calculus.hpp:79
cpg::numerical_analysis::parametric_derivative
auto parametric_derivative(cpt::type_container< CmdTypes... > cmd, FuncType &&f, std::tuple< VariableTypes... > const &vars, ParamTypes... ps)
Definition: cpg_calculus.hpp:1500
cpg::numerical_analysis::negative_zero
constexpr Type negative_zero
Definition: cpg_calculus.hpp:63
cpg::numerical_analysis::seven_point_stencil
auto seven_point_stencil(FuncType &&f, ArgType x) noexcept
Definition: cpg_calculus.hpp:1198
cpg::numerical_analysis::compute_delta
BoundType compute_delta(CountType count, std::array< BoundType, 2 > &bound) noexcept
Definition: cpg_calculus.hpp:874
cpg::numerical_analysis::create_command
constexpr auto create_command(SeqType, SeqTypes...) noexcept
Definition: cpg_calculus.hpp:919
cpg::numerical_analysis::adjust_zero
Type adjust_zero(auto arg) noexcept
Definition: cpg_calculus.hpp:898
cpg::numerical_analysis::integral
BoundType integral(FuncType &&f, std::tuple< BoundType, BoundType > bound)
Definition: cpg_calculus.hpp:700
cpg::numerical_analysis::negative_approach
constexpr Type negative_approach(Type value)
Definition: cpg_calculus.hpp:54
cpg::numerical_analysis::gaussian_quadrature
ValueType gaussian_quadrature(FunctionType &&f, ValueType x1, ValueType x2)
Definition: cpg_calculus.hpp:552
cpg::numerical_analysis::fix_variables_other_than_VarIndex_ed
auto fix_variables_other_than_VarIndex_ed(FuncType &&func, ArgTypes... args) noexcept
Definition: cpg_calculus.hpp:1341
cpg::numerical_analysis::gradients
auto gradients(FuncType &&func, std::array< BoundType, 2 > bound_x, std::array< BoundType, 2 > bound_y, std::array< BoundType, 2 > bound_z)
Definition: cpg_calculus.hpp:2234
cpg::numerical_analysis::adaptive_gaussian_quadrature
BoundType adaptive_gaussian_quadrature(FuncType &&f, BoundType a, BoundType b)
Definition: cpg_calculus.hpp:572
cpg::numerical_analysis::gradient
auto gradient(FuncType &&function, ArgTypes... args) noexcept
Definition: cpg_calculus.hpp:1563
cpg::numerical_analysis::simpson_rule
BoundType simpson_rule(FuncType &&f, BoundType a, BoundType b) noexcept
Definition: cpg_calculus.hpp:2436
cpg::numerical_analysis::divergence
auto divergence(FuncTypeX &&func_x, FuncTypeY &&func_y, FuncTypeZ &&func_z, ArgTypes... args) noexcept
Definition: cpg_calculus.hpp:1672
cpg::numerical_analysis::positive_approach
constexpr Type positive_approach(Type value)
Definition: cpg_calculus.hpp:47
cpg::numerical_analysis::curl
auto curl(FuncTypeX &&func_x, FuncTypeY &&func_y, FuncTypeZ &&func_z, const std::array< ArgType, N > &args) noexcept
Definition: cpg_calculus.hpp:1654
cpg::numerical_analysis::adaptive_simpson_quadrature
BoundType adaptive_simpson_quadrature(FuncType &&f, BoundType a, BoundType b)
Definition: cpg_calculus.hpp:2465
cpg::numerical_analysis::positive_infinity
constexpr Type positive_infinity(Type v)
Definition: cpg_calculus.hpp:34
cpg::numerical_analysis::negative_infinity
constexpr Type negative_infinity(Type v)
Definition: cpg_calculus.hpp:40
cpg::numerical_analysis::partial_derivative
auto partial_derivative(cpt::sequence< VarIndex, Order >, FuncType &&func, ArgTypes... args) noexcept
Definition: cpg_calculus.hpp:1376
cpg::numerical_analysis::curl
auto curl(FuncTypeX &&func_x, FuncTypeY &&func_y, FuncTypeZ &&func_z, ArgTypes... args) noexcept
Definition: cpg_calculus.hpp:1613
cpg::numerical_analysis::five_point_stencil
ArgType five_point_stencil(FunctionType &&f, ArgType x) noexcept
Definition: cpg_calculus.hpp:1282
cpg::numerical_analysis::evaluate
auto evaluate(FuncType &&f, ArgType arg1)
Definition: cpg_calculus.hpp:614
cpg::numerical_analysis::directional_derivative
auto directional_derivative(FuncType &&f, std::tuple< VariableTypes... > const &vars, ParamType param)
Definition: cpg_calculus.hpp:1480
cpg::numerical_analysis::divergence
auto divergence(FuncTypeX &&func_x, FuncTypeY &&func_y, FuncTypeZ &&func_z, const std::array< ArgType, N > &args) noexcept
Definition: cpg_calculus.hpp:1708
cpg::numerical_analysis::adjust_integer
Type adjust_integer(auto arg) noexcept
Definition: cpg_calculus.hpp:891
cpg::numerical_analysis::positive_zero
constexpr Type positive_zero
Definition: cpg_calculus.hpp:60
cpg::numerical_analysis::gradient
auto gradient(FuncType &&function, const std::array< ArgType, N > &args) noexcept
Definition: cpg_calculus.hpp:1596
cpg::numerical_analysis::get_delta
constexpr auto get_delta(delta< DeltaType > del) noexcept
Definition: cpg_calculus.hpp:849
cpg::types::cast_ref
auto & cast_ref(std::unique_ptr< T[], Deleter > &uptr) noexcept
Definition: cpg_types.hpp:1290
cpg::types::sequence
std::integer_sequence< std::common_type_t< std::remove_cvref_t< decltype(Indices)>... >, Indices... > sequence
Definition: cpg_types.hpp:2665
cpg
Includes subnamespace conversion.
Definition: cpg_3dvector.hpp:11
cpg::endl
constexpr auto endl
Definition: cpg_types.hpp:213
stl_extensions::apply
constexpr decltype(auto) apply(F &&f, T(&&c_array)[N])
tpf::split_range::range
std::pair< T, T > range
Definition: tpf_split_ranges.hpp:41
tpf::types::hidden::common_type_t
typename st_common_type_solver< Types... >::type common_type_t
Definition: tpf_types.hpp:4119
tpf::types::remove_cvref_t
remove_cv_ref_t< Type > remove_cvref_t
Definition: tpf_types.hpp:298
tpf::types::is_invocable_v
constexpr bool is_invocable_v
Definition: tpf_types.hpp:6063
cpg::numerical_analysis::delta
Definition: cpg_calculus.hpp:837
cpg::types::type_container
Definition: cpg_types.hpp:827